© The Authors 2024

1 Introduction

Massive stars, including O-type objects and their evolved Wolf-Rayet (WR) counterparts, are known to produce strong stellar winds through the line-driving mechanism (Castor et al. 1975; Puls et al. 2008). These winds are made of a plasma at a temperature of a few 10⁴ K, actively producing thermal Bremsstrahlung in the radio domain. The thermal radio spectrum of massive star winds is known to be optically thick, with a dependence of the flux density (S_ν) as a function of frequency (ν) such as S_ν ∝ ν^α, and a spectral index α close to 0.6 (Wright & Barlow 1975; Panagia & Felli 1975). However, some deviations with respect to this behaviour have pointed to the existence of a complementary emission mechanism for a significant number of massive stars. This additional process is synchrotron radiation and it is characterised by a spectral index, α, that is significantly lower than the thermal value and which may even be definitely negative. Synchrotron radiation requires the existence of a magnetic field (of stellar origin) and a population of relativistic electrons. The latter is most likely produced by diffusive shock acceleration (DSA, Drury 1983). The occurrence of this process requires the presence of strong magneto-hydrodynamic shocks, and these are provided by the strong wind collisions that occur when two (or more) massive stars are part of a multiple system. Provided the stellar separation is high enough, the stellar winds collide at their terminal speeds, leading to velocity jumps of the order 2000–3000 km s⁻¹. The subset of colliding-wind binaries able to accelerate particles up to relativistic velocities is referred to as particle-accelerating colliding-wind binaries (PACWBs). To date, about 50 PACWBs have been identified, mainly thanks to the signature of synchrotron radio emission (De Becker & Raucq 2013)¹.

The interest in studying PACWBs is mainly twofold. First, it allows us to investigate shocks physics, including particle acceleration and non-thermal emission processes, in astrophysical laboratories different from supernova remnants (SNRs). The shock speeds, local magnetic field, and geometry are different, but the physical processes are the same. Second, particle acceleration in the pre-supernova massive star environment opens up the possibility of a moderate contribution to the production of cosmic rays. Even though SNRs are very likely to be the major contributors to galactic cosmic rays, we cannot reject the idea that a small contribution may come from alternative sources. The latter point calls upon an investigation of the population of particle accelerators, aimed at exploring the question of the fraction of PACWBs among massive binaries (De Becker et al. 2017, 2024).

As PACWBs are composite emitters (with thermal emission coming from individual winds and non-thermal emission from the colliding-wind region), it is not straightforward to disentangle their emission components. Most radio interferometers operating at centimetric wavelengths can reach angular resolutions of a few arcseconds, but that is not enough to resolve the wind collision from the inner parts of the winds. The only way to spatially resolve the synchrotron emission from the thermal emission is the use of very long baseline interferometry (VLBI), with baselines of thousands of kilometers. The angular resolution at centimetric wavelengths can be as good as a few milli-arcseconds (mas), thereby allowing for binary systems with orbital period of at least a few years at distances of a few kpc to be spatially resolved.

To date, only a handful of PACWBs have benefitted from such high angular-resolution measurements. In particular, we note WR 147 (Williams et al. 1997), WR 140 (Dougherty et al. 2005), WR 146 (O'Connor et al. 2005), Cyg OB2 #9 (Dougherty & Pittard 2006), Cyg OB2 #5 (Ortiz-Leόn et al. 2011), HD 93129A (Benaglia et al. 2015), and, more recently, Apep, the most powerful synchrotron emitter in this category (Marcote et al. 2021).

This study is aimed at investigating the synchrotron emission region from two PACWBs: HD 167971, and HD 168112. The goal is to make use of the European VLBI Network (EVN) to obtain valuable radio images of both systems, with the aim to characterise their synchrotron emission at the mas scale. Even though VLBI techniques have already been applied to HD 167971 (Sanchez-Bermudez et al. 2019), this paper reports on the very first high-resolution radio imaging of HD 168112. The paper is organised as follows. The two targets are introduced in Sect. 2. The data acquisition and processing are described in Sect. 3. Our results are presented in Sect. 4 and discussed in Sect. 5. Finally, we present our summary and conclusions in Sect. 6.

2 Targets: Two well-established PACWBs

The two systems investigated in this paper are located in the NGC 6604 open cluster (Reipurth 2008), at a distance of about 1.7 kpc. Both are part of the catalogue of PACWBs (De Becker & Raucq 2013) and have already been targets of several observation campaigns across the electromagnetic spectrum.

HD 167971 (BD-12 4980) is a hierarchical triple system made of a 3.32 d period system (O6-7V + O6-7V) with a third object (O8I) on a wider (close to 20-yr) orbit (Ibanoglu et al. 2013). Adopting the naming convention proposed by De Becker (2015), the short period system is made up of components Aa and Ab, while the third star is component B. After a first confirmation of an orbital motion between components A and B using infrared interferometry (De Becker et al. 2012), the very first three-dimensional (3D) characterisation of the wide orbit has been achieved thanks to a combination of interferometry and spectroscopy by Le Bouquin et al. (2017). The period of the wide orbit is 7803 ± 540 d, with a significant eccentricity (e = 0.44 ± 0.02).

The radio emission at centimetre wavelengths was investigated in detail by Blomme et al. (2007), confirming the synchrotron emitter status through all main usual criteria (significant variability, non-thermal spectral index, and high brightness temperature). HD 167971 is the brightest synchrotron emitter among O-type PACWBs. The light curves present clear evidence for a strong variability of the radio emission compatible with the long period in the system. The maximum flux density was measured in 1988, which is coincident with the expected periastron passage according to the ephemeris published by Le Bouquin et al. (2017). This is in agreement with the expectation that the maximum of the synchrotron emission should occur close to periastron passage (see Sect. 5.1.1). More recently, Sanchez-Bermudez et al. (2019) published the results of the first VLBI observations of HD 167971 at two epochs: 2006 and 2016.

HD 168112 (BD-12 4988) displays strong hints for binarity, as revealed by radio and X-ray measurements (De Becker et al. 2004), even though no orbital solution exists at the time of writing. Even the period of the system is still completely unknown; however, high angular-resolution imaging in the visible could identify the presence of a visual companion at an angular distance of about 3.3 mas (Sana et al. 2014). The flux ratio between the primary (O5.5III, component A) and the secondary allowed De Becker (2015) to propose an O5.5–7.5III spectral type for the secondary (component B). On the basis of the variation of the X-ray emission, produced in part by the wind-wind interaction region, a conservative lower limit on the eccentricity of 0.38 could also be determined (De Becker 2015).

The radio emission from HD 168112 shows clear evidence for synchrotron radiation, namely, non-thermal spectral index and variability (De Becker et al. 2004; Blomme et al. 2005). However, this system has not been the target of high angular resolution imaging measurements prior to the present study.

3 Observations and data reduction

Our observations with the EVN were performed on 5 November 2019 10:00–17:30 UTC at 18 cm (1.67 GHz) with project code EM131F (PI: B. Marcote), with the following 15 participating stations: Lovell Telescope, Effelsberg, Westerbork single dish, Medicina, Noto, Onsala 32-m, Tianma, Urumqi, Toruń, Hartebeesthoek, Svetloe, Zelenchuskaya, Badary, Irbene, and Sardinia. The data were correlated with the SFXC software correlator (Keimpema et al. 2015) at the Joint Institute for VLBI ERIC (JIVE; The Netherlands) with the total bandwidth of 256 MHz divided in eight 32-MHz subbands, with 64 spectral channels each, full circular polarisation, and 2-s time integration.

The source J1743–0350 was used as fringe finder and bandpass calibrator. J1825–0737 was used as phase calibrator for both target sources, HD 167971 and HD 168112, in a phase-referencing cycle of 3.5 min on target and 1.5 min on the calibrator. The targets were alternated on each phase-referencing cycle. As a result, HD 167971 and HD 168112 were observed for a total of ≈2 h each. We note that the angular separation between the phase calibrator and the targets is 5.0° and 4.8°, respectively.

The EVN data have been reduced in AIPS² (Greisen 2003) and DIFMAP (Shepherd et al. 1994) following standard procedures. An a priori amplitude calibration was performed using the known gain curves and system temperature measurements recorded on each station during the observation as provided by the EVN Pipeline. We corrected dispersive delays for the ionosphere using the maps of total electron content provided by the Global Positioning System satellites via the TECOR task in AIPS. We manually removed bad data (mainly frequencies and times affected by radio frequency interference). We first corrected for the instrumental delays and bandpass calibration using the fringe-finder source and thereafter conducted the global fringe-fit of the data using all calibrator sources. The phase calibrator was then imaged and self-calibrated to improve the final calibration of the data in DIFMAP. The obtained solutions were transferred to the target sources, which were then imaged.

The cleaned images were obtained by model fitting a 2D Gaussian to the (u, v) data. Due to the limited signal-to-noise ratio (S/N) in the images and sparse (u, v) coverage, this approach produced the optimal characterisation of the sources while retaining a minimum number of free parameters. The resulting images then allowed us to characterise the geometry of the wind-collision region for each target.

Fig. 1 Field of HD 167971 (left) and HD 168112 (right) as seen by the EVN on 5 November 2019 at 18 cm. Both sources are detected in our data with a significance of 7σ and 9σ above the rms noise level of 140 and 70 μJy beam−1, respectively. Contours start at 3σ the rms noise level and increase by factors of $\sqrt{2}$. The synthesised beam size of 12 × 22 mas2, PA = 3°, is represented at the bottom left corner of each image as a gray ellipse.

4 Results

Both sources, HD 167971 and HD 168112, were detected as compact sources at the mas scale in our observation on 5 November 2019. In the following, we present the results obtained from the aforementioned EVN data, independently for each target source. These results are also summarised in Table 1 and the resulting images of both targets are shown in Fig. 1.

HD 167971 is detected as a compact radio source at the position (J2000) of α = 18^h18^m5.89300^s ± 1.1 mas, δ = −12°14′33.3672″ ± 1.1 mas. We note that the quoted uncertainties represent the 1σ confidence interval and take into account the statistical uncertainties of the position measured in our data, the uncertainty in the position of the phase calibrator source, J1825–0737, within the international celestial reference frame (ICRF)³, and the estimated uncertainties associated with the phase-referencing technique (Kirsten et al. 2015).

The radio source was detected with a peak brightness of 0.97 ± 0.14 mJy beam⁻¹. A Gaussian fit to the emission reports a total flux density of 3.42 ± 0.14 mJy extended over a linear extent of 34.8 ± 0.3 mas with a position angle of 78 ± 7°. The 2D Gaussian fit converted to a 1D Gaussian; namely, the source size in one of the dimensions is negligible when compared to the resolution of the images (12 × 22 mas², PA = −2.7°). We note that the absolute amplitude calibration of the EVN data can exhibit uncertainties of ~15% due to the intrinsic calibration procedures related to VLBI arrays. These uncertainties were added in quadrature to the aforementioned value and the final flux density measurement is quoted in Table 1.

HD 168112 is detected as a radio source with a flux density of 1, 65 ± 0.12 mJy and peak brightness of 0.65 ± 0.07 mJy beam⁻¹ (prior absolute uncertainties in the absolute flux density scale) at the position (J2000) of α = 18^h18^m40.86587^s ± 1.0 mas, δ = −12°6′23.4233″ ± 1.0 mas. A Gaussian fit to this source provided a measurement of the source size of 28.3 ± 0.6 mas, elongated towards the direction 61.7 ± 1.5° (measured from the northern direction to the east). The source size in the orthogonal direction was negligible when compared to the synthesised beam size.

Figure 1 shows the obtained images for both HD 167971 and HD 168112. In the case of HD 167971, we detected a compact radio source at the center of the map, while the elongated structure observed at its north-eastern part is ambiguously connected to the source. Indeed, the 3–4σ blobs detected in the South left and right corners of the images (similar in position for both images) are indicative of some phase losses and possible flux from the target that got smeared into the map. These structures actually disappear when performing a self-calibration of the data based on the fitted Gaussian model, indicative of the fact that they are not astrophysical signals. This is not surprising in the context of these data due to the following factors: the targets exhibit a low declination that made them visible only at low elevations for most of the EVN array. During the observation, both targets were observed at an average elevation of ~20°. Furthermore, the phase calibrator is located at a significant separation with respect to the targets (~5°), making the phase-referencing technique less accurate when transferring the calibration from the phase calibrator to the targets (i.e. the assumption that all sources require the same corrections because the angular separation between them is negligible is slightly deviated). Finally, this observation was performed at a data rate of 2 Gbps (total bandwidth of 256 MHz) after the increase of bandwidth in the EVN. While the larger bandwidth is indeed beneficial to improve sensitivity, at the observed frequency (1.67 GHz), the presence of radio frequency interference (RFI) has significantly increased in the last years. At this band, several of the subbands were strongly affected by RFI, producing major changes in the gains of the telescopes⁴. This resulted into a larger scatter in amplitudes than usual, producing the aforementioned artifacts in the images and a potential loss of part of the source fluxes. The final flux density values quoted in Table 1 take this into account. We clarify that these unexpected features, namely, the north-easterly elongated structure, along with the southern ones, do not appear at all in the radio images of HD 167971 published by Sanchez-Bermudez et al. (2019) on the basis of their VLBA data. Furthermore, these structures completely disappear when we perform a self-calibration on the data, confirming that they are originated by phase losses during the calibration process. As a result of these artifacts in the data, the geometry in the source structures at the 3–4σ levels are not fully reliable. This is the reason why instead of conducting a fit assuming a bow-shaped structure as expected in these sources (see e.g. Benaglia et al. 2015; Marcote et al. 2021), we limited the analysis to a 2D Gaussian. This approach allows us to constrain the position angle of the emission, without being biased to possible curvature in the source structure that is not real.

5 Discussion

5.1 Synchrotron emission

5.1.1 HD 167971

In this triple system, two colliding-wind regions are active: between components Aa and Ab, and between components A and B, respectively. Given the size of the orbit of the AaAb binary, it is very likely that any putative synchrotron emission from there would be completely suppressed by free-free absorption (FFA). However, the synchrotron emission that is measured certainly comes from the wide orbit, in agreement with the light curve obtained by Blomme et al. (2007). On the basis of the ephemeris published by Le Bouquin et al. (2017) for the wide orbit, the orbital phase of our EVE measurement is ϕ = 0.52. Our observation thus occurred very close to apastron. Not counting the effect of FFA, we would expect the synchrotron emission to be at minimum at apastron. The reason is twofold. First, the local magnetic field in the emission region (of stellar origin) decreases as the distance between the stars increases (Usov & Melrose 1992). Second, when the separation is longer the density of the colliding flows is lower, very likely leading to a drop in the injection rate of particles into the DSA process. We can thus expect our measurement to be at the minimum of the light curve. This is indeed confirmed when checking the light curve published by Blomme et al. (2007), which presents a minimum compatible with an expected time of apastron passage. As a result, the flux density from the synchrotron emission region at any other orbital phase should be greater than our 2019 measurement.

Our radio data alone are not enough to determine the position of the stars relative to that of the detected source. However, we can discuss some expectations based on the wind parameters of the two components. According to Le Bouquin et al. (2017), the projected semi-major axis of the 21-yr orbit is about 18.15 mas. With an eccentricity of 0.443, this leads to a stellar separation at apastron of about 26 mas, with a line of centers very close to the east-westerly direction. Based on the wind parameters quoted in Table 2, the wind momentum rate ratio $\left(\eta ={\dot{M}}_{\text{B}}{v}_{\infty ,B}{\dot{M}}_{\text{A}}^{-1}{v}_{\infty ,\text{A}}^{-1}\right)$ is about 0.38. This leads to an angular separation between the wind collision and component A (resp. B) of about 10 mas (resp. 16 mas). These numbers are a bit below and above the east-west width of the synthesised beam (12 mas), respectively. In Fig. 1 (left), this translates to a potential position component A a bit beyond the left (eastern) limit of the second contour of the radio source, and component B should be a bit further on the right (western) of the source. However, we recommend some caution as these considerations are based on assumed wind parameters and not on a direct astrometric determination of the position of the stars.

A comparison with the measurements published by Sanchez-Bermudez et al. (2019) is also relevant. We have to extrapolate flux densities at the same frequencies to proceed with a valid comparison. The 2006 and 2016 measurements are summarised in Table 3. However, it is important to note that the spectral indices at both epochs point to a surprisingly steep spectrum, much steeper than the expectation for a standard optically thin synchrotron spectrum from a population of relativistic electrons accelerated through DSA in high Mach number shocks, in the test-particle regime (α = −0.5). On the basis of the spectral index map, Sanchez-Bermudez et al. (2019) considered that an α value of about −1.1 is more typical of a significant part of the synchrotron emission region. They interpreted the steepness of the synchrotron spectrum as a likely consequence of efficient inverse Compton cooling, leading to a softening of the relativistic electron spectrum.

We also want to stress that such a steep spectral index could also be a signature of shock modification. Basically, the backstreaming of high-energy particles upstream provides a contribution to the compression, leading to a drop in upstream velocity right behind the shock front. This creates a shock precursor that constitutes a main component of the modified shock structure. When relativistic electrons are backscattered by scattering centers located in the precursor and not in the far-upstream region, they experience a lower velocity jump that translates into a lower compression ratio, as compared to high Mach number adiabatic shocks for monoatomic gas (χ < 4). The mean free path of relativistic particles diffusing in the magnetised plasma is proportional to the gyroradius, depending itself on particle energy. As a result, lower energy relativistic electrons diffuse across shorter distance upstream and do not go beyond the shock precursor. As the electron index (p) depends on the compression ratio, $$p=\frac{\chi +2}{\chi -1},$$(1)

lower energy electrons are characterised by a greater p value, leading to a relativistic electron population that is steeper than expected from linear DSA by unmodified shocks. In a given magnetic field, the typical synchrotron photon frequency is directly proportional to the relativistic electron energy. This back reaction of relativistic particles on the shock structure is thus expected to lead to a steeper synchrotron spectrum at low frequencies. This is typically what is observed in the case of young supernova remnants where particle acceleration efficiency is high enough to push DSA into the non-linear regime, leading to steep synchrotron radio spectra (Dubner & Giacani 2015).

Assuming a spectral index of −1.1, we can extrapolate the flux density measured with the EVN to the bands used by Sanchez-Bermudez et al. (2019). We obtain flux densities of ~1.0 and ~0.58 mJy at 5.0 and 8.4 GHz, respectively. These values are indeed significantly lower than measured in 2006 and 2016 epochs, in agreement with the notion that our measurement occurs close to apastron, when the synchrotron emission should be minimum. We may consider the fact that our measurement close to apastron provides us with the bottom flux value of the full radio light curve.

Our radio measurement ought to be compared to the expected thermal free-free emission from the winds of the system. We adopted the approach described, for instance, by Wright & Barlow (1975) to calculate reasonable values of the thermal flux density for the three stellar winds contributing to this triple system. In agreement with Leitherer et al. (1995), we assumed values of the molecular weight, the mean electron number, and the RMS ionic charge of 1.8, 1.0, and 1.0, respectively. Focusing on the outer parts of the wind where the measurable thermal radio emission is expected to come from, we assumed a clumping factor of the order of 4.0 (Runacres & Owocki 2002) and an electron temperature equal to 30% of the stellar effective temperature (Drew 1990). All useful parameters are summarised in Table 2. We thus predict thermal flux densities at the μJy level for all component at 18 cm. These values are about three orders of magnitude lower than our actual measurement. Let’s also mention that in principle some thermal emission is also expected to arise from the colliding-wind region (Pittard 2010). The shocked wind material is a hot plasma that is likely to lead to some free-free emission. A significant thermal contribution from this region would require sufficiently high densities (emission process proportional to the square of the density) that could appear in very short binaries. However, for long-period systems, the winds are colliding at a much longer distances, leading to a severe drop in the density that is scaling as 1/r², with r being the distance from the star to the wind collision. According to the simulations by Pittard (2010), we could expect a significant contribution from colliding wind regions of radiative shocks that occur in quite short period systems, especially at radio and sub-millimeter frequencies well above our measurement frequency. On top of that, we cannot call upon any free-free emission enhancement due to the clumpiness of the emitting material, as most likely density structures do not survive the shocks and dissipate in the post-shock region of adiabatic shocked winds (Pittard 2007). At 1.67 GHz, for a wide system producing adiabatic shocks such HD 167971, the thermal emission arising from the colliding-wind region is certainly even lower than that coming from the individual winds. The radio source we measured with the EVN is thus by no means compatible with thermal emission. This lends further support to the idea that our EVN source is attributable to the synchrotron emission region.

5.1.2 HD 168112

The lack of existing orbital solution for this system prevents us from interpreting our radio measurement as a function of the orbital phase. However, the striking and important result is the very first imaging of the synchrotron emission region. In the absence of any other VLBI measurement at other epochs and at other frequencies for this source, we must be careful with our interpretation with respect to the nature of the measured radio emission. In other words, we have to be sure that we are not dealing with thermal radio emission.

First of all, we predicted the expected thermal free-free emission flux density, adopting the same approach as for HD 167971 above. Our predictions lead to values of the order of ~10 μJy (see Table 2). These values are much lower (almost a factor 200) than our measurement of 1.7 ± 0.4 mJy. The predicted thermal emission is clearly below our detection threshold. As in the case of HD 167971, we are dealing with a rather long period system producing adiabiatic shocks, leading therefore to a negligible putative thermal emission contribution from the stellar winds, especially at a frequency as low as 1.67 GHz. Second, the morphology of the source is significantly elongated, with the elongation of the emission region not coincident with the major axis of the synthesised beam. This is in agreement with expectation that the synchrotron emission region should be somewhat extended, coincident with the colliding-wind region. Once again, this measured morphology is not at all compliant with the expectation of a purely unresolved point source for thermal stellar winds. These two facts can fully be interpreted in terms of a synchrotron emission region at the limit to be resolved. These features are not at all compatible with thermal emission from the winds.

The measured position angle for this elongation (61.7 ± 1.5°) constrains the expected position for the two stars of the system. We would thus expect these stars to be placed within a position angle (at the epoch of the observation, 5 November 2019) of ≈152°. Unfortunately, the angular separation between the stars cannot be constrained with these data.

5.2 Energy budget

Given the negligible contribution of thermal emission at our measurement frequency, the radio flux densities reported here can be seen as clean and direct measurement of the synchrotron emission from the two targets. The synchrotron emission discussed in Sect. 5.1 actually results from a chain of energy conversion processes, on top of which lies the kinetic power of the stellar winds involved in the shock physics responsible for particle acceleration. The kinetic power, expressing the rate of transfer of mechanical energy by an individual stellar wind, can be expressed as: $$P_{\text{kin }}=\frac{1}{2}\dot{M}{v}_{\infty }^2=3.155\times {10}^{35}{\dot{M}}_{-6}{v}_{\infty ,8}^2\left({\text{ergs}}^{-1}\right)\text{, }$$(2)

where ${\dot{M}}_{-6}$ is the mass loss rate in units of 10⁻⁶ M_⊙ yr⁻¹, and v_∞,8 is the terminal velocity in units of 10⁸ cm s⁻¹. Using the wind parameters quoted in Table 2, we determined the total kinetic power for both systems using Eq. (2), resulting from the sum of individual contributions from all components.

The synchrotron luminosity results from the integration of the flux density in a given spectral domain: $$\begin{array}{cc}{L}_{\text{synch}}& =4\pi d^2\underset{v_{\min }}{\overset{v_{man}}{{\displaystyle \int }}}{S}_v\text{d}v\\ & =4\pi d^2{S}_{v,0}\frac{v_0^{-\alpha }}{\alpha +1}\left[v_{\max }^{\alpha +1}-v_{\min }^{\alpha +1}\right]\left({\text{ergs}}^{-1}\right),\end{array}$$(3)

where the flux density dependence on the frequency has been defined as S_ν ∝ v^α, while d is the distance, and ν₀ is the frequency of our measurement. The lower integration boundary can be set to 0.1 GHz as it is very likely that the synchrotron emission will be self-absorbed below that value. For the upper boundary, a conservative value of 100 GHz was adopted, as above that value the synchrotron flux is not expected to contribute much, given the negative index of the spectrum.

Assuming α = −0.5 for a standard population of relativistic electrons accelerated by DSA in high Mach number shocks, we computed L_synch using Eq. (3) for both systems and obtained values of the order of a few times 10²⁹ erg s⁻¹ (see Table 4). We can then determine the radio synchrotron efficiency (RSE) defined as the fraction of the wind kinetic power converted into synchrotron radiation, namely, RSE = L_synch/P_kin. Our results lie between the boundaries presented by De Becker et al. (2017) and quoted as RSE_lower and RSE_upper in Table 4, in full agreement with expectations for the radio emission from PACWBs.

In the case of HD 167971, as our measurement occurred close to apastron, the L_synch and RSE values quoted here should be considered as minimum values. In addition, in that part of the orbit the separation is so large that it is unlikely that the current measurement is affected by significant FFA. As a result, our energy budget estimate should be considered as representative of the actual synchrotron production at that orbital phase. For HD 168112, the lack of any clear information on the orbit prevents us from interpreting further our results. Depending on the size and phase of the orbit (that is known to be quite eccentric, see Sect. 2), we may expect some free-free absorption to occur, especially close to periastron. Thus, we cannot claim that our estimate of L_synch (and RSE) is typical of the actual synchrotron radiation production. Our measured quantities should thus be considered as lower limits on the synchrotron radiation production rate at that specific epoch.

For any specific PACWB at any orbital phase, RSE depends on several energy conversion factors: (i) the fraction of the kinetic power actually invested in shock physics, (ii) the fraction of that energy injected into particle acceleration, (iii) the fraction of energy in non-thermal particles that is going to relativistic electrons, and, finally, (iv) the fraction of relativistic electron energy that is radiated as synchrotron radiation. Even though it sounds reasonable to consider that most of these energy conversion factors may be similar from one system to the other, they cover a parameter space that is large enough to allow RSE to be significantly different among the population of PACWBs. A deeper investigation of these energetic considerations is required for a larger sample of systems, at various orbital phases and on the basis of high quality radio measurements to better quantify their potential to convert wind kinetic energy into synchrotron radiation.

5.3 The synchrotron-binarity correlation

Although this study is dedicated to the study of synchrotron emission from massive stellar systems, it is important to place it in a broader context. It is clear from the catalogue of PACWBs (De Becker & Raucq 2013, including its updated on-line version), that there is a strong correlation between the PACWB status (mainly based on the identification of synchrotron emission) and the binary (or higher multiplicity) status. On the one hand, this constitutes an incentive to search for new PACWBs focusing on already identified binaries. On the other hand, however, this opens up the possibility to make use of tracers of synchrotron radiation to establish the binary status of massive stars that are still elusive from the perspective of their multiplicity study using usual techniques. In other words, given that the so-called synchrotron-binarity correlation is quite well-established, synchrotron emission should also be considered as a tracer of multiplicity among massive stars.

In this framework, HD 168112 is a quite instructive example. Despite some hints for binarity from its behaviour in X-rays (De Becker 2015) and the identification of an astrometric companion with no monitoring to trace a potential orbital motion to date (Sana et al. 2014), this object has never revealed any radial velocity variation suggestive of an orbital motion. Adding the radio variation to this picture, we can talk about converging hints pointing to a binary status. The high angular resolution imaging we propose here, as a unique snapshot measurement, is actually enough to ascertain the binary status of the object.

This opens up the possibility to use the method presented in this paper to identify binaries among massive stars. However, we must clarify that this approach is valid, provided the synchrotron emission is not overly free-free absorbed by the stellar winds. This puts some significant constraints on the binary parameter space likely to be investigated. Basically, only sufficiently long period systems are relevant, namely, with periods of at least several months, or even years in the case of thicker stellar winds of Wolf-Rayet stars responsible for FFA up to longer distances. This fact is also quite interesting, given that longest period systems are in general more difficult to identify (a requirement for a long monitoring is to identify variations, eccentric orbits showing significant radial velocity variations in a narrow range of orbital phases, with slower orbital motion reducing the amplitude of the radial velocity curve, etc.). Thus, provided that the high angular resolution radio measurement does not occur at a phase of strong FFA (i.e. periastron), it is likely to reveal features directly pointing to a binary status. Of course, a unique measurement is not enough to access a full description of the orbit, but it is worthwile to use this approach to first identify as-yet-unrevealed binaries, before organising more focused campaigns using spectroscopic or astrometric techniques with the aim of characterising the orbit.

6 Summary and conclusions

We reported on the results of VLBI measurements of two well-established PACWBs, members of the NGC 6604 open cluster: HD 167971 and HD 168112. We used the EVN at 18 cm to image the radio emission from the two systems at one epoch, 5 November 2019. In both cases, we confidently detected a radio source, which we attribute to synchrotron radiation from the colliding-wind region.

In the case of HD 167971, our measurement is coincident with apastron, where the synchrotron radiation is at its lowest level. The source is unresolved and too bright to be explained by thermal free-free emission from the stellar winds of the components of the system. For HD 168112, we obtained the very first image of the synchrotron emission region. The radio source appears a bit elongated with a position angle that is significantly different from that of the synthesised beam, lending support to a slightly resolved synchrotron source. This makes HD 168112 the eighth known PACWB with a synchrotron emission region imaged using high angular resolution radio measurements, joining the ranks of: WR147, WR140, WR146, Cyg OB2 #9, Cyg OB2 #5, HD 93129A, HD 167971, and Apep.

Based on energy budget considerations, we quantified the synchrotron emission in terms of the fraction of the wind kinetic power effectively converted into synchrotron radiation. In the case of HD 167971, the lack of significant FFA makes our estimate typical of the actual synchrotron radiation production rate. However, in the case of HD 168112, the lack of information on the size and phase of the orbit prevents us from being more specific in our interpretation. In both cases, the mechanical-to-synchrotron energy conversion is fully compliant with expectations for PACWBs.

Finally, we stress the relevance and importance of snapshot high angular resolution radio imaging in the identification of synchrotron radio emitters. This is particularly true of HD 168112, which displays signs of binarity, although it has not yet provided us with the ultimate proof. The clear detection of the radio source with a flux density much greater than expected from the thermal emission from stellar winds, along with a morphology deviating from that of an unresolved point source, is enough to confirm the synchrotron nature of the emission; that, in turn, is a clear sign a colliding-winds in a massive binary system. As a consequence, on top of being an efficient tool to identify synchrotron radio emission from massive stars, high angular resolution radio imaging constitutes a highly relevant method for unveiling binary systems that are still elusive from the perspective of usual multiplicity study techniques. The latter point is especially relevant for the identification of long-period binaries.

Acknowledgements

The authors would like to thank the referee for a fair report and constructive comments that helped to improve the paper. This research is part of the PANTERA-Stars collaboration, an initiative aimed at fostering research activities on the topic of particle acceleration associated to stellar sources (https://www.astro.uliege.be/~debecker/pantera). The European VLBI Network is a joint facility of independent European, African, Asian, and North American radio astronomy institutes. Scientific results from data presented in this publication are derived from the following EVN project code: EM131. The SIMBAD database was used for the bibliography. B.M. acknowledges financial support from the State Agency for Research of the Spanish Ministry of Science and Innovation under grant PID2019-105510GB-C31/AEI/10.13039/501100011033 and through the Unit of Excellence Maria de Maeztu 2020-2023 award to the Institute of Cosmos Sciences (CEX2019-000918-M).

References

All Tables

All Figures

Fig. 1 Field of HD 167971 (left) and HD 168112 (right) as seen by the EVN on 5 November 2019 at 18 cm. Both sources are detected in our data with a significance of 7σ and 9σ above the rms noise level of 140 and 70 μJy beam−1, respectively. Contours start at 3σ the rms noise level and increase by factors of $\sqrt{2}$. The synthesised beam size of 12 × 22 mas2, PA = 3°, is represented at the bottom left corner of each image as a gray ellipse.

In the text