A VLBI study of the wind-wind collision region in the massive multiple HD 167971

Context . Colliding winds in massive binaries are able to accelerate particles up to relativistic speeds as the result of the interaction between the winds of the di ﬀ erent stellar components. HD167971 exhibits this phenomenon which makes it a strong radio source. Aims . We aim at characterizing the morphology of the radio emission and its dependence on the orbital motion, traced independently by near-infrared (NIR) interferometry of both the spectroscopic binary and the tertiary component comprising HD167971. Methods . We analyze 2006 and 2016 very long baseline interferometric data at C and X bands. We complement our analysis with a geometrical model of the wind-wind collision region and an astrometric description of the system. Results . We conﬁrm that the detected nonthermal radio emission is associated with the wind-wind collision region of the spectroscopic binary and the tertiary component in HD167971. The wind-wind collision region changes orientation in agreement with the orbital motion of the tertiary around the spectroscopic binary. The total intensity also changes between the two observing epochs in a way that is inversely proportional to the separation between the two components, with a negative-steep spectral index typical of an optically thin synchrotron emission possibly steepened by an inverse Compton cooling e ﬀ ect. The wind-wind collision bow-shock shape and its position with respect to the stars indicates that the wind momentum from the spectroscopic binary is stronger than that of the tertiary. Finally, the astrometric solution derived for the stellar system and the wind-wind collision region is consistent with independent Gaia data.


Introduction
Massive binaries (and higher degree multiple systems) made of O-type and/or Wolf-Rayet (WR) stars are known to produce stellar wind collisions and have shown their capability to accelerate particles up to relativistic speeds.The individual stellar winds collide, and form shocks that define the wind-wind collision region (WWCR).Within this region, electrons are accelerated to relativistic velocities, emitting synchrotron radiation which is detected as non-thermal radio emission.This emission is characterized by a negative spectral index 1 and a high brightness temperature (see e.g., Eichler & Usov 1993;Williams et al. 1990;White & Becker 1995;Dougherty et al. 1996;Dougherty & Williams 2000;Dougherty et al. 2003;Pittard & Dougherty 2006;Pittard et al. 2006;Blomme et al. 2007;Montes et al. 2009Montes et al. , 2015)).Most of the stellar non-thermal radio emitters are multiple 1 S ν ∝ ν α ; It is assumed that the flux density, S ν , is proportional to the frequency of observation, ν, scaling with a power-law defined by the spectral index α systems with at least one WR-star.There are only a few known non-thermal systems composed exclusively by O-stars (see e.g., Rauw 2004;Benaglia & Koribalski 2007).The study of these systems is necessary to understand the properties of their windwind interactions.During the last decade several individual studies have been done in stellar radio sources like 9 Sgr (Blomme & Volpi 2014), HD 168112 (De Becker et al. 2004;Blomme et al. 2005), or Cyg OB2 No. 9 (van Loo et al. 2008) among others, confirming that characterizing the radio emission in terms of the stellar parameters of the sources.O7.5III (Aa), O9.5III (Ab) and O8I (B)2 .Near-Infrared interferometric observations (obtained with AMBER, PIONIER and GRAVITY at the Very Large Telescope Interferometer -VLTIand reported by De Becker et al. 2012;Le Bouquin et al. 2017) demonstrated that the SB and T are gravitationally bound, with a projected separation that varies from 8 to 15 milliarcseconds (mas).The angular separation between the components of the SB is of the order of 0.1 mas, and cannot be spatially resolved directly with the VLTI nor with radio very long baseline interferometry (VLBI).Therefore, HD 167971 appears as a binary system (SB-T) in the near-infrared interferometric data.
Figure 1 presents the best-fit orbit of T around the SB reported by Le Bouquin et al. (2017).This solution is in agreement with previous estimates reported by the light-curve analyses of Blomme et al. (2007) and Ibanoglu et al. (2013).The VLA and ATCA radio light curves studied by Blomme et al. (2007) presented a periodicity consistent with the reported orbital period of the tertiary, and the emission shown a negative spectral index, thus confirming its non-thermal nature.These findings supported the hyphothesis that the non-thermal radio emission is produced in the wide orbit, probably, in the adiabatic wind colliding region between the SB and T. The difference of the light-curve profiles between 6 cm (5 GHz) and 20 cm (1.5 GHz), suggests that the emitting region is quite extended.
Additionally to the observed radio emission, De Becker et al. (2005) reported X-Ray (thermal) emission at energies of 2-4 keV that correspond to a high temperature plasma (2.3 -4.6 × 10 7 K), which is typical of pre-shock winds near their terminal velocity.This fact is in agreement with a wind-wind interaction scenario with strong adiabatic shocks, and a phase variability scaling with the inverse of the distance between the SB and T.These characteristics make HD 167971 a unique target to test the physics of wind-wind collisions in O-stars.
Here, we report new HD 167971 X-(3.5 cm; 8.4 GHz) and C-band (6.0 cm; 5 GHz) observations with the Very Long Baseline Array (VLBA) obtained in 2016, complemented with archival VLBA data, at the same frequencies, from 2006.Sect. 2 presents our observations and data reduction.In Sect. 3 the analysis of the radio emission and our results are presented.Then, in Sect. 4 we present the conclusions of our work.

Observations and data reduction
We performed VLBA observations of HD 167971 on August 13th, 2016 at C-and on August 14th, 2016 at X-band.We used a sustained data recording rate of 2 Gbit/s in two-bit sampling.Each frequency band was split into eight intermediate frequencies (IFs) of 32 MHz bandwidth each, for a total synthesized bandwidth of 256 MHz.Each IF was in turn split into 32 channels of 1.0 MHz bandwidth.The data were corrected from Earthorientation and ionospheric effects.The source J1751+0939 was used as bandpass calibrator and, the observations were phasereferenced to the source J1818-1108.Since this phase calibrator is not compact, all the fringe-fitting solutions were corrected from its structure.
The data were correlated at the NRAO data processor using an averaging time of 2s.We performed standard a-priori gain calibration using the measured gains and system temperatures of each antenna.This calibration, as well as the data inspection and editing, were done within the NRAO Astronomical Image Processing System (AIPS).No self-calibration was performed on the data since the peaks of emission were too faint.We also used AIPS to produce the images of both the reference source and target.Additionally, we analyzed archival VLBA observations performed on August 4th, 2006, half the orbital period of the system earlier.Then, a total synthesized bandwidth of 32 MHz was used, hence making those observations of lower sensitivity than those of 2016.The same calibrator source was used, and the fringe-fitting solutions were also corrected from the calibrator structure.Figure 2 shows the images of both epochs that were recovered by applying natural weighting to the data.Table 1 displays the main characteristics of the reconstructed maps.The synthesized beam, coordinates of the field center, peak-flux, and contour levels are also displayed in each one of the panels.

Radio Maps Analysis
The reported VLBA observations image the non-thermal radio emission of the WWCR. Figure 2 displays the radio maps, which show that the emission is resolved and consists of an elongated structure that changes depending on the observing epoch.The observing epochs correspond to half an orbital period, and show radio structures that have rotated by almost 90 • .Considering the orbital solution from Le Bouquin et al. ( 2017), and the radio images, we confirm that the long-axis of the radio emission -for every frequency and epoch-is perpendicular to the line that connects the SB with T. Our observations, thus, support the hypoth-esis that the observed radio emission corresponds to the windwind collision region between the SB and T. The following structural properties are highlighted: -The interacting point of the two winds could be characterized assuming that the WWCR is driven by the wind momenta ratio of the interacting stars (β= ṀT ν ∞, T / ṀSB ν ∞, SB ).In this case, a bowshock-like structure is formed at the shock front between the two winds.This bowshock has its tails pointing towards the star with the smaller wind momentum.In the 2006 epoch, such bow-shock (arc-like) feature is observed with small tails pointing towards T, implying that the wind momentum of the SB is larger.
-The fact that the bowshock is better defined in the 2006 map is consistent with the projected orbital position of T at the predicted orbital phase, φ.In the 2006 data3 (φ 2006 =-0.1), we are observing the parabolic shock profile with a small inclination angle4 (i=14.5 • ) from the plane of the sky, while in the 2016 epoch (φ 2016 =0.37 ; where T is behind the SB in our line-of-sight) we are observing the projected apex of the parabolic bowshock pointing towards us (i=-152.8• ) and, thus, its morphology resembles an elongated Gaussian instead of a parabolic arc (see Fig. 1 for a visual identification of the position of T at the two analyzed epochs).-In the X-band maps, we identified an additional compact structure to the bowshock.In the 2006 epoch it can be observed as a small, elongated appendix to the northwest of the central emission peak.In the 2016 map, a similar appendix is also observed to the northeast of the central peak, and the peak itself appears to be shifted from the center of the emission.Possible explanations for this include (i) the interaction of non-spherical/clumpy winds and, (ii) the possibility of the interaction of wind-momenta of similar magnitude that could disturb the profile from the parabolic bow-shock solution (see Fig. 2 in Canto et al. 1996).-We obtained spectral index maps of HD 167971 between 8.4 and 5 GHz (see Appendix A).The spectral index corresponds to non-thermal emission, with a value of α ∼-1.1 quite constant along the whole shocked region, suggesting that the Fermi acceleration mechanism is efficient.Dougherty et al. (2003) show that a synchrotron spectrum, without losses, goes as ∝ ν −0.5 .However, HD 167971 has a steeper spectral index, which suggest a possible attenuation of the synchrotron emission due to several physical mechanisms (see below).Blomme et al. (2007) determined, based on typical stellar parameters for O-stars, that the radio photosphere of optical depth unity at 5 GHz is ≤ 2000R (at 8.4 GHz, it is smaller), shorter than the minimum linear separation between the eclipsing binary and the tertiary component of the system (∼ 4000 R ).Therefore, the observed WWCR should be outside the radio photosphere, with a negligible free-free absorption, showing the intrinsic optically thin synchrotron emission at all orbital phases.It is predicted that at frequencies lower than 5 GHz (e.g., at 1.4 GHz), the radio photosphere is much more extended.Therefore, the radio emission of the WWCR is embedded in the radio photosphere for a large fraction of the orbit, resulting in spectral index and turnover frequency changes of the radio emission along the orbit, which are not detectable at 5/8.4 GHz.Under these assumptions, attenuation of the synchrotron emission at the observed frequencies due to the Razin effect (Dougherty et al. 2003) is also negligible.We suspect that the steepness of α is caused by inverse Compton cooling, which is more efficient for electrons with higher energies.According to the model presented by Pittard et al. (2006), inverse Compton decreases rapidly the energy of electrons close to the central volume of the WWCR, confining only a portion or relatively energetic electrons in layers near the shocks, which, in turn, results in a steeper spectrum at high frequencies.

Bow-shock Model
Since the projected bow-shock profile in the 2006 epoch does not exhibit long tails, despite the small inclination angle from the plane of the sky, we suspect that the wind-momenta ratio of the SB and T are quite similar.To estimate the separation of the WWCR and β, we used the algebraic solution of Canto et al. (1996) to define the shock profile of two interacting winds: where, in our system, R corresponds to the bow-shock distance from T, and D is the separation between the SB and T. The distance from T to the apex of the bowshock is known as the stand-off distance R 0 .The angles θ T and θ SB are the different half-opening angles between the line that connects the SB and T with the profile of the wind-wind shock front.The half-opening angles depend on β, as follows: Since the wind-wind shock profile is more homogeneous in the C-band, we used the 2006 map to model the bow-shock profile with the previous geometrical solution through the following method: -We mapped the position of the maximum of the source's brightness distribution for each column of pixels, perpendicular to the bow-shock structure, in the 2006 C-band image.
A prior rotation according to the orbital position, and image centering in a common reference frame were applied to the original map.We restricted our analysis to the structure with a brightness larger than 20% of the flux peak.To map the bow-shock profile, a Gaussian was fit to each non-zero element per column of pixels in the image.To estimate the errors in the position of the bow-shock profile, the standarddeviation of the Gaussian fit was used with a weight proportional to the number of pixel elements fitted.-Once the profile was estimated, it was compared with the projected canonical solution of the wind-wind bowshock obtained from a Monte-Carlo Markov-Chain (MCMC) model using emcee (Foreman-Mackey et al. 2013).Our MCMC model evaluates values of β between 0.3 and 0.9, smaller β values would produce a very close bowshock with long tails, while higher ones would produce a completely open bowshock.The model was run with ten independent chains with 1000 steps each one.For every iteration the distance between T and the shock front was adjusted for both the model and data, according to the following expression: Although the stand-off distance is variable in this procedure since the half-opening angle of the bowshock depends on the wind momenta ratio, we found the most probable model at a β and a stand-off distance that match better the half-opening angle of the extracted profile.The probability density function of the MCMC exhibits a clear peak at β=0.48±0.07(see Fig 3), which locates the shock front at a projected R 0 =0.41D±0.02Dfrom T, being D the separation between the SB and T across the orbital path.The β value corresponds, thus, to a half-opening angle θ T = 103 • ± 3 • of the bow-shock profile.
The used model is a solution for a radiative wind-wind collision shock.However, it is quite probable for the shock in HD 167971 to be adiabatic.Therefore, we compared the model solutions for adiabatic shocks described in Gayley (2009) (see also Pittard & Dawson 2018).The model for an adiabatic shock without mixing requires a β ∼ 0.39 (notice that this value is within 2σ of the reported β value for the radiative shock) to generate a half-opening angle of θ T ∼103 • .However, the situation changes when there is an adiabatic shock with mixing.In this case, the value of β depends on the ratio of the terminal velocity of the winds, u=ν P /ν S .The bigger the value of u is, smaller values of β are required to keep the required half-oppening angle (see Fig. B.1).In turn, the position of the shock relative to the two stellar components would be modified considerably.It is beyond the scope of this work to explore in detail these models.However, they should be taken into account for future studies.

Astrometric solution
Our VLBI radio observations provide us with the absolute astrometric position of the bowshock and our model allowed us to locate the projected position of the binary referenced to the radio emission.In Table 2, the positions of the emission's peaks of the radio maps are reported, together with the relative position of the SB from its respective radio maps (as computed from our model).We confirmed our 2016 astrometric solution by comparing the position of the system reported in the Gaia Data Release 2 archive (Gaia Collaboration et al. 2016a,b) 2).These new coordinates set the Gaia astrometric solution in between the line that connects the SB and T, confirming the astrometry derived from our VLTI binary fit and VLBA radio emission (see Fig. 4).

Conclusions
In this paper we have shown that the non-thermal radio emission observed in the C and X bands is in agreement with the windwind collision region between the spectroscopic binary and, the more distant, tertiary component in HD 167971.The morphology of the emission changes in accordance with the predicted orbital motion of the components as determined from NIR interferometric observations, with the semi-major axis perpendicular to the projected line that connects the SB and T across the orbital path.The total intensity also changes accordingly between the two observing epochs in a way inversely proportional to the separation of the stellar components in the long-period system, with a steep spectral index of α ∼ -1.1.Since the WWCR is outside the radio photosphere of the SB and T for all the orbital path, free-free absorption and the Razin effect are excluded as attenuating mechanisms of the synchrotron emission at the observed frequencies.In contrast inverse Compton cooling appears to play an important role in the steepness of α.Our model allowed us to locate the orbital components relative to the bow-shock position, finding an absolute astrometric solution for the system, which, in turn, was confirmed with an independent Gaia measurement.
This work presents the first step for understanding the physics of the observed synchrotron radio emission of HD 167971.It has shown the importance of VLBI observations to characterize the wind-wind interaction region in the system.Future multi-wavelength VLBI radio observations are necessary to fully characterize the emission across the different orbital phases.Monitoring the system would be important to constrain The reported errorbars have in consideration the error reported for the position of the binary system with respect to the radio emission derived from our model, and the intrinsic astrometric error of the radio maps.
b The reported errorbars include the Gaia uncertainty of the 2015.5 position and the uncertainties in the proper motions.
the stellar parameters of the components and to predict, through radiative transfer modelling, its evolution.Complementary observations with the high-resolution (R∼4000) mode of GRAV-ITY/VLTI are also envisioned to characterize possible infrared shock tracers (like Brγ) correlated with the observed WWCR at radio frequencies (see e.g., Sanchez-Bermudez et al. 2017;Gravity Collaboration et al. 2018).HD 167971 is a unique target that can bring us new insights about the synchrotron emission in multiple O-star systems.Understanding the properties of the stars, could help us to better characterize more distant targets in high-mass star forming regions (e.g., the stellar radio emitters in W51; Ginsburg et al. 2016), which in turn would provide important constraints on the radiation and chemical feedback of massive stars into the inter-stellar medium.

Fig. 1 .
Fig. 1.Left panel: Orbit of the SB-T system projected in the plane of the sky.Right panel: Orbit of the SB-T system in a plane orthogonal to the plane of the sky, formed by the observer's line-of-sight and the direction North.The orbital solution was obtained from Le Bouquin et al. (2017) and the mean orbital parameters are labelled on the image.The phases of T, at the two VLBI epochs reported in this work, are shown with orange and green dots in the orbital solution.The black-arrow shows the projected clock-wise motion of the tertiary component.The periastron and apastron are also displayed in the panels.The portion of the orbit in which T is in front of the SB, in our line-of-sight, is displayed in blue; while the portion of the orbit in which T is behind the SB is labelled in red.

Fig. 2 .
Fig. 2. Maps of the non-thermal radio emission observed in HD 167971.The radio-bands and epochs are shown in each one of the image panels.The synthesized beam, coordinates of the field center, peak-flux, and contour levels are also displayed in each one of the panels.

Fig. 3 .
Fig. 3. Probability density distribution of β.The values were obtained from the MCMC minimization of the geometrical modelling applied to the wind-wind bow-shock profile of HD 167971.

Fig. 4 .
Fig. 4. Astrometric position of the SB-T system and the C-band radio emission of the 2006 and 2016 epochs.The image is centered at R.A.= 274.5245616,Dec.= -12.2425868.The different stellar components and the radio emission are represented in black for the 2006 epoch, and in white for the 2016 one.The SB is represented by a cross and T with a dot.The projected orbital path is shown with a solid line.The 2006 Cband radio-emission is in colour scale contoured in black (representing 30, 50, 70, and 90% of the peak), while the 2016 one is only displayed in white contours.The thick-solid lines represent the 2D bow-shock profile from our geometrical model.The red diamond shows the Gaia astrometric measurement for the system in the 2016 epoch.