Open Access
Issue
A&A
Volume 624, April 2019
Article Number A55
Number of page(s) 7
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/201834659
Published online 09 April 2019

© J. Sanchez-Bermudez et al. 2019

Licence Creative Commons
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Open Access funding provided by Max Planck Society.

1. 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 nonthermal radio emission. This emission is characterized by a negative spectral index1 and a high brightness temperature (see e.g., Eichler & Usov 1993; Williams et al. 1990; White & Becker 1995; Dougherty et al. 1996, 2003; Dougherty & Williams 2000; Pittard & Dougherty 2006; Pittard et al. 2006; Blomme et al. 2007; Montes et al. 2009, 2015.) Most of the stellar nonthermal radio emitters are multiple systems with at least one WR star. There are only a few known nonthermal 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 wind-wind interactions. In the last decade several individual studies have been carried out on 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, characterizing the radio emission in terms of the stellar parameters of the sources, and of the wind-wind collision region.

HD 167971, at a distance of 1.8 kpc (De Becker et al. 2012), is one of the brightest synchrotron stellar radio emitters (Bieging et al. 1989). It is a hierarchical triple system (Leitherer et al. 1987), which consists of a spectroscopic binary (SB), with a period of ∼3.3 days, and a tertiary (T) component moving on a wider orbit, with a period of ∼21.4 years. The spectral types of the stars are O7.5III (Aa), O9.5III (Ab), and O8I (B)2. Near-infrared (NIR) interferometric observations obtained with AMBER, PIONIER and GRAVITY at the Very Large Telescope Interferometer (VLTI) and reported by De Becker et al. (2012) and Le Bouquin et al. (2017) demonstrated that the SB and T are gravitationally bound, with a projected separation that varies from 8 to 15 milliarcsec (mas). The angular separation between the components of the SB is of the order of 0.1 mas, and cannot be spatially resolved 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 showed a negative spectral index, thereby confirming its nonthermal nature. These findings supported the hyphothesis that the nonthermal 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 relatively extended.

thumbnail 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 north. The orbital solution was obtained from Le Bouquin et al. (2017) and the mean orbital parameters are labeled on the image. The phases of the 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 T. The periastron and apastron are also displayed in the panels. The portion of the orbit in which the T is in front of the SB, in our line-of-sight, is displayed in blue, while the portion of the orbit in which the T is behind the SB is labeled in red.

Open with DEXTER

In addition 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 × 107 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 the 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. Section 2 presents our observations and data reduction. In Sect. 3 the analysis of the radio emission and our results are presented. Subsequently, in Sect. 4 we present the conclusions of our work.

2. Observations and data reduction

We performed VLBA observations of HD 167971 on August 13, 2016, at C-band and on August 14, 2016 at X-band. We used a sustained data-recording rate of 2 Gbit s−1 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 for Earth-orientation and ionospheric effects. The source J1751+0939 was used as bandpass calibrator and, the observations were phase-referenced 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 2 s. 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 4, 2006, half the orbital period of the system earlier. A total synthesized bandwidth of 32 MHz was used, hence rendering the observations of lower sensitivity than those of 2016. The same calibrator source was used, and the fringe-fitting solutions were also corrected for 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.

thumbnail Fig. 2.

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

Open with DEXTER

Table 1.

VLBA observations.

3. Analysis and results

3.1. Radio map analysis

The reported VLBA observations image the nonthermal 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 the T. Our observations therefore support the hypothesis that the observed radio emission corresponds to the wind-wind collision region between the SB and the 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, a similar bow-shock (arc-like) feature is observed with small tails pointing towards the 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 the 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 the 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 the 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 nonspherical/clumpy winds or (ii) 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 nonthermal 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 suggests a possible attenuation of the synchrotron emission due to several physical mechanisms (see below).

Based on typical stellar parameters for O stars, Blomme et al. (2007) determined that the radio photosphere of optical depth unity at 5 GHz is ≤2000 R (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 or 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 rapidly decreases the energy of electrons close to the central volume of the WWCR, confining only a portion of relatively energetic electrons in layers near the shocks, which in turn results in a steeper spectrum at high frequencies.

3.2. 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:

(1)

where in our system R corresponds to the bow-shock distance from the T, and D is the separation between the SB and the T. The distance from the T to the apex of the bowshock is known as the stand-off distance R0. The angles θT and θSB are the different half-opening angles between the line that connects the SB and the T with the profile of the wind-wind shock front. The half-opening angles depend on β, as follows:

(2)

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 brightness distribution of the source 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 nonzero element per column of pixels in the image. To estimate the errors in the position of the bow-shock profile, the standard-deviation 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 of 1000 steps each. For every iteration, the distance between the T and the shock front was adjusted for both the model and data, according to the following expression.

(3)

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 the half-opening angle of the extracted profile better. 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 R0 = 0.41D ± 0.02D from the T, with D being the separation between the SB and the T across the orbital path. The β value therefore corresponds to a half-opening angle θT = 103 ° ±3° of the bow-shock profile.

thumbnail Fig. 3.

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

Open with DEXTER

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) and Pittard & Dawson (2018). The model for an adiabatic shock without mixing requires a β of ∼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 = νT/νSB. The bigger the value of u, the smaller the values of β required to keep the observed half-opening 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 these models in detail. However, they should be taken into account for future studies.

3.3. 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 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 2016a,b). The Gaia 2015.5 epoch reported a RA = 274.5245621 and Dec = −12.2425954 (with uncertainties around 0.1 mas) for the position of HD 167971. Correcting by the barycentric proper motions (PM = −0.54 ± 0.19 mas yr−1; PM = −1.25 ± 0.17 mas yr−1), we derived a mean Gaia astrometric position for the 2016.61 radio-epoch (see Table 2). These new coordinates set the Gaia astrometric solution in between the line that connects the SB and the T, confirming the astrometry derived from our VLTI binary fit and VLBA radio emission (see Fig. 4).

thumbnail 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 RA = 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, the T with a dot. The projected orbital path is shown with a solid line. The 2006 C-band 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.

Open with DEXTER

Table 2.

HD 167971 astrometry.

4. Conclusions

In this paper we have shown that the nonthermal radio emission observed in the C and X bands is in agreement with the wind-wind collision region between the SB and the more distant T 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 the 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 both 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 a first step towards understanding the physics of the observed synchrotron radio emission of HD 167971. It shows 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. A monitoring campaign is needed to constrain the stellar parameters of the components of HD 167971 and to predict, through radiative transfer modelling, its evolution. Complementary observations with the high-resolution (R ∼ 4000) mode of GRAVITY/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 2018). HD 167971 is a unique target that can bring us new insights into 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.


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 α.

2

The labels Aa, Ab and B follows the IAU component designation. However, in this work we would refer to the system Aa-Ab as the spectroscopic binary, SB, and the component B as the tertiary, T.

3

For the current work, we assume epoch T0 = 2454736.5 JD as reference for the orbital phases.

4

Inclination angle, i, is measured between the bowshock (rotational) symmetry axis to the plane of the sky. This angle increases clockwise in a plane orthogonal to the plane of the sky in direction of the observer’s line-of-sight.

Acknowledgments

We thank the anonymous referee for the very useful comments. J.S.B acknowledges support from the ESO Fellowship program. This work 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. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement no. [614922]. A.A., M.P.-T. acknowledge support from the Spanish MINECO through grant AYA2015-63939-C2-1P, cofunded with FEDER funds. J.C.G., and J.M.M. were partially supported by the Spanish MINECO project AYA2015-63939-C2-2-P. R.G.M. and J.S.B. acknowledge support from UNAM-PAPIIT Programme IN104319.

References

Appendix A: HD 167971 spectral index map

thumbnail Fig. A.1.

Spectral index of HD 167971 between X-band and C-band of the 2016 epoch is represented in colour. In contours (representing 30, 50, 70, and 90% of the peak) the image displays the structure of HD 167971 observed at the X-band. We note how the index across the whole structure is quite consistent around α = −1.1, confirming the nonthermal emission of the observed radio emission.

Open with DEXTER

Appendix B: Bowshock profiles

thumbnail Fig. B.1.

Value of β as a function of the half-opening angle (θT) of the bowshock profile for different types of wind-wind collision shocks according to Canto et al. (1996) (radiative) and Gayley (2009) (adiabatic). The labels on the frame describe each of the different shock types plotted. In the case of the adiabatic shock with mixing, three different cases are shown for different values of terminal velocity ratio, u. The vertical black solid line shows the corresponding β values θT = 103.4°, which is the best-fit value for the shock profile obtained from our 2006 C-band image.

Open with DEXTER

All Tables

Table 1.

VLBA observations.

Table 2.

HD 167971 astrometry.

All Figures

thumbnail 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 north. The orbital solution was obtained from Le Bouquin et al. (2017) and the mean orbital parameters are labeled on the image. The phases of the 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 T. The periastron and apastron are also displayed in the panels. The portion of the orbit in which the T is in front of the SB, in our line-of-sight, is displayed in blue, while the portion of the orbit in which the T is behind the SB is labeled in red.

Open with DEXTER
In the text
thumbnail Fig. 2.

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

Open with DEXTER
In the text
thumbnail Fig. 3.

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

Open with DEXTER
In the text
thumbnail 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 RA = 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, the T with a dot. The projected orbital path is shown with a solid line. The 2006 C-band 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.

Open with DEXTER
In the text
thumbnail Fig. A.1.

Spectral index of HD 167971 between X-band and C-band of the 2016 epoch is represented in colour. In contours (representing 30, 50, 70, and 90% of the peak) the image displays the structure of HD 167971 observed at the X-band. We note how the index across the whole structure is quite consistent around α = −1.1, confirming the nonthermal emission of the observed radio emission.

Open with DEXTER
In the text
thumbnail Fig. B.1.

Value of β as a function of the half-opening angle (θT) of the bowshock profile for different types of wind-wind collision shocks according to Canto et al. (1996) (radiative) and Gayley (2009) (adiabatic). The labels on the frame describe each of the different shock types plotted. In the case of the adiabatic shock with mixing, three different cases are shown for different values of terminal velocity ratio, u. The vertical black solid line shows the corresponding β values θT = 103.4°, which is the best-fit value for the shock profile obtained from our 2006 C-band image.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.