© The Authors 2025

1 Introduction

Runaway stars move with a high peculiar velocity relative to their environment (see Blaauw 1961 for a seminal work, and Carretero-Castrillo et al. 2023 for an updated reference). Aside from the studies of the origin of their runaway status, there is a whole area of research related to the phenomena that occur when a runaway star, through its strong stellar wind, interacts with the interstellar medium (ISM).

Stars of spectral types O and early B are among those with the most intense stellar winds. About 30% of O stars and 5–10% of B stars have been found to be runaway stars (Blaauw 1961; Stone 1979; Carretero-Castrillo et al. 2023). Detailed studies to identify runaway stars have been performed by several authors (Tetzlaff et al. 2011; Boubert & Evans 2018; Maíz Apellániz et al. 2018; Kobulnicky & Chick 2022; Carretero-Castrillo et al. 2023). In this context, a stellar bow shock is formed by the accumulation of ISM matter as a runaway star with a powerful wind passes through the ISM (first findings, e.g. Gull & Sofia 1979). If the motion is supersonic (i.e., the peculiar velocity of the star is greater than the speed of sound in the environment that is typically of some kilometers per second) strong shocks take place around a discontinuity surface separating the two media (see for instance Mackey 2023, and references therein).

The swept-up matter, formed by warm dust and heated by stellar ultraviolet photons, shines at infrared (IR) wavelengths. The bow shocks usually show an arc-shaped structure but sometimes even a bubble surrounding the star, where part of the edge can be enhanced (that in the direction of the stellar motion). Initial searches for nearby runaway stars were carried out by Noriega-Crespo et al. (1997) (see also references therein). The authors studied high-resolution IRAS satellite images around runaway OB stars and published the first catalog of 21 possible bow shocks and 4 possible bubbles. One of them, later named EB27 and associated with the star BD+43 3654, was well resolved with IR data from the MSX satellite by Comerón & Pasquali (2007). More recently, four relevant catalogs of stellar bow shocks with IR emission have been compiled. The first one compiled by Peri et al. (2012, 2015), called the Extensive stellar Bow Shock Survey (EBOSS), presented 73 object candidates in WISE images. Almost half of them had an associated runaway star. For the rest, there was either no information on the velocity of the star (i.e., its runaway status) that could be associated with the bow shock, or, in the case of a few of them, no associated star. The second catalog was compiled by Kobulnicky et al. (2016), with 709 candidates detected using WISE and Spitzer images. In Kobulnicky et al. (2017), the team publish the photometric data, and in Chick et al. (2020), they study the 84 – out of 709 – cases in which stars have been identified supporting the bowshock nebulae. The third one, Bodensteiner et al. (2018), was compiled from different O, B, Be, and A-type star catalogs, but for general IR nebulae (including bow shocks) around massive stars. In this work, they test if the IR nebulae information could be a binary indicator, and confirm that for about 29% of their sample. The fourth one, Jayasinghe et al. (2019), contains 1394 bubbles and 453 bow shocks, compiled through the Milky Way Project using Spitzer data.

In terms of theoretical developments, the shape, dynamics, and emission of bow shocks have been modeled. Wilkin (1996) presented one of the first solution methods that provided, among other parameters, the shape of the shell, the mass column density, and the velocity distribution of the shocked gas within it. The spectral energy distribution (SED) down to gamma rays of stellar bow shocks formed by O-type runaway stars was computed by del Valle & Romero (2012) using ρ Oph as a test case. del Palacio et al. (2018) developed a multi-zone stellar bow shock model and SED for different scenarios, along with synthetic radio maps to reproduce the morphology of EB27 and predicted the high-energy (HE) emission. Using deep Very Large Array (VLA) low-frequency observations of the same target, Benaglia et al. (2021) model the intensity and morphology of the radio emission.

Strong bow shocks are promising places for the acceleration of particles up to relativistic velocities. Relativistic particles radiate across the electromagnetic spectrum, particularly in the radio, X-rays, and gamma rays (del Palacio et al. 2018; Martinez et al. 2023). In this context, Binder et al. (2019) searched for X- ray emission after stacking archival Chandra observations from 60 IR-bright Galactic bow shocks. Although there was no detection, the authors could establish a constraining upper limit on the X-ray luminosity of IR-detected bow shocks of <2 × 10²⁹ erg s⁻¹, which is on the same order as model predictions. The question of whether particles can be accelerated up to relativistic energies at the shocks of these structures was investigated in the radio domain by Benaglia et al. (2010). They performed L- and C-band radio observations of the source EB27 with the VLA, to look for nonthermal radio emission. Their ultimate goal was to determine whether bow shocks could be the counterparts of HE gammaray sources, especially of the unassociated ones. Radio emission coinciding with the IR extension of the EB27 bow shock was detected in both bands, and a central part of it showed a spectral index, α (with S_ν ~ ν^α), close to −0.5, which is characteristic of nonthermal emission. Follow-up JVLA observations in the S-band in full polarization mode showed no polarized emission above the 1% degree level (Benaglia et al. 2021). Martinez et al. (2023) have studied the nature of the radio emission for EB27 using a multi-zone model. Compared with radio observations, they conclude that the emission of EB27 can be explained through nonthermal electrons following a hard energy distribution below ~1 GeV. Nevertheless, the emission from EB27 together with its spectral index and polarization absence needs further investigation.

In recent years, the commissioning of interferometers in the southern hemisphere has allowed several radio detections of stellar bow shocks and spectral index studies. van den Eijnden et al. (2022a) detected the Vela X-1 bow shock in radio with MeerKAT and found that its multiwavelength emission is more compatible with thermal free-free emission than with nonthermal synchrotron emission. van den Eijnden et al. (2022b) examined the Rapid ASKAP Continuum Survey (RACS, McConnell et al. 2020) at the positions of 50 EBOSS candidates and detected 10 structures, classifying them as 3 radio bow shocks, 3 probable radio bow shocks, 1 to be confirmed, and 3 unclear (probably not radio bow shocks). The analysis of the first seven structures yielded results consistent with nonthermal emission for two of them. In addition, Moutzouri et al. (2022) measured a spectral index below −0.5 for the northern of the Bubble Nebula (NGC 7635, most likely formed by the runaway star BD+60 2522), which is the second bow shock with clear evidence for nonthermal emission. However, the question of bow shocks capable of accelerating relativistic particles and producing HE emission remains open.

The runaway stars associated with the bow shocks may also produce radio emission. In binary systems containing a massive O or Be star and a compact object, nonthermal radio emission can be produced by synchrotron-emitting relativistic electrons. There are two main types of these binary systems: (a) X-ray binaries with a radio jet or microquasars (Mirabel & Rodríguez 1999; Fender & Muñoz-Darias 2016); and (b) gammaray binaries (Dubus 2006; Bosch-Ramon et al. 2012; Dubus 2013). For the latter, three of the nine confirmed TeV-emitting gamma-ray binaries are runaways: LS 5039, PSR B1259-63 and 1FGL J1018.6-5856 (Ribó et al. 2002; Miller-Jones et al. 2018; Marcote et al. 2018). Therefore, new gamma-ray binaries could be found by searching for massive runaway stars and their associated multiwavelength emission. This strategy could provide new members of this reduced sample, eventually allowing for the population studies that are needed to answer some of the many open questions in this area of HE astrophysics (see Bordas 2023, and references therein).

For all of the above reasons, any study of stellar bow shocks and their radio counterparts can add information on a wide variety of phenomena. This is especially true if the runaway star that produces the bow shock is identified. In this work, we aim to detect bow shocks around new runaway stars and describe their characteristics. In Sect. 2, we present the data used in this work. In Sect. 3, we describe the methodology of the search, and the identification and characterization of the objects. In Sect. 4, we present the results for the different types of objects we found. We discuss these results in Sect. 5. We close with a summary and prospects for further studies in Sect. 6.

2 Data

2.1 Runaway star catalogs

Bow shocks are expected around massive runaway stars. Therefore, we used the catalogs of runaway O and Be stars obtained by Carretero-Castrillo et al. (2023), which were identified using Gaia DR3 astrometric data. In this work, the authors crossmatched the GOSC (Maíz Apellániz et al. 2013) and the BeSS (Neiner et al. 2011) catalogs with Gaia DR3 data. After applying several quality cuts, they produced the GOSC-Gaia DR3 catalog with 417 O-type stars, and the BeSS-Gaia DR3 catalog with 1335 Be-type stars. In the GOSC-Gaia DR3 catalog, they found 106 O runaway stars, with 42 new identifications as runaways. In the BeSS-Gaia DR3 catalog they found 69 Be runaway stars, with 47 new identifications as runaways. The runaway percentages found with respect to the catalogs are 25.4% and 5.2% for the O and the Be runaway stars, respectively. The runaway stars were identified through an E value (see Eq. (3) in Carretero-Castrillo et al. 2023), which indicates the normalized significance of the runaway nature. Therefore, higher E values denote more confident runaway star identifications. The first ten rows of data corresponding to the runaways with the higher E values, of the GOSC and BeSS-Gaia DR3 runaway catalogs, respectively, are shown in Tables 1 and 2. These tables include the runaway id from the original work, the star names, the Gaia DR3 ids, coordinates (J2000), distances, corrected proper motions (see Sect. 3), G magnitudes, spectral types, the twodimensional peculiar velocities computed in Carretero-Castrillo et al. (2023), E values, the runaway references (in case the star had already been reported as runaway), and a last column indicating the reference in the cases that they had already been reported to contain bow shocks. As for the runaway id, GR (GOSC Runaways) refers to runaways identified from the GOSC catalog, while BR (BeSS Runaways) to those from the BeSS catalog.

2.2 Infrared data: WISE

The Wide-field Infrared Survey Explorer (WISE, Wright et al. 2010) was a mission funded by NASA that mapped the whole sky at mid-IR wavelengths from December 2009 to February 2011. It had a much better sensitivity than that of earlier IR survey missions. The different survey bands and resolutions are shown in Table 3.

Bow shocks are related to the warm dust emission, and this emission is mostly expected in the W4 band. Therefore, we used the W4 band to look for bow shock or bubble-like structures around the runaway stars. WISE images are available through the WISE Image Service¹ of the NASA/IPAC Infrared Science Archive.

3 Methodology

3.1 Bow shock searches and identification

We performed and cross-checked a visual inspection of the W4 images for all the 106 and 69 runaways in Tables 1 and 2, respectively. We used the astronomical data visualization tool SAOImageDS9 (Joye & Mandel 2003) to process the WISE images. We used WISE images of size 20′ × 20′ to make sure that we did not lose part of the possible bow shock structures. We looked for bow shock or bubble-like structures around the position of the runaway stars. We also found some features that we would classify not as bow shocks or as bubble-like, but as miscellaneous structures near the runaway stars. These cases will be treated separately. The criterion defined to classify a bow shock or a bubble candidate around a runaway star was based on two aspects:

1. The presence of a bow shock or a bubble-like structure around the runaway star.

2. The arc-shaped structure, or (part of) the rim of the bubble, must be in the direction of the runaway star’s corrected proper motion.

To verify the second point, we used the Gaia DR3 proper motions corrected for the ISM motion caused by Galactic rotation. For this, we followed the prescription given in Comerόn & Pasquali (2007) (originally from Scheffler & Elsӓsser 1987), in which it is assumed that the ISM at any given position moves in a circular orbit around the Galactic center (GC) following a rotation curve. This prescription uses the Oort constants (see Oort 1927; Ogrodnikoff 1932), which can be computed for any given Galactic rotation curve. As was done in Carretero-Castrillo et al. (2023), we used the Galactic rotation curve provided by the A5 fit of Reid et al. (2019) and fit the original data to compute the rotation curve slope. The values we used for the velocity of the Sun, its distance to the GC, the circular rotation velocity at the position of the Sun, and its slope (for distances to the GC larger than 5 kpc) are: U_⊙ = 10.6 ±1.2 km s⁻¹, V_⊙ = 10.7 ± 6.0 km s⁻¹, W_⊙ = 7.6 ± 0.7 km s⁻¹, R_⊙ = 8.15 ± 0.15 kpc, Θ_⊙ = 236 ± 7 km s⁻¹, and dΘ/dR = −1.4 ± 0.3 km s⁻¹ kpc⁻¹. These are the same values used in all computations presented in Carretero-Castrillo et al. (2023), though some of them are incorrectly quoted in the original manuscript. We could then compute the Oort constants needed to compute the ISM proper motion, and found them to be A = 15.7 ± 0.5 km s⁻¹ and B = −13.8 ± 0.5 km s⁻¹. We finally subtracted the local ISM proper motions to the star’s proper motions, and found the proper motions of the star with respect to its local ISM (i.e., the corrected proper motions). The corrected proper motions associated with the runaway stars are listed in Tables 1 and 2. We also computed the propagation of uncertainties for these variables. From the corrected proper motions and their uncertainties, we computed for each star the position angle (PA) from north to east of the corrected proper motions and its uncertainty, ∆PA. For the O-runaway stars, ∆PA ranges from 0.7 to 15.5°, with a median of 3.9°. For the Be-runaway stars, ∆PA ranges from 1.1 to 12.4°, with a median of 5.1°. Thus, ∆PA is usually only a few degrees; it has a tendency to increase as the peculiar velocity of the star decreases, since the uncertainties of the corrected proper motions are larger compared to the modulus of the corrected proper motion itself. We included the corrected proper motions in α, in δ, and its modulus and PA, together with their corresponding uncertainties, in the electronic versions of Tables 1 and 2. Finally, to check for the second item above, we overlapped in the WISE images arrows with 2′ length in the direction of the corrected proper motion.

3.2 Bow shock characterization

Once the bow shock and bubble-like candidates were identified, we proceeded to their geometrical characterization. We measured the sizes of the IR geometrical structures in the W4 band of the bow shock and bubble candidates projected on the plane of the sky. For this purpose, we used the region tools provided by SAOImageDS9. For the bubble candidates, only the projected width, w – the diameter of the bubble structure – was measured. For the bow shock structures, we also measured the projected standoff distance, R, which is the distance between the bow shock apex and the star, and the projected length, l, of the bow shock. The shape of a bow shock structure can be approximated by a partial ellipse (see Fig. 1). Therefore, we used an ellipse shape to parameterize the bow shock in IR and measure its arc length. For each bow shock, we superimposed an ellipse to the W4 image and obtained its semimajor and semiminor axes. We determined the initial and final angles that define the arc length of the IR emitting region, measured from the center of the ellipse and defined counterclockwise from the semimajor axis. With all these parameters, we can compute the partial perimeter of the ellipse – that is, the length, l (continuous yellow arc in Fig. 1) – using the differential arc length equation for a curve defined by an ellipse. We note that these measurements were computed more accurately than the ones presented in Peri et al. (2012) and Peri et al. (2015). The results of the geometrical measures are presented in Sect. 4. The values given in linear units (parsec) were computed using the distances provided in Tables 1 and 2.

Fig. 1 Image of a known bow shock illustrating the geometric characterization used in this work. R is the projected standoff distance, l is the projected length of the bow shock, and w is its projected width. The arc of the ellipse shown with a continuous yellow line was used to derive l (see text for details).

3.3 Interstellar medium ambient density

We estimated the ISM ambient density, n_ISM, at the bow shock position using Eq. (1) of Wilkin (1996): $$R_0=\sqrt{\frac{\dot{M}{v}_{\infty }}{4\pi {\rho }_a{V}_{\text{PEC}}^2}},$$(1)

where R₀ is the standoff distance (not the projected one that we measure, R), Ṁ is the stellar mass-loss rate, υ_∞ is the wind terminal velocity, V_PEC is the peculiar velocity of the star, and ρ_a is the ambient medium density. From Eq. (1), we estimated ρ_a, and from ρ_a = μn_ISM, we finally derived n_ISM. In this equation, we assumed that R₀ ~ R, which will be the case if the motion occurs in the plane of the sky. We note that n_ISM is the volume density in H atoms, and µ is the mass of the ISM per H atom. We assumed that the fractional abundances for the hydrogen and the helium are X = 0.9, and Y = 0.1, respectively. Therefore, the mass of the ISM per H atom is then μ = μ_H + 0.1μ_He = [0.9 · 1.67 + 0.1 · 6.6] × 10⁻²⁴ g = 2.16 × 10⁻²⁴ g.

The stellar parameters, Ṁ and υ_∞, for each of our stars were derived from the tables in Vink & Sander (2021)², assuming solar metallicity for our runaway stars, and using the effective temperatures, masses, and luminosities from Martins et al. (2005). Runaway stars without a luminosity class (GR100 and BR65) were not considered given the large uncertainty in the determination of their temperatures, masses, and luminosities. Given the linear dependence of υ_∞ on the T_eff of the star, we computed a linear regression to the values quoted in Vink & Sander (2021), and used it for our stars. Owing to the variability of the mass-loss rates as a function of T_eff, we preferred to compute a linear interpolation to assign mass-loss rates to our stars.

For the values of V_PEC, we considered either the 2D or 3D peculiar velocities. The first ones are those from Carretero- Castrillo et al. (2023). To compute the 3D peculiar velocities of the runaway stars, we searched for their heliocentric radial velocities in the literature. For stars with measured values, we computed the regional standard of rest (RSR) velocities in the manner described in Carretero-Castrillo et al. (2023), but using real radial velocities instead of the theoretical ones used there. We then computed the 3D peculiar velocities, />, as the moduli of the RSR velocities. We highlight here that we used the 2D and 3D peculiar velocities to compute $n_{\text{ISM}}^{2\text{D}}$ and $n_{\text{ISM}}^{3\text{D}}$, respectively.

However, since the 3D peculiar velocities might be contaminated by radial velocity variations in the case of binary systems, $n_{\text{ISM}}^{2\text{D}}$ might be seen as more reliable, although contamination from binarity is unlikely to exceed a few tens of kilometers per second in the case of massive stars. We note that the obtained $n_{\text{ISM}}^{2\text{D}}$ values will always be higher than the $n_{\text{ISM}}^{3\text{D}}$ ones, as lower velocities $\left(V_{\text{PEC}}^{2\text{D}}\text{VS}.V_{\text{PEC}}^{3\text{D}}\right)$ require higher ISM densities to form a bow shock, assuming all other parameters remain the same in Eq. (1).

The stellar parameters, velocities, and computed ISM densities are shown in Tables 4 and 5. We caveat that the values of n_ISM should be taken as estimates, given the large uncertainties in the determinations of the stellar parameters, and the fact that we used the projected distance, R, as the R₀ value in Eq. (1).

4 Results

During the examination of the W4 images, we found a wide variety of IR structures around the runaway stars: bow shock and bubble-like candidates, intermediate structures³, and miscellaneous structures. There were also some bow shocks or bubbles that were already known. Therefore, we grouped the results in three different categories: new, known, and miscellaneous IR structures. We finish this section by summarizing the search for radio emission from the bow shocks and bubbles discussed in this work.

4.1 The new candidates for stellar bow shocks or bubbles

We identified nine new bow shock candidates, three bubble-like candidates, and one intermediate structure for the runaway stars in the GOSC- and BeSS-Gaia DR3 runaway catalogs. These 13 sources are presented in Table 4. Ten of these structures are associated with O-type runaway stars, and only three of them with Be-type runaway stars. Table 4 presents their associated runaway star, with stellar parameters and velocities, the type id (bow shock (BS) or bubble (BU)), and the geometrical measurements, R, l, and w, of the IR structures in the W4 images. We also provide the eccentricity, e, of the ellipse used to measure the length of the bow shocks. The last two columns present the estimated n_ISM 2D and 3D values for the bow shock candidates for which we have all the necessary data (see Eq. (1), where we used the projected value, R, as the real value, R₀). For the new bow shock candidates, the 2D ISM density, $n_{\text{ISM}}^{2\text{D}}$, ranges from 0.8 to ~220 cm⁻³, with a median of 6.4 cm⁻³. In the 3D case, $n_{\text{ISM}}^{3\text{D}}$ ranges from 0.6 to ~200 cm⁻³, with 3.8 cm⁻³ as the median⁴. As is mentioned in Sect. 3, $n_{\text{ISM}}^{2\text{D}}$ values are always higher than the $n_{\text{ISM}}^{3\text{D}}$ ones. BS-GR72 has the largest ISM densities probably due to the high mass-loss rate and wind terminal velocity of its runaway star, combined with a small value for the R distance.

For each candidate, we present the corresponding RGB image in W4+W3+W2 bands in Figs. 2 and A.1. The dashed white arrows indicate the directions of Gaia DR3 proper motions and the solid cyan arrows the directions of the corrected proper motions. We note that in some cases the corrected proper motion direction of the runaway star is not well aligned with the bow shock direction. We tabulate and comment on these differences in Appendix B. We note that some of these cases could be examples of in situ bow shocks (Povich et al. 2008; Kobulnicky et al. 2016), which result from the interaction of outflows from star-forming regions with the stellar wind of massive stars. The orientation of these bow shocks is the result of the flowing ISM and is thus not related to the runaway star motion. Thus, we further inspected the new bow shock fields that showed the largest misalignments, those with angular differences (ADs) greater than 5σ (see Appendix B), to determine whether they are in fact in situ bow shocks. For this, we searched for H II regions or molecular clouds in the Simbad database around the position of these bow shocks. We inspected larger WISE fields centered at the bow shock positions, and finally, we looked for counterparts in H II region catalogs such as the Green Bank Telescope H II Region Discovery Survey (Anderson et al. 2012) and the WISE Catalog of Galactic H II Regions (Anderson et al. 2014).

We shall comment now on the three cases with AD above 5σ: BS-GR71, BS-GR72, and BS-GR93. For BS-GR71, there were no H II regions found on the field. Therefore, this cannot be a case of an in situ bow shock and we have kept it as a bow shock candidate associated with the runaway star GR71. In contrast, BS-GR72 is coincident with the evolved H II region S206 (or NGC 1491), located at a distance of ${2.96}_{-0.15}^{+0.17}$ kpc (Méndez-Delgado et al. 2022). This distance is compatible with the one of the runaway star GR72, 2.94 ± 0.3 kpc (Carretero-Castrillo et al. 2023). Therefore, this bow shock could be indeed an in situ bow shock. In addition, we note that this H II region S206 presents radio emission (see Appendix C). BS-GR93 has the molecular cloud [DB2002b] G307.25+2.86 (Hartley et al. 1986) at ~10′, related to the globule BHR 87, whose distance appears to have different values in the literature: 400 pc (Launhardt et al. 2010) and 1 kpc (Bourke et al. 2005). However, none of these distances is compatible with the one of the runaway star GR93, 3.85 ± 0.35 kpc. Therefore, we cannot claim that BS-GR93 is an in situ bow shock candidate.

We shall now comment on two particularly interesting cases, BS-GR35, and BU-BR46. While working with the RGB image of BS-GR35, we discovered another bow shock candidate in the field, ~2′ south from the runaway GR35 (see Fig. 2). We called it BS J075339−2616.5 because of the coordinates of its center: RA = 07h 53m 39.1s and Dec = −26° 16′31″. It has a perfect bow structure in the W4 and W3 images. We show an RGB image centered and zoomed on this source in Fig. 3. We note that it cannot be associated with our runaway star GR35 because of its position and orientation with respect to the star and its proper motions. We looked for possible Gaia stars nearby BS J075339−2616.5 but did not find a clear candidate. We also looked for possible associated sources in the region, but there is only one IRAS source, IRAS 07515−2608. Notably, we found a radio detection for this new bow shock candidate in RACS-low and RACS-mid images, but it does not appear to be a clear counterpart (see Appendix C.2). As for BU-BR46, we originally proposed it as a bow shock candidate. However, after a detailed examination of the IR structure, we classified it as a bubble. We found that the bow shape was actually caused by another star inside the W4 region of emission, to the southwest of BU-BR46 (see Fig. A.1). It is a long-period variable star, 2MASS J19493059+3257068, that is not physically related to our runaway.

Fig. 2 WISE RGB images in W4+W3+W2 for three examples of the new bow shock and bubble candidates presented in Table 4 in equatorial coordinates, with north up and east to the left. Dashed white arrows indicate the directions of the original proper motions from Gaia DR3, and solid cyan arrows indicate the directions of the proper motions corrected for the ISM motion caused by Galactic rotation. Each field has a different size, but the arrows are fixed to 2′. Each panel contains the name of the star provided in the GOSC or the BeSS catalogs, and the runaway star id. The remaining RGB images of the sources from Table 4 are presented in Appendix A.

Fig. 3 WISE RGB image in W4+W3+W2 for the bow shock candidate BS J075339−2616.5 in equatorial coordinates, with north up and east to the left. The white cross indicates the center of the bow shock emission, at coordinates RA = 07h 53m 39.1s and Dec = −26°16′31″. This bow shock appeared in the WISE field of BS-GR35 (see Fig. 2), and has no associated runaway star.

4.2 Known candidates for stellar bow shocks or bubbles

There are a total of 17 known bow shocks or bubble-like candidates among the GOSC-Gaia DR3 runaway stars. There is only one known candidate among the BeSS-Gaia DR3 runaways, BS-BR14, which is indeed coincident with BS-GR42 because it is a runaway star common to the two catalogs⁵. Their bow shock or bubble references are from different catalog searches (Peri et al. 2012, 2015; Kobulnicky et al. 2016; Maíz Apellániz et al. 2018; Moutzouri et al. 2022), and are included in Tables 1 and 2. However, some of them have not been geometrically characterized in previous works. In this section, we focus on these cases. We selected those bow shocks and characterized them using the same procedure as for the new bow shocks discovered in this work. Table 5 presents their associated runaway stars, with stellar parameters and velocities, the type id, the geometrical measurements of the IR structures in W4, and the eccentricity, e. Again, we estimated the ISM densities for the known bow shock candidates. In this case, $n_{\text{ISM}}^{2\text{D}}$ ranges from 0.4 to ~200 cm⁻³, and $n_{\text{ISM}}^{3\text{D}}$ ranges from 0.3 to ~190 cm⁻³, with BS-GR83 presenting the largest value due to its short R distance.

As for the new candidates, we also present the corresponding WISE RGB images in W4+W3+W2 bands in Fig. 4. We shall comment here on two special cases: BS/BU-GR51 and EB43 of Peri et al. (2015), coincident with GR20.

BS/BU-GR51 is a known IR structure identified in Maíz Apellániz et al. (2018). In that work, they introduce it as an asymmetric region of hot dust around the runaway star rather than a bow shock. The likely existence of an IR bubble around GR51 was already identified in Simpson et al. (2012). However, given the asymmetric morphology, we have classified this candidate as an intermediate structure, since a bubble should be more spherical. The RGB image in W4+W3+W2 can be seen in Fig. 4.

We also comment on the particular case of EB43, a bow shock candidate of Peri et al. (2015). The position of this bow shock is coincident to the runaway GR20. However, the strange shape of EB43, cut in half just in the direction in which the bow shock is expected, does not meet the criterion presented in Sect. 3 and has made us discard it as a bow shock candidate. We note that it was already proposed as a doubtful candidate on the basis of morphology in Peri et al. (2015). Therefore, although the number of known bow shocks within our sample was initially 17, after discarding this one we considered only 16 cases.

4.3 Miscellaneous infrared structures

In this section, we present different types of IR structures found in the W4 images around the runaway stars. We cannot classify these various structures as bow shock or bubble candidates with certainty, but they present a particular morphology that makes them interesting in their own right. Therefore, we decided to group them in a category we called miscellaneous IR structures, and dedicate this section to them. All the sources belonging to this category are presented in Table 6, along with a morphological classification. In what follows, we comment on a particular subgroup that we call mini-bubbles and discuss the possible origin of the W4 emission.

When examining the W4 images, we found some runaway stars that exhibit relatively compact emission: 35 and 22 sources among the GOSC- and the BeSS-Gaia DR3 runaways, respectively. We note that all of them look point-like in the other WISE bands. We applied Guassian fits to the W4 images at the position of these sources. We found that most of them are compatible with being point-like. However, we found two particular Be-runaway stars, BR10 and BR15, that have excesses spanning around two times the WISE PSF in W4 (which is 12″). We then classified the objects in this subgroup as mini-bubbles, since they are not as extended and big as typical bubbles (the minimum size among our bubbles is 69″), but they are not compatible with point-like source emission⁶.

We looked for explanations for the W4 excesses found in point-like sources and mini-bubbles in the literature, but did not find a clear pattern. In the case of Be-type stars, this IR emission could be related to the presence of a circumstellar disk (Rivinius et al. 2013). We found a runaway star in the literature, V452 Sct, with an IR excess in W4 that is visually consistent with what we have called mini-bubbles. Its RGB image is presented in Fig. 6 of Maíz Apellániz et al. (2018). In this work, they present the IR emission as an excess in W4 consistent with it being a point source centered on the star, but not a bow shock. They also found that this excess was studied in Miroshnichenko et al. (2000), which proposed that it was likely caused by an LBV-type outburst in the past. We applied a Guassian fit to the W4 excess, and found that the size was larger than twice the WISE PSF in W4. Therefore, we also consider this runaway star to be a mini-bubble.

We also inspected the point-like sources and mini-bubbles with the VizieR Photometry viewer tool, which allows for easy visualization of photometry-enabled catalogs in VizieR. We used this tool to search for components in addition to the expected blackbody stellar component. For the case of V452 Sct, the minibubble of Maíz Apellániz et al. (2018), two components can be clearly distinguished. However, we only found another clear bump in BR29, and some hints in GR41, BR10, BR15, BR28, and BR34. Therefore, it is difficult to associate the mini-bubble nature with any additional component in their SEDs. We present in Appendix D the RGB images for the seven miscellaneous sources from Table 6 that are not point-like.

Fig. 4 WISE RGB images in W4+W3+W2 for the known bow shock and bubble candidates presented in Table 5 in equatorial coordinates, with north up and east to the left. Dashed white arrows indicate the directions of the original proper motions from Gaia DR3, and solid cyan arrows indicate the directions of the proper motions corrected for the ISM motion caused by Galactic rotation. Each field has a different size, but the arrows are fixed to 2′. Each panel contains the name of the star as provided in the GOSC or the BeSS catalogs, and the runaway star id.

4.4 Search for radio counterparts

We conducted a search for radio emission among the newly detected IR structures. Details of this search can be found in Appendix C, with data presented in Appendix C.1 and the obtained results shown in Appendix C.2. We did not find pointlike radio emission coincident with the positions of the runaway stars. We did not find any clear radio counterparts of the IR sources either, but we found radio emission partially overlapping them in six particular cases that are discussed in Appendix C.2: BS-GR9, BS-GR72, BS-GR75, BS J075339-2616.5, BS/BU- GR51, and BS-GR83.

We also investigated a thermal or nonthermal origin of the possible radio emission for the bow shock candidates, based on their IR characterization and the RACS radio data used in this work (see Appendix C.3). We conclude that the sources BS-GR75 and BS-GR83 should be detectable in radio for standard parameters. Future, more sensitive radio observations could reveal radio emission in these sources, which are more likely to show nonthermal than thermal radio emission (see details in Appendix C.3).

5 Discussion

In this section, we compare our results with those of previous works in terms of bow shock searches, the geometrical characterization of the IR structures, and the estimated ISM densities. We found nine new stellar bow shock candidates, three new bubble candidates, and one new intermediate structure candidate within the catalogs of O and Be-runaway stars compiled in Carretero-Castrillo et al. (2023). Among the catalogs, there were also 17 known bow shock candidates (see Bow. Ref. column in Tables 1 and 2). The mean $V_{\text{PEC}}^{2\text{D}}$ are 34 km s⁻¹, and 42 km s⁻¹, for the runaways associated with the new and the known IR structures, respectively. The mean E values are 1.53 and 1.84 for the runaway stars associated with the new and the known IR structures, respectively. The runaways associated with the known bow shocks have higher velocities and E values because they were the most obvious and easiest to find. In Carretero-Castrillo et al. (2023), we discover new runaways that are closer to the runaway threshold (E ~ 1, not as evident as the already known ones), and that thus have slower velocities. More than half of our new discoveries in the IR are associated with the new runaway discoveries in Carretero-Castrillo et al. (2023), and they are closest to the runaway threshold. Therefore, more detailed searches for runaway stars, such as the one presented in Carretero-Castrillo et al. (2023), working with space velocities and the accurate Gaia DR3 data, open the door to identifying both new runaway stars and bow shock structures around them.

Among the 17 known bow shocks found, we discarded the one associated with GR20 (EB43 in Peri et al. 2015), and reclassified the bubble of GR51 as an intermediate structure. In this work, we thus confirm 16 known candidates as bow shocks. The known bow shocks were identified for the first time in the following works: Peri et al. (2012), Peri et al. (2015), Kobulnicky et al. (2016), Maíz Apellániz et al. (2018), and Moutzouri et al. (2022).

We note that among the 106 O-type runaway stars we have found 9 new bow shocks plus 16 known ones, representing 23.6%. This is very similar to the percentage one can obtain from Peri et al. (2015) considering the sample of O-type runaways of Tetzlaff et al. (2011): 11 bow shocks among 54 runaways representing 20.4%. In contrast, for Be-type runaway stars we obtain 2 bow shocks among the 69 runaway stars in our sample, representing ~3%. By comparison, among the 12 Be runaway stars from Tetzlaff et al. (2011), there were no bow shocks identified by Peri et al. (2015).

From the references listed in the previous paragraph, only Peri et al. (2015) and Moutzouri et al. (2022) provide a geometrical characterization of the bow shocks in W4-band WISE images. In this work, we also include the geometrical measurements for both the bow shock and the bubble candidates. We also characterize those known bow shocks of our runaway catalogs that are not geometrically characterized in Kobulnicky et al. (2016), Maíz Apellániz et al. (2018), Bodensteiner et al. (2018), and Jayasinghe et al. (2019). The order of magnitude for R and w is less than 1 pc, while for l it is a few parsecs. We conclude that both new and known bow shocks have similar and compatible R, l, and w to those presented in Peri et al. However, in this work we have improved the procedure of measurement, using the SAOIm- ageDS9 software to measure l through ellipse shapes. We have also provided the eccentricity of the ellipses used to measure the bow shocks. We highlight the importance of the geometrical characterization of these objects, as it can be used to put constraints on the underlying radio emission mechanism of these sources (van den Eijnden et al. 2022b).

The surrounding ISM densities, $n_{\text{ISM}}^{2\text{D}}\text{ and }{n}_{\text{ISM}}^{3\text{D}}$, were estimated for 2D and 3D peculiar velocities, respectively, using Eq. (1), and are presented in Tables 4 and 5. The medians obtained for $n_{\text{ISM}}^{2\text{D}}\text{ and }{n}_{\text{ISM}}^{3\text{D}}$ are 6.4 and 3.8 cm⁻³, respectively. We note that these values were computed with our reduced sample of 12 bow shocks, for which we could estimate ISM densities. In Peri et al. (2012) and Peri et al. (2015), they computed n_ISM in 3D for 28 and 31 sources, respectively. The median n_ISM for all their 59 sources is ~2 cm⁻³. Despite our reduced sample, our 3D median and theirs are quite similar. Their sample is much larger, so their median should be more reliable than ours. However, we remark that our geometrical measures are more accurate than theirs (see Sect. 3). The typical densities of the cold neutral ISM have a wide dispersion, in a range from 1–100 cm⁻³ (see Kulkarni & Heiles 1988 for a review and the introduction of the recent paper Patra & Roy 2024). Our values fit well within this range, with two particular exceptions with densities of ≃200 cm⁻³, BS-GR72 and BS-GR83. We note that small R measures result in large n_ISM estimates. Therefore, projection effects may overestimate ISM densities. This is probably the case for BS-GR83. However, for BS-GR72 it seems to be a combination between its small R, high mass-loss rate, and wind terminal velocity. However, if BS-GR72 were an in situ bow shock, the n_ISM estimate would not be so relevant in this particular case. We emphasize again that the individual results in ISM density estimates should be used with caution because they are derived from a combination of Ṁ, $v_{\infty },V_{\text{PEC}}^{2\text{D}}\text{ or }{V}_{\text{PEC}}^{3\text{D}}$, and R values, for which: i) R is used instead of R₀ in Eq. (1), though in particular cases they might differ significantly for very small values of R, and ii) it is difficult to obtain accurate estimates of Ṁ and υ_∞ values.

6 Summary and future work

In this work, we have searched for stellar bow shocks around the new sample of O and Be runaway stars provided in Carretero- Castrillo et al. (2023), using WISE IR images, Gaia DR3 astrometric data, and radio data. The main conclusions of these work are summarized below.

While searching for stellar bow shocks around runaway stars, we found nine new stellar bow shock candidates, three new bubble candidates, and one new intermediate structure. Of these, ten are associated with O-type stars, and only three with Be-type stars. We also visually found a bow shock candidate in the BS- GR35 field, with no associated runaway star, which we called BS J075339–2616.5. We investigated and quantified the possible misalignment between the runaway star corrected proper motions and bow shock directions, and we found BS-GR72 as an in situ bow shock candidate related the H II region S206. We identified 17 known bow shocks, all of which are associated with O-type stars, but discarded the one associated with GR20 (EB43 in Peri et al. 2015). We identified 62 miscellaneous IR structures, 40 and 22 of which are associated with O- and Be-type stars, respectively. Among them, we introduced a new designation for bubble-like sources with around twice the W4 PSF: mini-bubble sources. We found two, BR10 and BR15, which are Be-type stars. We include under this new designation sources such as V452 Sct, of Maíz Apellániz et al. (2018). More detailed studies of these sources could reveal the nature of the IR emission around these runaways.

On the other hand, we conducted a geometrical characterization of the sources through IR measures in W4. For bow shocks, we measured the projected standoff distance, R, length, l, and width, w, while for bubbles we only provided the width. The values of R and w are always similar, and their order of magnitude is less than one parsec. The values of l are typically a few parsecs. The geometrical IR measures, together with mass-loss rates and stellar wind velocities of the runaway stars, allowed us to estimate the ISM density of the medium in a 2D- and 3Dpeculiar velocity case. We found $n_{\text{ISM}}^{2\text{D}}$ in the range from 0.8 to ~220 cm⁻³, with a median of 6.4 cm⁻³, and $n_{\text{ISM}}^{3\text{D}}$ in the range from 0.3 to ~200 cm⁻³, with a median of 3.8 cm⁻³. Particularly high values could be overestimated due to projection effects in the R measurement. Both the geometrical IR measures and the n_ISM estimates obtained are in agreement with previous works (Peri et al. 2012, 2015; Patra & Roy 2024 and references therein). We note that the obtained values of n_ISM, even considering the limitations outlined in Sect. 5, are estimated at the position of our sources. Therefore, they could be useful for other studies at those positions, such as stellar formation studies or radiation processes in these regions.

The extensive bow shock and bubble searches performed in works such as Peri et al. (2012), Peri et al. (2015), Kobulnicky et al. (2016), Bodensteiner et al. (2018), and Jayasinghe et al. (2019) leave few options for finding new ones from such sources in the sky. However, deeper searches for runaway stars using accurate Gaia DR3 data (Carretero-Castrillo et al. 2023) allow new runaway stars to be found closer to runaway thresholds, and then a still-hidden population of bow shocks. Most of the new bow shock candidates identified correspond to new runaway discoveries in Carretero-Castrillo et al. (2023). They also have smaller $V_{\text{PEC}}^{2\text{D}}$ velocities than known bow shocks. The differences are, on average, about ∼7 km s⁻¹. As for the bow shock percentage in runaways by spectral type, we obtained 23.6% and ∼3% for O and Be runaway stars, respectively. On the other hand, the main difference with previous works of general bow shock searches is that here we provide the geometrical IR characterization of the sources identified. These measurements are very useful to assess the underlying radio emission mechanism in these sources (see van den Eijnden et al. 2022b).

In this context, we searched for radio emission from these bow shocks in archival radio data and we evaluated the possible thermal or nonthermal origin of the radio emission mechanisms for our sources. More sensitive observations are needed to detect bow shocks at radio wavelengths, uncover a previously undetected population of radio-emitting bow shocks, and shed light on their radio and HE emission connection.

Data availability

Full Tables 1, 2, 4, 5 are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/694/A250

Acknowledgements

We thank the anonymous referee for useful suggestions and comments that helped to improve the content of the manuscript. This work made used of the public JUPYTER notebook of van den Eijnden et al. (2022b) to reproduce Fig. C.4. We thank the authors of that work for the data availability. We thank T. Antoja for useful discussions. We acknowledge financial support from the State Agency for Research of the Spanish Ministry of Science and Innovation under grants PID2019-105510GB-C31/AEI/10.13039/501100011033, PID2019-104114RB-C33/AEI/10.13039/501100011033, PID2022-136828NB- C41/AEI/10.13039/501100011033/ERDF/EU, and PID2022-138172NB- C43/AEI/10.13039/501100011033/ERDF/EU, and through the Unit of Excellence María de Maeztu 2020–2023 award to the Institute of Cosmos Sciences (CEX2019-000918-M). We acknowledge financial support from Departament de Recerca i Universitats of Generalitat de Catalunya through grant 2021SGR00679. MC-C acknowledges the grant PRE2020-094140 funded by MCIN/AEI/10.13039/501100011033 and FSE/ESF funds.This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of NASA’s Astrophysics Data System. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.

Appendix A WISE RGB images for new bow shock and bubbles

Fig. A.1 WISE RGB images in W4+W3+W2 for the new bow shock and bubble candidates presented in Table 4 in equatorial coordinates with North up and East to the left. Dashed white arrows indicate the directions of the original proper motions from Gaia DR3, and solid cyan arrows indicate the directions of the proper motions corrected for the ISM motion caused by Galactic rotation. Each field has a different size, but the arrows are fixed to 2′ . Each panel contains the name of the star as provided in the GOSC or the BeSS catalogs, and the runaway star id. The RGB images of the first three sources from Table 4 are presented in Fig. 2.

Appendix B Bow shock direction differences

As can be seen in some of the bow shocks shown in Figs. 2, A.1, and 4, the corrected proper motion direction of the runaway star is not well aligned with the bow shock direction. We tabulate and comment on these differences here.

For each runaway with an associated bow shock, either new or known, we list in Table B.1 the position angle, PA, computed from north to east of the corrected proper motion of the runaway star, as well as its corresponding uncertainty ∆PA, as explained in Sect. 3. We also list the bow shock angle, BS Angle, that we estimate visually from the W4 images and indicates the direction of the bow shock in the plane of the sky, together with an assumed systematic uncertainty that we call ∆BS Angle of 5° at 1σ. We also list the AD, computed as the absolute value of the difference between PA and BS Angle for each bow shock, together with AD Significance, which is the significance in σ units of the AD considering the previously quoted uncertainties.

For comparison, we list these values both for the new and the known bow shocks identified in our samples. For AD above 3σ, there are 4 new and 5 known bow shocks, while above 5σ, there are 3 new and 3 known bow shocks. This indicates that some bow shock candidates identified in previous studies also exhibit misalignments between their proper motion and bow shock directions. For the new bow shocks that have AD above 5σ, we comment on their possible origin as in situ bow shocks in Sect 4.1.

Appendix C Search for radio emission

C.1 Radio data

Stellar bow shocks emit mainly at IR but there is some observational evidence of radio emission coming from these sources (Benaglia et al. 2010; van den Eijnden et al. 2022a,b; Moutzouri et al. 2022). For the IR sources identified, we searched for observations carried out with the Very Large Array⁷ (VLA), the Giant Metrewave Radio Telescope (GMRT) (Swarup et al. 1991), and the Australian Square Kilometre Array Pathfinder (ASKAP) (Schinckel et al. 2012). In particular, we searched for the images provided by the NRAO VLA Sky Survey (Condon et al. 1998, NVSS), the VLA Sky Survey (Lacy et al. 2020, VLASS), and the Rapid ASKAP Continuum Survey, in both RACS-low, and RACS-mid images (Hale et al. 2021; Duchesne et al. 2024, respectively). In addition, in the case of GMRT and VLA archival observations, we looked for projects that might have produced related papers and images. Table C.1 presents the main properties of the radio surveys used in this work. These include their wavelength, angular resolution and sky coverage. We note that not all our sources are visible for each of the mentioned radio interferometers due to declination limits.

C.2 Search for radio counterparts

We performed a search for radio observations for the 17 objects listed in Tables 4 and 5, together with the BS J075339-2616.5 candidate. In the various radio archival data presented in C.1, we looked for radio emission that could be superimposed in the plane of the sky, either on the runaway stars or on the bow shock and bubble candidates, and focusing in the direction of the movement of the runaway star. No discrete radio emission was detected at the specific positions of the runaway stars. We include a summary of the findings for each of these sources in Table. C.2. For a few of them, the images presented reduction artifacts, precluding any findings (ill-imaged). We note that the fact that some emission was found does not indicate that it is a radio counterpart of the IR source. Out of the 18 objects, we comment on 6 cases where we found radio emission partially overlapping the IR source at some frequency (even if not significant in some cases). For the most interesting cases, we include the corresponding RACS images since they have better angular resolution than those of NVSS, and VLASS images are optimized for discrete rather than extended emission. The sources are presented in the same order as in Table C.2.

C.2.1 BS-GR9 (HD 41997)

In NVSS and RACS, we did not find emission at the bow shock position, but ~ 4′ southward from it. We found that this emission is not associated with the bow shock but to the Lower’s Nebula Sh2-261, which is an H II region (e.g., Wang et al. 2023, and references therein).

We highlight here that we found that BS-GR9 (and GR9) is contained in the uncertainty ellipse of the new Fermi source 4FGL J0609.1+1544 of the latest 4FGL-DR4 release (Abdollahi et al. 2022; Ballet et al. 2023). This High Energy (HE) gamma-ray source is not associated with any known source at other energies. Although this Fermi ellipse has a semimajor axis of 25′ and a semiminor axis of 18′, an association of this new HE source with the bow shock would be very interesting, because it could provide additional evidence of nonthermal processes in these sources (Benaglia et al. 2010). In this context, we note the attempts by H.E.S.S. to search for Very High Energy (VHE) gamma-ray emission from bow shocks of runaway stars (H. E. S. S. Collaboration 2018). Deeper radio observations than the available ones from RACS would be interesting in the search for nonthermal emission from this bow shock. Moreover, runaway stars can be found in binary systems containing compact objects (see Sect. 1), which can present nonthermal point-like radio emission and eventually HE and/or VHE gamma-ray emission. GR9 is not detected in the VLASS image with an RMS of ~0.2 mJy beam⁻¹. Thus, deeper observations could unveil the presence of a nonthermal radio source that could be connected to the HE gamma-ray one.

C.2.2 BS-GR72 (BD +50 886)

For the source BS-GR72 there are continuum images at 235, 325 and 610 MHz from GMRT data (Omar et al. 2002) showing a rather strong and extended source closely overlapping the extent of the WISE bow shock candidate. This radio source has been identified as the evolved H II region S206 (or NGC 1491), at a distance of ${2.96}_{-0.15}^{+0.17}$ kpc (Méndez-Delgado et al. 2022). At max imum angular resolution, its size is ~ 4′ or 740 pc. Omar et al. (2002) derived the electron temperature and emission measure of the source, and showed a spectral index (their Fig. 2) consistent with optically thick emission of such an object at the lower frequency bands, and a turnover around the largest frequency they observed, 610 MHz. The radio emission is also visible in the NVSS image. The VLASS images showed some emission at ~10″ scales superposed to the position of this bow shock, but with severe reduction artifacts. The runaway star GR72 is located at 2.94 ± 0.3 kpc (Carretero-Castrillo et al. 2023) in the same direction of S206. Therefore, this bow shock could be indeed an example of an in situ bow shock formed by ISM outflows from the H II region (see Sect. 4.1).

C.2.3 BS-GR75 (HDE 338 916)

For BS-GR75, the VLASS, RACS, and NVSS images show emission close to the bow shock position. However, it seems to come from a strong and extended source to the southeast, in the opposite direction of motion of the runaway. We then discarded that the mentioned emission can be a counterpart for our bow shock. That emission, in the case of VLASS, produces severe artifacts in the field, and through the bow shock position itself. The RACS fields are also quite noisy, so we did not include these radio images here. The RMS are 0.5 mJy beam⁻¹ and 0.2 mJy beam⁻¹ for RACS-low and -mid images, respectively. The mentioned strong radio emission hampers the possible detection of a weaker radio emission coming from our bow shock. Therefore, deeper radio observations, preferentially including also shorter baselines, could unveil a possible radio counterpart for BS-GR75 (see Sect. C.3).

C.2.4 BS J075339−2616.5 (−)

As explained in Sect. 4.1, this bow shock candidate was found in a different way than those presented in this work. It appeared in the BS-GR35 WISE field, at ~2′ south from the GR35 runaway. The RGB image in W4+W3+W2 centered and zoomed in BS J075339−2616.5 is presented in Fig. 3. The white cross in this figure denotes the center of the bow shock emission, whose ICRS coordinates, RA = 07h 53m 39.1s, DEC = −26° 16′ 31″, give it its name. No VLASS and only weak NVSS emission were detected at the position of this source. However, it has radio emission in both RACS−low and −mid images. Figure C.1 shows the RACS-low (top panel) and -mid (bottom panel) emission centered in BS J075339−2616.5. In RACS−low, the radio emission appears to be bright and extended, while in RACS-mid, the radio emission is rather compact and weak. The radio source in RACS-low is contained in the RACS Stokes I source catalog (Hale et al. 2021), which identifier is RACS_0741−25A_1914. The integrated flux is 6.8±2.1 mJy, and the root mean square (RMS) is ~0.3 mJy beam⁻¹. We computed the RACS-mid flux using a Gaussian function, and obtained 2.0 ± 0.5 mJy. We also estimated the RMS and found ~0.2 mJy beam⁻¹. Searching in SIMBAD at the position of BS J075339−2616.5 we only found the IRAS source IRAS 07515−2608. Therefore, there appear to be no known astrophysical sources in the region associated with this radio source. However, the peak position and the morphology of the radio source do not follow the ones of the IR emission, and we cannot confirm that it is associated with the BS J075339−2616.5 bow shock candidate.

C.2.5 BS/BU−GR51 (BD −08 4623)

Although there is no radio emission in VLASS, there is in NVSS and RACS images at the position of BS/BU-GR51. There are actually two sources, one at the position of the bow shock and another brighter one to the North; see the RACS-low and -mid images presented in top and bottom panels of Fig. C.2, respectively. The northern source appears in Hale et al. (2021) catalog, RACS_1824-06A_1653, and does not seem to be related to our bow shock. In RACS-low, there is radio emission close to the center of the bow shock. However, the morphology does not seem to follow that of the bow shock in W4. Therefore, we leave this source as a likely radio counterpart, but we cannot claim for a clear association. In any case, we measured the radio flux density and found 5.1 ± 0.1 mJy. In RACS-mid, there are just some peaks of radio emission at the position of the bow shock, so a possible radio counterpart to the bow shock is doubtful. The flux density measured in RACS-mid over the 3σ (RMS) contour is 0.9 ± 0.1 mJy (RMS = 0.3 mJy beam⁻¹). The fact that the emission decreases with increasing frequency could indicate a possible nonthermal origin of the radio emission.

Fig. C.1 WISE W4 image of the bow shock candidate BS J075339−2616.5. The color scale is WISE flux measured in profile-fit photometry, in units of Digital Numbers (DN; see https://wise2.ipac.caltech.edu/docs/release/allsky/expsup/sec1_4c.html#photcal). Top: RACS-low emission contours onto W4, levels at –0.7 (in white), 0.7, 2 and 3 mJy beam−1 (in black); synthesized beam of 25″ × 25″. Bottom: RACS-mid emission contours onto W4, levels at –0.5 (in white), 0.5 and 1 mJy beam−1 (in black); synthesized beam of 9.4″ × 8.3″.

C.2.6 BS-GR83 (HD 155775)

For BS-GR83, we did not find images from VLASS. In NVSS, there is no relevant emission superimposed to the bow shock. However, in RACS-low (see Fig. C.3 top panel), there is weak radio emission at the position of the bow shock. In RACS-mid (see Fig. C.3 bottom panel), there are only two very weak peaks. The measured RACS-low flux density is 2.5 ± 0.6 mJy; the RMS of RACS-mid image is 2 mJy beam⁻¹. We note that both radio maps are quite noisy. Therefore, we cannot claim a radio detection in this case, but there are some radio hints. As in BS-GR75, deeper radio observations could reveal radio emission, and, in this case, it may in fact be nonthermal (see Sect. C.3).

Fig. C.2 WISE W4 image of the bow shock candidate BS/BU-GR51. The color scale is WISE flux measured in profile-fit photometry, in units of Digital Numbers (DN; see https://wise2.ipac.caltech.edu/docs/release/allsky/expsup/sec1_4c.html#photcal). Top: RACS-low emission contours onto W4, with levels at –2.1 (in white), 2.1, 3, 3.5, 3.8 and 6 mJy beam−1 (in black); synthesized beam of 25″ × 25″. Bottom: RACS-mid emission contours onto W4, with levels at –0.5 (in white), 0.5 and 1 mJy beam−1 (in black); synthesized beam of 10.9″ × 8.2″. The cross indicates the position of the runaway star.

C.3 Testing radio emission mechanisms for bow shocks

Following van den Eijnden et al. (2022b), we can consider the nonthermal (synchrotron) and the thermal (free-free) scenarios to estimate the expected peak radio flux densities for the nondetected bow shocks. For this, we need the observables from radio images in conjunction with geometrical IR measurements (R, w), stellar parameters $\left(\dot{M},v_{\infty }\right)$, distances, and the estimated ISM density $n_{\text{ISM}}^{3\text{D}}$⁸. In our work, we already obtained all the necessary ingredients of the IR characterization for the bow shock candidates (see Tables 4 and 5).

To compute the radio observables, we used RACS-low DR1 images corresponding to a frequency of 887.5 MHz (34 cm). We accessed the images through the CSIRO Data Access Por- tal⁹. There are two types of RACS images: raw data and final data products. The raw data have reduced sensitivity towards the edge of the field, but provide higher resolution for any individual tile. The final data products are obtained by convolving to a common resolution of 25″. van den Eijnden et al. (2022b) worked with the raw data products because the final ones were not available at the time of writing. We worked with the two for comparison, but here we only present the raw data to be coherent with that work. The differences we found between the two data products are: elliptic and smaller beam sizes in raw images, and larger RMS values in final images. The RMS is always larger in final data products: depending on the field we found relative differences from 6 to 86%, with a mean value of 40%. In Table C.3, we present the radio observables needed from RACS-low (raw) images: the RMS measured in a circular region of 2′ centered at the bow shock position, the full image RMS, the beam major and beam minor axes, and the field name of the RACS image used.

Fig. C.3 WISE W4 image of the bow shock candidate BS-GR83. The color scale is WISE flux measured in profile-fit photometry, in units of Digital Numbers (DN; see https://wise2.ipac.caltech.edu/docs/release/allsky/expsup/sec1_4c.html#photcal). Top: RACS-low emission contours onto W4, with levels at –1.5 (in white), 1.5, 2.1 and 4 mJy beam−1 (in black); synthesized beam of 25″ × 25″. Bottom: RACS-mid emission contours onto W4, with levels at –0.5 (in white), 0.5 and 1 mJy beam−1 (in black); synthesized beam of 9.7″ × 8.5″.

In van den Eijnden et al. (2022b), a nonthermal scenario is set by assuming either a low magnetic field of B = 10 µG or the maximum magnetic field B_max so that the bow shock can form with a compressible stellar wind (see Eq. (2) in van den Eijnden et al. 2022b, and references therein). Other assumptions by these authors are: a high injection efficiency η_e of 10%, a maximum electron energy E_max = 10¹² eV, and a radio spectral index¹⁰ α = −0.5. van den Eijnden et al. (2022b) also considered a thermal scenario based on thin free-free emission. The electron number density n_e and temperature T of the electron population within a bow shock are the two primary physical properties that determine the thermal free-free emission from a shock. These parameters determine both the radio luminosity of the system and the surface brightness of H_α line emission. Using observational and literature constraints on these quantities, van den Eijnden et al. (2022a) inferred values on n_e and temperature T, which were compatible with the range of T = 6 × 10³ K or T = 1.4 × 10⁵ K found for a sample of H_α detected bow shocks in Brown & Bomans (2005). This was the temperature range employed by van den Eijnden et al. (2022b) and subsequently by us in this work. For n_e, it is expected a shock density enhancement in the bow shock. van den Eijnden et al. (2022b) assumed n_e = 4n_ISM, based on the 4 factor predicted by Rank- ine–Hugoniot equations (Landau & Lifshitz 1959) (for more details on the scenarios considered and/or the assumptions, see the original work by van den Eijnden et al. (2022b)).

Thanks to the data availability from van den Eijnden et al. (2022b), we can use their published code¹¹ to easily reproduce their results and integrate ours. To implement Fig. 10 of this work with our data, we only took the bow shock (not bubbles) candidates in Tables 4 and 5 that have estimated values of $n_{\text{ISM}}^{3\text{D}}$, RACS-low images available, and are non-detected bow shocks. They are BS-GR9, BS-GR35, BS-GR71, BS-GR75, BS-GR93; BS-GR40, BS-GR70, and BS-GR83. We excluded BS-GR51 since it has a clear radio detection, although it does not appear to be associated with the bow shock. We did not exclude BS-GR83 because its possible radio detection is not clear.

The predicted peak radio flux densities for our bow shocks in the nonthermal and thermal scenarios are presented in Fig. C.4. They are plotted as a function of three times the local RMS of each bow shock. Left panel corresponds to the nonthermal scenario, while right panel to the thermal one. The different markers show the different conditions in B or T of the scenarios, and also differentiate van den Eijnden et al. (2022b) sources of ours. There are two markers for each source: circles and squares representing the minimum and maximum considered values of B, respectively (left panel); and triangles and pentagons for the minimum and maximum considered values T, respectively (right panel). Therefore, for the nonthermal scenario, the minimum and maximum predicted peak radio flux density for each source is given by a circle and a square, respectively (except for 3 cases where van den Eijnden et al. 2022b found B_max < 10 µG). In the thermal scenario, the minimum and maximum predicted peak radio flux for each source is given by a pentagon and a triangle, respectively. In both panels the red line denotes where the two flux densities are equal, i.e., sources above that line are expected to be detected in RACS-low images. Thus, the figure gives an idea of how feasible it would be to detect a source in one scenario or the other. For the nonthermal scenario, 6 of our sources should be detectable for the maximum B (two yellow squares are overlapping), and two sources, BS-GR75 and BS- GR83, should be detectable for both magnetic fields. This is the main difference with respect to van den Eijnden et al. (2022b) sources, where they only found detectable sources for the maximum B. For the thermal scenario, we find the same two sources detectable for both electron temperatures, and BS-GR35 barely for the minimum one. The non-detection of BS-GR75 and BS- GR83 could be attributed to i) the significant fine tuning required for $\dot{M},v_{\infty }$, and n_ISM, to have a magnetic field close to B_max; ii) the lower limit of B = 10 µG is susceptible to small variations (for B = 8 µG the sources would not be detectable); iii) the possible overestimation of n_ISM that would leave the sources undetectable in the thermal scenario. Explanation iii) could play an important role when working with projected distances (as it is our case here), because of the inverse relation between R and n_ISM.

Fig. C.4 Predicted peak radio flux density as a function of three times the local RMS for the nonthermal (left) and thermal (right) scenarios. This plot contains the original data presented in van den Eijnden et al. (2022b), and our sources in yellow and magenta colors, with different symbols depending on the scenario conditions. The one-to-one line is shown by the red line; sources above it should, at minimum, have observable emission in a single beam.

The non-detection of BS-GR83 could be questionable, since it exhibits doubtful emission at the bow shock position (see Sect. C.2.6). A possible explanation for its unclear detection could be its noisy radio field (see Fig. C.3). In case the radio source is confirmed, it seems to fit better the nonthermal scenario, since the predicted flux density in the thermal scenario is strongly influenced by a very high $n_{\text{ISM}}^{3\text{D}}$ value. In the case of the non-detection of BS-GR75, the RACS-low image shows a strong radio source near its position. This bright source could hamper the detection of our source. Contrary to BS-GR83, in the case of BS-GR75 it is more difficult to choose one scenario or the other, since its $n_{\text{ISM}}^{3\text{D}}$ does not seem to be so clearly overestimated. More detailed radio observations on these two sources could probably lead to radio detections.

Lastly, we comment on BS-GR51. This source has a radio detection, and there are hints of nonthermal emission (see Sect. C.2.5). Moreover, its R is also small. The radio flux density obtained for this source is 5.1 mJy. Therefore, following the above work but in the inverse sense, and assuming a nonthermal scenario, we can estimate as 4.3 µG the magnetic field required to explain the BS-GR51 flux density. This is a reasonable magnetic field of the order of the typical value of 10 µG (Gabici et al. 2007; Martinez et al. 2023), and therefore the nonthermal scenario could explain the observed flux density. However, we note here again that it is not clear if this radio source is associated with the bow shock.

C.4 Summary of radio emission searches

We searched in different archives for observational evidence of radio emission coming from our sources. In particular, we used data from NVSS, VLASS and RACS radio surveys. For BS J075339−2616.5, and BS-GR51, we found radio emission coincident to the position of the bow shock. However, the morphology of the radio emission differs to some extent from that of the WISE bow shock. Therefore, we cannot claim that they have actual radio-emitting counterparts. For BS-GR75 and BS-GR83, there are some hints of radio emission, but their noisy fields make it difficult to assess whether they have real radio detections. BS-GR72 presented radio emission coincident to the bow shock position. However, it is associated with the H II region S206, and it is indeed thermal radio emission (Omar et al. 2002). Finally, for BS-GR9 we found radio emission but at ~ 4′ south from it. This radio emission is also associated with a H II region, in this case, the Lower’s Nebula Sh2-261 (Wang et al. 2023). We also investigated a thermal (free-free) and nonthermal (synchrotron) scenario for the origin of the radio emission. We computed a predicted radio flux density of our non-detected bow shocks. We found that two sources, BS-GR75 and BS-GR83, should be detectable in both scenarios and under all conditions considered. Their non-detection could be explained by the difficulty to reach B_max, small variations in the lower limit of B, and a possible overestimation of $n_{\text{ISM}}^{3\text{D}}$ . However, taking into account the geometrical characterization of these sources (small R), and the observational radio results, they could still be compatible with nonthermal radio-emitting sources.

All these leave us with a clear conclusion: more sensitive observations are needed to detect bow shocks at radio wavelengths, and clarify the physical mechanisms that could give rise to the radio emission. Deep radio observations with more sensitive instruments as VLA could uncover a new population of as yet non-detected radio bow shocks, and shed light on their radio and HE emission connection.

Appendix D Miscellaneous IR structures

The miscellaneous IR structures found in W4 around some runaway stars were introduced in Sect. 4.3. In this appendix we include the WISE RGB images in W4+W3+W2 of the sources presented in Table 6, except for the point-like sources. The images are presented in Fig. D.1. Some of the sources present complex IR fields (GR47, GR68), doubtful IR structures (GR38, GR101), some are mini-bubbles (BR10, BR15), and another presents a kind of bubble followed by a trail (GR81). Detailed multiwavelength studies on these sources could unveil the nature of the IR emission.

On the other hand, from the 62 miscellaneous sources we found, 55 presented nearly point-source IR emission. This emission could be associated with an IR excess from the circumstellar disk, especially in the Be-type stars (Rivinius et al. 2013).

Fig. D.1 WISE RGB images in W4+W3+W2 for the miscellaneous IR structures presented in Table 6 (except for the point-like sources) in equatorial coordinates with North up and East to the left. Dashed white arrows indicate the directions of the original proper motions from Gaia DR3, and solid cyan arrows indicate the directions of the proper motions corrected for the ISM motion caused by Galactic rotation. Each field has a different size, but the arrows are fixed to 2′ . Each panel contains the name of the star as provided in the GOSC or the BeSS catalogs, and the runaway star id.

References

All Tables

All Figures

	Fig. 1 Image of a known bow shock illustrating the geometric characterization used in this work. R is the projected standoff distance, l is the projected length of the bow shock, and w is its projected width. The arc of the ellipse shown with a continuous yellow line was used to derive l (see text for details).
In the text

Fig. 2 WISE RGB images in W4+W3+W2 for three examples of the new bow shock and bubble candidates presented in Table 4 in equatorial coordinates, with north up and east to the left. Dashed white arrows indicate the directions of the original proper motions from Gaia DR3, and solid cyan arrows indicate the directions of the proper motions corrected for the ISM motion caused by Galactic rotation. Each field has a different size, but the arrows are fixed to 2′. Each panel contains the name of the star provided in the GOSC or the BeSS catalogs, and the runaway star id. The remaining RGB images of the sources from Table 4 are presented in Appendix A.

In the text

	Fig. 3 WISE RGB image in W4+W3+W2 for the bow shock candidate BS J075339−2616.5 in equatorial coordinates, with north up and east to the left. The white cross indicates the center of the bow shock emission, at coordinates RA = 07h 53m 39.1s and Dec = −26°16′31″. This bow shock appeared in the WISE field of BS-GR35 (see Fig. 2), and has no associated runaway star.
In the text

Fig. 4 WISE RGB images in W4+W3+W2 for the known bow shock and bubble candidates presented in Table 5 in equatorial coordinates, with north up and east to the left. Dashed white arrows indicate the directions of the original proper motions from Gaia DR3, and solid cyan arrows indicate the directions of the proper motions corrected for the ISM motion caused by Galactic rotation. Each field has a different size, but the arrows are fixed to 2′. Each panel contains the name of the star as provided in the GOSC or the BeSS catalogs, and the runaway star id.

In the text

Fig. A.1 WISE RGB images in W4+W3+W2 for the new bow shock and bubble candidates presented in Table 4 in equatorial coordinates with North up and East to the left. Dashed white arrows indicate the directions of the original proper motions from Gaia DR3, and solid cyan arrows indicate the directions of the proper motions corrected for the ISM motion caused by Galactic rotation. Each field has a different size, but the arrows are fixed to 2′ . Each panel contains the name of the star as provided in the GOSC or the BeSS catalogs, and the runaway star id. The RGB images of the first three sources from Table 4 are presented in Fig. 2.

In the text

Fig. C.1 WISE W4 image of the bow shock candidate BS J075339−2616.5. The color scale is WISE flux measured in profile-fit photometry, in units of Digital Numbers (DN; see https://wise2.ipac.caltech.edu/docs/release/allsky/expsup/sec1_4c.html#photcal). Top: RACS-low emission contours onto W4, levels at –0.7 (in white), 0.7, 2 and 3 mJy beam−1 (in black); synthesized beam of 25″ × 25″. Bottom: RACS-mid emission contours onto W4, levels at –0.5 (in white), 0.5 and 1 mJy beam−1 (in black); synthesized beam of 9.4″ × 8.3″.

In the text

Fig. C.2 WISE W4 image of the bow shock candidate BS/BU-GR51. The color scale is WISE flux measured in profile-fit photometry, in units of Digital Numbers (DN; see https://wise2.ipac.caltech.edu/docs/release/allsky/expsup/sec1_4c.html#photcal). Top: RACS-low emission contours onto W4, with levels at –2.1 (in white), 2.1, 3, 3.5, 3.8 and 6 mJy beam−1 (in black); synthesized beam of 25″ × 25″. Bottom: RACS-mid emission contours onto W4, with levels at –0.5 (in white), 0.5 and 1 mJy beam−1 (in black); synthesized beam of 10.9″ × 8.2″. The cross indicates the position of the runaway star.

In the text

Fig. C.3 WISE W4 image of the bow shock candidate BS-GR83. The color scale is WISE flux measured in profile-fit photometry, in units of Digital Numbers (DN; see https://wise2.ipac.caltech.edu/docs/release/allsky/expsup/sec1_4c.html#photcal). Top: RACS-low emission contours onto W4, with levels at –1.5 (in white), 1.5, 2.1 and 4 mJy beam−1 (in black); synthesized beam of 25″ × 25″. Bottom: RACS-mid emission contours onto W4, with levels at –0.5 (in white), 0.5 and 1 mJy beam−1 (in black); synthesized beam of 9.7″ × 8.5″.

In the text

Fig. C.4 Predicted peak radio flux density as a function of three times the local RMS for the nonthermal (left) and thermal (right) scenarios. This plot contains the original data presented in van den Eijnden et al. (2022b), and our sources in yellow and magenta colors, with different symbols depending on the scenario conditions. The one-to-one line is shown by the red line; sources above it should, at minimum, have observable emission in a single beam.

In the text

Fig. D.1 WISE RGB images in W4+W3+W2 for the miscellaneous IR structures presented in Table 6 (except for the point-like sources) in equatorial coordinates with North up and East to the left. Dashed white arrows indicate the directions of the original proper motions from Gaia DR3, and solid cyan arrows indicate the directions of the proper motions corrected for the ISM motion caused by Galactic rotation. Each field has a different size, but the arrows are fixed to 2′ . Each panel contains the name of the star as provided in the GOSC or the BeSS catalogs, and the runaway star id.

In the text