© ESO 2019

1 Introduction

The radio investigation of massive stars, that is, OB-type and Wolf-Rayet (WR) objects, revealed both thermal emission and non-thermal emission spectra (see e.g., Abbott et al. 1984, 1986; Montes et al. 2009). The thermal emission is free–free radiation from the optically thick stellar winds as described by Panagia & Felli (1975) and Wright & Barlow (1975), characterized by a spectral index α ~ 0.6 (for a frequency dependence of the flux density defined as S_ν ∝ ν^α). The non-thermal emission is synchrotron radiation produced by relativistic electrons in the presence of a magnetic field (White 1985), characterized by a negative spectral index. In this context, flatter spectra suggest a combination of unresolved thermal and non-thermal sources, the so-called composite spectra. The particle acceleration process responsible for the presence of relativistic electrons is believed to be diffusive shock acceleration (Bell 1978; Drury 1983), in the presence of the strong shocks produced by the collision of stellar winds in binary system (Eichler & Usov 1993; De Becker 2007). The subset of colliding-wind binaries known to accelerate particles is now referred to as the class of particle-accelerating colliding-wind binaries (PACWBs). The importance of binarity (or even higher multiplicity) for particle acceleration is emphasized in the catalogue of PACWB published by De Becker & Raucq (2013), where almost all systems show compelling evidence for multiplicity. To date, the detection of synchrotron radiation is the most efficient tracer of particle acceleration in massive binaries.

One striking fact, revealed by the inspection of the catalogue of PACWBs, is the spread of stellar and orbital parameters, suggesting the usual parameter space populated by a large fraction of massive binaries could be adequate for particle acceleration (De Becker et al. 2017). It is instructive to investigate the apparent limits of this parameter space, in order to achieve a better view of its extension. In this context, the lower period systems carry a significant information. First of all, the synchrotron emission is very likely to be significantly absorbed by the optically thick wind material in small size systems, preventing their identification as particle accelerators. Second, the efficient cooling of relativistic electrons by inverse Compton scattering constitutes a strong inhibition for electron acceleration, and this effect is enhanced in short period systems where the colliding wind region is closer to the stellar photospheres (see De Becker et al. 2017, for a discussion). A better view of the circumstances likely to allow for particle acceleration in shorter period systems is strongly needed to improve our understanding of the non-thermal physics in massive stars.

Among O-type stars, the shortest period system with known synchrotron radiation coming from a well-identified colliding-wind region is Cyg OB2 #8A, with a period of about 22 days (De Becker et al. 2004b; Blomme et al. 2010). Among WR systems, the shortest period referenced in the catalogue of PACWBs is WR 11, with a period of about 80 days (North et al. 2007). However, its non-thermal radio emitter status has recently been questioned by the recent radio investigation by Benaglia (2016), though it still may be a particle accelerator provided the Fermi γ-ray source reported by Pshirkov (2016) is physically related to WR 11.

The present study is dedicated to the next system in order of increasing period in the catalogue: WR 133. The paper is organized asfollows. In Sect. 2, we briefly describe the target and the reasons justifying a deeper investigation of its particle accelerator status. Section 3 details the processing and analysis of the data obtained with the Karl G. Jansky VeryLarge Array (JVLA). A general discussion is given in Sect. 4 and the conclusions are presented in Sect. 5.

2 The massive binary WR 133

WR 133 (HD 190918) consists of a WN5 star with an O9I companion. Underhill & Hill (1994) reported on an eccentric (e = 0.39) orbit with a period of 112.4 d, which is much shorter than most orbital periods identified for WR systems in the catalogue of PACWBs. The minimum masses and the size of the semi-major axis for both components of the system can easily be estimated on the basis of the results published by Underhill & Hill (1994) and are given in Table 1. On the basis of the calibration of stellar parameters for O9I stars given by Martins et al. (2005), the expected mass of the O-type component should be close to 30 M_⊙. Comparing this quantity to the minimum mass value given in Table 1 one estimates the inclination angle (i) to be on the order of 18°. This angle allows to estimate the absolute semi-major axis starting from the projected value published by Underhill & Hill (1994). In the absence of astrometric measurement of the orbit, this rough approach allows to achieve a reasonable view of its actual size.

The mass loss rate and the terminal velocity for the O9I star were retrieved from Muijres et al. (2012). For the WR component, we use the typical parameters for a WN5 star given by Crowther (2007). Considering the wind parameters given in Table 1, it is clear that the wind kinetic power of the WN wind dominates the wind collision that is located closer to the O-star. The wind temperature (T) is assumed to be 0.5 ×T_eff (Drew 1989), and will be used in Sect. 4.2.1.

Previous radio measurements presented by Montes et al. (2009) revealed a likely composite spectrum, that is, made of thermal and non-thermal emission without any strong and clear signature for a radio emission dominated by synchrotron radiation. This is not surprising as the small size of the system barely allows synchrotron photons to escape efficiently from the combined optically thick winds in the system, especially considering the WN component. The flux densities reported by Montes et al. (2009) are 0.31 ± 0.03 mJy and 0.57 ± 0.07 mJy at 8.4 and 23 GHz, respectively, with an upper limit of 0.41 mJy at 4.8 GHz. The spectral index is equal to 0.6, that is, the expected value for a purely thermal emission from a massive star wind (Panagia & Felli 1975; Wright & Barlow 1975). According to the ephemeris published by Underhill & Hill (1994), the observation date (6 May 2007) corresponds to orbital phase ϕ = 0.55, shortly after apastron passage. However, measurements at 4.8 and 8.4 GHz (0.38 ± 0.08 and 0.27 ± 0.03 mJy, respectively) performed on 31 May 1993 (ϕ = 0.27) suggest a spectral index of –0.61 ± 0.43, typical of many non-thermal sources. Though the error bar on the spectral index is large, the confidence interval still indicates at least a composite spectrum, and rejects a pure thermal emission. Measurements by Montes et al. (2009) are represented in Fig. 1 as a function of frequency. The different spectral nature roughly below and above 8.4 GHz could be interpreted in terms of composite spectrum, with the lower frequencies revealing some synchrotron emission, and the higher frequency domain dominated by thermal emission mainly from the WN stellar wind. One also has to caution that both epochs correspond to different orbital phases, and the difference may also be attributed to an orbital phase-locked variation.

On the other hand, the time interval of more than 20 yr between both observations raises the question of a potential longer term variability of the radio emission, independently of the expected variation along the 112 days orbit. The only point of comparison is the measurement at 8.4 GHz obtained at both epochs. The flux densities are equivalent within error bars, but they correspond to different orbital phases. One cannot therefore conclude on the long term behavior of the radio spectrum of WR 133 on the basis of these previous observations only. Moreover, one has also to envisage the possibility that the uncertainties on the ephemeris published by Underhill & Hill (1994) propagate significantly in the determination of the orbital phases of the observations. More than 25 yr elapsed since the reference time of the ephemeris at the time of our new JVLA observations, corresponding to more than 80 orbits of the system. One cannot thus reject the idea that the orbital phases are not correct. Without any more accurate determination of the ephemeris, one can only speculate on the potential impact of this error propagation. However, the spread of our new time series (see Table 2) is enough to sample the 112.4-d orbit, whatever the impact of this ephemeris error propagation. A revision of the orbital solution based on additional spectroscopic observations to be used altogether with the Underhill & Hill (1994) is thus necessary. In particular, a significantly much longer time basis of radial velocity measurements would allow to reduce significantly the error on the period, and reduce the uncertainties on the orbital phases of the JVLA observations as well.

Our motivation is to monitor WR 133 to connect the previous indication of non-thermal radio emission to its orbital phase-dependent behavior, as expected for an eccentric binary system. This would allow in principle to clarify the role of free–free absorption as a function of the stellar separation along the orbit, which is anticipated to be quite strong in this system (strong stellar winds, rather short period). A phase-locked variation would constitute a strong support for the colliding-wind origin of the synchrotron radiation in WR 133. Any deviation from this anticipated behavior would trigger new questions for our understanding of this system, and by extension open prospects for future investigations.

Fig. 1 Radio measurements published by Montes et al. (2009). Measurements at both epochs are overplotted.

3 Radio observations of WR 133

At our request, WR 133 was observed with the JVLA in configuration A at C and X bands, at various epochs (Obs. Id 14A-074 and 15A-023; see Table 2), respectively, in 2014 and 2015. The orbital phases of the observations determined on the basis of the ephemeris of Underhill & Hill (1994) are given in Table 2.

The observations of 2015 were made at the frequency bands centered at 5.5 GHz (band 1) and 9 GHz (band 2), along a 2 GHz bandwidth, while 2014 observations were made in band 1 only. The data were reduced with the Common Astronomy Software Applications package (CASA, McMullin et al. 2007). The source 3C295 was used as primary calibrator, with fluxes of 5.582 Jy at 5.488 GHz (C band) and 3.077 at 8.988 GHz (X band; Perley & Butler 2013). J2015 +3710 served as the phase calibrator. At each frequency band, the bandwidth comprised 16 spectral windows of 64 channels each. The time on source at each band and epoch was 16.5 min in 2015, and 5.2 min in 2014 in band 1. The theoretical image noise was 15 microJy beam⁻¹ in 2014 and 8 microJy beam⁻¹ in 2015.

Bad data were flagged, and the images were built using robust weightings¹. The final images show an rms between10 and 20 microJy beam⁻¹. WR 133 is detected in all images. The synthesized beams and the rms attained of each image are given in Table 2. The July image obtained at 9 GHz is displayed in Fig. 2. Two point sources are detected in the field, our target (WR 133) and an additional source northeast from WR 133. Since the synthesized beam is larger than the true dimensions of the sources in both bands, the size of the imaged sources is that of the synthesized beam. The flux density of WR 133 was estimated by fitting Gaussians, and the measured values are quoted in Table 2. For data analysis we also used the Miriad package (Sault et al. 1995) and the karma software (Gooch 1996). The spectral index for the 2014 and 2015 observations was determined using the equation α = log (S_ν,1∕S_ν,2)∕log (ν₁∕ν₂), and the error bar was estimated on the basis of the following equation: $$\Delta \alpha =\sqrt{{\left[\frac{1}{\ln ({\nu }_1/{\nu }_2)}\right]}^2\left[{\left(\frac{\Delta S_{\nu ,1}}{S_{\nu ,1}}\right)}^2+{\left(\frac{\Delta S_{\nu ,2}}{S_{\nu ,2}}\right)}^2\right]}.$$(1)

The spectral index measurements we determined here are thus inter-band values. Considering the signal-to-noise ratio on the data, any measurement of in-band spectral index would lack the required accuracy to carry any relevant information. In addition, no specific in-band trend different from the inter-band one is expected.

We also compared our new data to observations obtained earlier by Montes et al. (2009), for which the VLA interferometer bandwidth was 50 MHz. To do so, we selected the spectral window (spw) that fell near or within the same central frequencies observed by Montes et al. (2009), which were 4.8 and 8.4 GHz. For the 5.5 GHz frequency data, we chose spw = 1–3 (4616–5000 MHz, tagged subb1) near 4.8 GHz. For the 9 GHz band, we chose spw = 2–4 (8244–8628 MHz, tagged subb2) near 8.4 GHz. The images were obtained using CASA. Figure 2 shows the image for the Jul. 2015 epoch data, that showed lower noise compared to the Sep. 2015 ones. The flux densities of WR 133 were estimated using the MIRIAD package. The flux of WR 133, for Jul. 2015, at subb2 is 0.27 ± 0.03 mJy, in accordance with the corresponding values reported by Montes et al. (2009). The measurement at subb1 band yielded a target flux of 0.21 ± 0.02 mJy. Although this result agrees with the flux upper limit derived from 2007 data presented by Montes et al. (2009, and also with the upper limit presented by Abbott et al. 1986), it is remarkably lower than the flux reported from the 1993 data also by Montes et al. (2009).

We also note that we reprocessed the 1993 data in order to achieve the full confirmation of the non-thermal emitter status of WR 133 at that epoch, as this constitutes a pivotal element in our understanding and interpretation of the radio properties of this system. Our analysis led to flux densities of 0.38 and 0.28 mJy, respectively, at 4.8 and 8.4 GHz. These values are in perfect agreement with the results published by Montes et al. (2009), and plotted in Fig. 1. The important information that will be considered in this study is the indication for non-thermal radio emission previously reported by Montes et al. (2009) for WR 133, on the basis of the 1993 measurements, which appears as a highly valid and consolidated fact.

Fig. 2 Radio image showing the two point sources detected in the JVLA data at 9 GHz, in July 2015, WR 133 (bottom) and the NEsource (top). The black ellipse at the bottom left corner illustrates the beam size. The contour levels are −0.025, 0.025, 0.08, 0.16, 0.24 and 0.32 mJy beam−1.

4 Discussion

4.1 The orbit of WR133

The radio measurements should be adequately interpreted in the context of the orbital configuration. Figure 3 represents the orbit of WR 133 based on the parameters given by Underhill & Hill (1994) and shows the relative position of the stars at the phases of our JVLA observations. For the sake of completeness, we also show the stellar positions at the epoch of the Montes et al. (2009) observations. One should clarify that the reference time for periastron passage determined by Underhill & Hill (1994) is located in 1987. As a result the 1993 measurements (where the non-thermal emission is measured, at ϕ = 0.27) is affected by the lowest uncertainty on the orbital phase determination, while our most recent measurements can in principle be more affected. Our JVLA measurements were indeed obtained about 25 yr after the reference time, translating into about 80 full elapsed orbital periods (whose error bar is about 0.2 d according to Underhill & Hill 1994). In the worst (an unrealistic) situation it may lead to an uncertainty on the orbital phase determination of several percent. Anyway, as stated in Sect. 2, these uncertainties do not question the nature of our sampling of the 112.4-d orbit with our new measurements. On the basis of the values given in Table 1, the absolute median separation is a = a_O + a_WR = 104.9 R_⊙, with extreme values a (1 − e) and a (1 + e) respectively at periastron and apastron.

As a first approach, let us make use of the wind parameters quoted in Table 1. We notethat the uncertainties on the terminal velocities should be on the order of 5–10% at most, if one considers an uncertainty on the stellar classification of one subtype or one luminosity class (according to the references quoted in the table). The same is true for the O-star mass loss rate. For the WN star, the uncertainty on the mass loss rate may reach a few tens of percent. Such uncertainties are inherent to usual stellar parameters of stellar winds of massive stars, with no significant impact on our conclusions. These parameters allow in particular to compute the wind momentum rate ratio (η) according to the following relation: $\eta ={(\dot{M}{V}_{\infty })}_{\text{O}}/{(\dot{M}{V}_{\infty })}_{\text{WR}}$, leading to the value η = 0.13. The position of the wind collision region (WCR) measured from the center of mass of the two stars can then be calculated: $r_{\text{WR}}=a{r}_{\text{rel}}/(1+\sqrt{\eta })$ and $r_{\text{O}}=a{r}_{\text{rel}}\sqrt{\eta }/(1+\sqrt{\eta })$, where r_rel = (1 − e)∕(1 + e cosυ) and υ being the true anomaly. These positions are also represented in the plane of the orbit in Fig. 3.

The size of the system is rather small, that is, less than an astronomical unit, placing the wind collision region in the wind acceleration zone of the O-type star. A kinematic adjustment should thus be considered, taking into account the “actual” pre-shock velocity instead of assuming terminal values, especially close to periastron passage. Such a kinematic adjustment leads the WCR toget closer to the O-star. Unstable solutions are expected as the distance between the O-star and the WCR calculated following this simple kinematic approach can even become smaller than the expected stellar radius, that is, about 23 R_⊙ (Martins et al. 2005). Beside, the estimated position of the WCR is sensitive to the adopted wind parameters, in particular those of the WR component which are quite uncertain. Even lowering its terminal velocity down to 1000 km s⁻¹ does not solve this issue. On the other hand, a higher mass loss rate for the WR could not be completely rejected (see for instance WN wind parameters given by Hamann & Koesterke 1998), widening further the orbital phase interval around periastron where the WR wind seems at first sight to crash onto the companion’s surface. One should however note that this approach does not take into account radiative effects such as radiative inhibition (Stevens & Pollock 1994) and sudden radiative breaking (Gayley et al. 1997), likely to prevent the WR wind to crash onto the O star surface. In addition, the simple description provided here does not take into account the orbital asymmetry due to the orbital motion that distorts the WCR. The full radiative and hydrodynamic treatment of these issues is beyond the scope of the present work. The values quoted above for the position of the wind–wind interaction region should thus be considered as reasonable guess values adequate for the present qualitative discussion.

Fig. 3 Illustration of the orbit of WR 133. The open hexagons and triangles represent the position of the O-star, respectively during two different orbits. The filled symbol stands for the expected position of the wind collision region. The long-dashed and dotted lines represent the stellar separation, in 2014 and 2015 respectively. The orbital phases of the Montes et al. (2009) observations are illustrated by the two open star symbols.

4.2 Radio emission from WR 133

The flux densities quoted in Table 2 suggest a constant radio emission along the 112.4-d orbit. The spectrum seems to be purely thermal. These results are at odd with those published by Montes et al. (2009) from VLA data, presenting the signature of a non-thermal emission component (see Sect. 2). We note that the observations presented here were accomplished with the (upgraded) Jansky Very Large Array. The flux models and scale of the primary calibrator used in the present investigation have been widely extended in frequency and rebuilt (Perley & Butler 2013). In the light of these changes, a meaningful comparison between our values and the ones in the literature is doubtful. The fluxes derived prior to 2010 corresponded to a bandwidth of 50 MHz, while the ones presented here were taken along a 2-GHz band. In addition, one also has to take into account the expected strong free–free absorption in that system, especially considering the dense WR wind likely to attenuate dramatically any radio emission produced close to the wind-wind interaction region. For these reasons, the interpretation of the radio behavior is not straightforward and it is worth considering two distinct potential interpretation frameworks, in other words the binary system and the triple system scenarios. These scenarios are critically discussed below, and summarized in Table 3.

4.2.1 The binary system scenario

In this scenario, the thermal and non-thermal emissions reported by Montes et al. (2009) and by us are produced in the WCR in the 112.4-d orbit. The rather small size of the orbit constitutes a severe issue considering the significant impact of free-free absorption by the stellar winds material. As a rough estimate of the impact of the free-free absorption by the stellar winds in the system, one can use the formalism developed for instance by Wright & Barlow (1975) and Leitherer et al. (1995) to estimate the typical radius of the radio photosphere of a given stellar wind. The radius of the radio photosphere is defined for an optical depth equal to one, and can be computed following the approach already applied in previous studies (e.g., De Becker et al. 2004a; Mahy et al. 2012).

Using the parameters given in Table 1, we derived curves for the radio photosphere radius as a function of wavelength for both stars in the system (see Fig. 4). The radio photosphere radii at 5.5 and 9 GHz are significantly larger than the size of the orbit, confirming the a priori idea that a substantial amount of the radio photons produced atthe position of the WCR should be absorbed. This is especially true for the WN5 curve, lending support to the idea that the WR stellar wind dominates the free–free absorption in this system. We note that within the uncertainties on the stellar wind parameters, in particular for the WN star (see Sect. 4.1), predictions on the radio photosphere radii would still be too much large to admit a potential escape of a measurable fraction of the putative synchrotron emission from the colliding-wind region in the 112.4-d orbit. For instance, lowering the mass loss rate of the WN star by a factor 2 would lead to a drop of the WN curve plotted in Fig. 4 by a factor of about 1.6 only. A much more severe change of the mass loss rate would be required, and as stated previously this is very unlikely that such a change in mass loss rate occurred between 1993 and the epoch of our new observations. One should keep in mind that the optical depth considered here is radial (i.e., centered on the star), and not measured along the line of sight toward the synchrotron emission region. The measurement of the line-of-sight optical depth would require a much more sophisticated treatment based on detailed hydrodynamic and radiative codes, out of the scope of the present study. One has however to mention the case of another system where the free–free absorption is a priori strong enough to prevent the detection of synchrotron emission: Cyg OB2 #8A. The latter system presents a clear non-thermal radio signature phase-locked with the 22-days orbital period (De Becker et al. 2004b; Blomme et al. 2010), providing evidence that the radio photosphere extension beyond the orbit size is certainly a pessimistic criterion. However, the comparison of the numbers plotted in Fig. 4 to the size of the orbit (at most ~105 R_⊙) clearly emphasizes the crucial impact of free–free absorption in this system, which is undoubtedly the dominant turn-over process preventing any obvious identification of synchrotron radiation from WR 133. In addition, the free–free absorption of any putative synchrotron emission component will be dependent on the orientation of the system. We are lacking the adequate information to achieve a detailed three-dimensional view of the orbit (i.e., the longitude of ascending node), but we estimate the inclination of the system to be on the order of 18° (see Sect. 2). It is not sure that such an inclination is low enough to favor the potential detection of a fraction of the synchrotron emission.

The most obvious source of variation in an eccentric system is the changing stellar separation along the orbit. One could apriori expect the detection of a non-thermal emission close to apastron, where a larger fraction of the synchrotron emission region may emerge from the opaque WR wind. However, even at orbital phase ϕ = 0.46 (i.e., the closest to apastron in our time series), we do not find any indication of non-thermal emission. Our sampling of the radio emission in two distinct orbits suggest a pure thermal emission along the orbit, in agreement with the expectation from the extension of the radio photosphere. As a result, we do not see how to reconcile our new measurements with the composite or non-thermal emission reported by Montes et al. (2009), in the context of the binary scenario with the 112.4-d period.

Fig. 4 Radio photosphere radius as a function of wavelength for both star winds in the WR 133 system. Vertical lines are located at wavelengths corresponding to our two observation bands.

4.2.2 The triple system scenario

The catalogue of PACWBs (De Becker & Raucq 2013) includes several hierarchical triple systems, with the measured non-thermal radio emission coming from the WCR of the outer, wider orbit. It is worth to investigate how such a scenario may be compatible with the results of radio observations of WR 133, though no hint for a third object exists to date at other wavebands. In such a scenario, one may distinguish two orbits: the so-called inner orbit as illustrated in Fig. 3, and the hypothetical outer orbit. Both are characterized by their own period and eccentricity, and are associated to a WCR. If the size of the inner orbit is too small to warrant a detection of synchrotron emission from its WCR (see the above discussion on the free-free absorption), the wind collision in the outer orbit may be located far enough from the inner parts of the WR wind to be significantly unveiled in a significant part of the orbit. In this scenario, the lack of detection of non-thermal radio emission in 2015 may be explained by a system close enough to the outer periastron passage, when the outer WCR reaches deeper layers of the opaque stellar winds of the inner system.

The existence of such a hypothetical triple system raises the question of its dynamical stability. To address this question, we made use of the stability criterion proposed by Eggleton & Kiseleva (1995). This criterion expresses the critical, minimal orbital period ratio (outer/inner) allowing for the dynamical stability of a triple system, as a function of the mass ratio of the inner and outer systems, and of their respective eccentricities. Using Eqs. (1)–(3) in Eggleton & Kiseleva (1995), we computed the critical period of the outer orbit for a range of mass ratio of the outer system (q_out), assuming different arbitrary values for the eccentricity of the outer orbit (e_out = 0.25, 0.50 and 0.75). The results are plotted in Fig. 5. For q_out, we considered somewhat conservative values. The minimum boundary is 1, which is quite optimistic. A third object as massive as the WN5 + O9Ipair would be bright enough to have been easily revealed by previous spectroscopic observations. For the maximum boundary, a value of 4 would be adequate for instance if the third component is an early B object (m₃ ~ 11–12 M_⊙). Such a star would still produce a stellar wind able to sustain a significant wind-wind interaction for non-thermal physics, and may be faint enough as compared to the two other components to remain undetected so far. Intermediate values of q_out could stand for late-type O-stars, especially on the main-sequence for high values of the mass ratio. One sees that outer periods in such a hypothetical triple system should be of at least several hundreds of days to warrant dynamical stability. Provided the outer orbit is long enough there is a priori no physical reason to consider such a triple system is unlikely.

The outer orbit must be wide enough to warrant a significant detection of non-thermal emission, but in addition it should be significantly eccentric to explain the very different behaviors reported between the epochs of the Montes et al. (2009) measurements (i.e., year 1993) and our observations. Close to periastron, the wind–wind interaction region (where electrons are certainly accelerated) is located closer to the inner pair. The synchrotron emission region is thus more deeply embedded into the dense WN5 + O9I winds, leading to a more severe free-free absorption of the synchrotron spectrum. At orbital phases far enough from periastron, the wind-wind interaction will be mainly located outside the radio photosphereand some composite radio emission (thermal + non-thermal) may be detected, as was the case for the 1993 measurements.

Finally, let us emphasize that any direct indication of the presence of a third component in the system is still lacking to date. At this stage, the triple system scenario is envisaged for the sake of completeness, and it may be supported or rejected according to future observations of WR 133. One has however to remind that the late detection of a third component orbiting a previously known massive binary system is not a scarce event (some of these examples are included in the catalogue of PACWBs, De Becker & Raucq 2013). For a quite eccentric orbit, the amplitude of the radial velocity variations may be quite low along a large part of the orbit, with measurable changes while approaching periastron passage. A low inclination of the outer orbit would also inhibit the detection of the putative third object. Let’s note also that the time sampling of the radial velocity measurements used by Underhill & Hill (1994) to compute the orbital solution is quite sparse, for a time basisof about six years. Such a sampling leaves significant room for a still undetected companion on a wide eccentric orbit. This is especially true considering the difficulty to identify the spectroscopic signature of companions in multiple systems harboring a Wolf-Rayet star. In this context, especially considering the likely long period of the putative third component, long baseline interferometry constitutes a very adequate technique, as demonstrated for instance by Le Bouquin et al. (2017) for a sample O-type binaries. This would allow us to search for an astrometric companion, and to monitor the expected orbital motion over long time scales (i.e., up to several years).

Fig. 5 Estimate of the critical minimum outer period of WR 133 allowing for a dynamical stability of a hierarchical triple system. The critical period (Pout) is plotted as function of the outer mass ratio (qout), assuming three values for the outer orbit eccentricity.

4.3 On the identification of synchrotron radio emitters

The switch between a thermal and a composite radio spectrum is especially relevant for the discussion of the detectability of synchrotron radiation in massive colliding-wind binaries. At best, a non-thermal/composite spectrum is identified at one epoch. Any observation of that system when the stellar separation is short enough to bury deeper the synchrotron emission region in the combined stellar winds will only reveal a purely thermal spectrum. In the hypothetical triple system scenario, the orbital phase interval where non-thermal emission should be detectable is wider, but observations reported here show that it is not detected in July 2014 and September 2015. This clearly indicates how critical observational biases affect our capability to detect non-thermal radio emission in these systems, and therefore to attribute a particle accelerator status.

The issue of observational biases is critical for the identification of PACWBs, as discussed by De Becker et al. (2017). Several systems already observed with radio observatories and displaying a thermal spectrum may reveal a non-thermal emission component at other orbital phases. On top of that, the inclination of the system and the longitude of ascending node will dictate which stellar wind contributes the most to the line-of-sight absorbing material at any orbital phase. This orientation effect is especially sensitive to dense and opaque winds of Wolf-Rayet stars, along with early supergiant winds with enhanced mass loss. This has two essential consequences. First of all, the quest for new particle accelerators to be included in the PACWB catalogue should proceed through multiepoch campaigns, improving the probability to observe a system at an orbital phase where the synchrotron emission may be detected. For systems with a well-characterized orbit, radio observations may focus on a few selected orbital phases. Second, this observational bias strongly suggests that the fraction of PACWBs among massive binaries is significantly underestimated by the present census.

4.4 Radio emission from a nearby source

The analysis of the radio images revealed the presence of a point source at about 10 arcsec northeast from WR 133, at coordinates [20:05:57.70; 35:47:26.05]. The closest counterpart in the 2MASS catalogue (Skrutskie et al. 2006) is 2MASS J20055779+3547287, at about 3 arcsec of the NE source. According to Mason et al. (1998), a few visual companions close to WR 133 are identified in the Washington Double Star catalogue (WDS), at angular distances ranging between 6 and 36 arcsec. None of these visual companions appear at a position angle of 30–40° with a separation of about 10 arcsec².

Following the same procedure as for WR 133, we determined the flux density in both bands and we derived the spectral index (see Table 4). This source seems to be constant and presents a non-thermal spectrum at both epochs, strongly suggesting the emission is significantly made of synchrotron radiation.

Assuming the source is of Galactic origin, one can only speculate on its nature. One may mention a couple of scenarios which require additional information to be checked consistently. The source could either be a massive binary presenting the same behavior as PACWBs (e.g., De Becker & Raucq 2013), or a young stellar object presenting non-thermal emission associated with jets interacting with the interstellar medium (e.g., Rodríguez-Kamenetzky et al. 2017). Negative spectral indices are also reminiscent of supernova remnants (SNRs) known as emblematic Galactic synchrotron sources (e.g., Dubner & Giacani 2015). Even though most SNRs are characterized by spectral indices close to − 0.5, some of them display values up to − 0.3. Such flatter spectra are sometimes measured for older SNRs, even though the physical justification for such a value is not straightforward (Reynolds et al. 2012). However, if known SNRs quoted in Green’s catalogue were moved to distances such that their flux density (actually, of the order of 1–10 Jy) would be as low as a fraction of mJy, one would have to push them beyond the limit of the Milky Way. This does not argue in favor of a Galactic SNR for this unidentified source.

5 Summary and conclusions

We investigated the radio emission from the WN + O binary WR 133 on the basis of dedicated JVLA observations at 5.5 and 9 GHz, at various epochs spread over two different orbits. Our measurements reveal a thermal spectrum, in apparent contradiction with the detection of a non-thermal/composite spectrum reported by Montes et al. (2009).

In the context of a simple binary system scenario, our monitoring of the 112.4-d orbit rejects the idea that a non-thermal emission component could be detected. This can easily be explained by the strong free–free absorption due to the dense stellar winds, constituting a severe obstacle for the emergence of a measurable non-thermal radio contribution from the wind collision region. However, this does not provide any satisfactory interpretation context for the non-thermal emission reported previously. Alternatively, one may think about a triple hierarchical system scenario. In this interpretation framework, the hypothetical third object would move along on a wide outer orbit. The outer wind collision region would be much less affected by the free–free absorption, opening the possibility that a non-thermal emission component may be detectable at some epochs provided the outer orbit is sufficiently eccentric. According to simple dynamical stability considerations, the outer orbital period should be at least several hundred days. The main drawback of the triple system interpretation is the present lack of any hint for a third object in this system, though (i) previous investigations of other systems illustrated the difficulty to identify high multiplicity systems, and (ii) previous spectroscopic studies of WR 133 are not sufficient to completely reject the triple system scenario. The search for a putative third object on a wider orbit would require either dedicated spectroscopic measurements over a long time basis, or more likely visible/infrared interferometric observations using long baseline interferometers to search for an astrometric companion.

Beside the question of the multiplicity (binary or triple) of WR 133, the change in the spectral shape (with a composite spectrum at one epoch) provides evidence for the existence of a synchrotron source, despite the strong free–free absorption at work in this system. The epoch-dependent detection of the non-thermal emission constitutes a compelling illustration of the observational biases that severely affect the detectability of synchrotron sources in colliding-wind massive binaries. The depth of the location of the synchrotron emission region in the combined stellar winds and orientation effects can represent significant obstacles for the detection. A system may present the properties of a purely thermal emitter in a significant part of its orbit, with only hints for synchrotron emission in another orbital phase interval. One can thus anticipate that the actual population of particle accelerators among colliding-wind binaries is certainly underestimated. Observation strategies aimed at unveiling new particle accelerators should take this issue into account notably through an adequate selection of orbital phases (when feasible), or at least by organizing multiepoch observations to sample the orbit.

Finally, our JVLA observations revealed the existence of non-thermal point source, about 10 arcsec northeast from the position of WR 133. The source seems constant, as far as our time sampling allows to say. The nature of the object is still undetermined and it may be related to a Galactic stellar source such as another colliding-wind binary or a young stellar object.

Acknowledgements

The authors want to express their gratitude to the NRAO staff for the scheduling of the VLA observations, to Drew Medlin, Rick Perley and the NRAO helpdesk. Our acknowledgements go also to Dr. Juan-Ramon Sanchez-Sutil, for his previous help in the preparing of the telescope time proposals and observations and to C. Hales for fruitful discussions on data reduction techniques. The authors also thank the anonymous referee for his/her constructive report that allows to significantly strengthen our conclusions. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. M.D.B. acknowledges the financial support from a ULiège “Crédit classique (DARA)”. The SIMBAD database was used for the bibliography. This research has made use of the Washington Double Star Catalog maintained at the U.S. Naval Observatory.

References

All Tables

All Figures

	Fig. 1 Radio measurements published by Montes et al. (2009). Measurements at both epochs are overplotted.
In the text

	Fig. 2 Radio image showing the two point sources detected in the JVLA data at 9 GHz, in July 2015, WR 133 (bottom) and the NEsource (top). The black ellipse at the bottom left corner illustrates the beam size. The contour levels are −0.025, 0.025, 0.08, 0.16, 0.24 and 0.32 mJy beam−1.
In the text

Fig. 3 Illustration of the orbit of WR 133. The open hexagons and triangles represent the position of the O-star, respectively during two different orbits. The filled symbol stands for the expected position of the wind collision region. The long-dashed and dotted lines represent the stellar separation, in 2014 and 2015 respectively. The orbital phases of the Montes et al. (2009) observations are illustrated by the two open star symbols.

In the text

	Fig. 4 Radio photosphere radius as a function of wavelength for both star winds in the WR 133 system. Vertical lines are located at wavelengths corresponding to our two observation bands.
In the text

	Fig. 5 Estimate of the critical minimum outer period of WR 133 allowing for a dynamical stability of a hierarchical triple system. The critical period (Pout) is plotted as function of the outer mass ratio (qout), assuming three values for the outer orbit eccentricity.
In the text