© The Authors 2024

1 Introduction

The first survey of radio emission in hot stars was done by Bieging et al. (1989). Among these stars, a good number of thermal radio emitters were detected, where the emission is due to free-free processes in the ionized stellar wind. The measured radio flux of these stars can be used to determine the mass-loss rate (Wright & Barlow 1975; Panagia & Felli 1975).

However, the survey also unveiled a number of non-thermal radio emitters, as recognized by the spectral index of their emission¹. They also have an excess flux compared to that expected from the free-free wind emission, a high brightness temperature, and a variable flux. Since that first survey, substantially more non-thermal radio emitters have been found (for a review and a catalogue, see e.g. De Becker 2007; De Becker & Raucq 2013).

For Wolf-Rayet stars, it quickly became clear that the non-thermal emission is mostly due to colliding-wind binaries (Dougherty & Williams 2000). In a massive-star binary system, the two stellar winds collide and shocks form on either side of the collision region. Around these shocks, a fraction of electrons can be accelerated by the Fermi mechanism to relativistic speeds (Bell 1978; Reitberger et al. 2014; Pittard et al. 2021). By spiralling around in the magnetic field, these electrons then emit synchrotron radiation, which we detect as non-thermal emission (Eichler & Usov 1993).

The situation for the O-type non-thermal emitters was not so clear. A number of sources were known to be binaries: HD 15558 (De Becker et al. 2006, and references therein), HD 167971 (a triple system – Leitherer et al. 1987, and references therein). Over time, most of the other sources classified by Bieging et al. (1989) as ‘definite’ or ‘probable’ non-thermal emitters were also found to be binaries: 9 Sgr (Rauw et al. 2016), Cyg OB2 #8A (De Becker et al. 2004b), and Cyg OB2 #9 (Nazé et al. 2008).

HD 168112 was also listed as a non-thermal radio emitter by Bieging et al. (1989), but for a long time there was only indirect evidence of it being a binary. De Becker et al. (2004a) confirmed the non-thermal radio emission and showed that the X-ray emission is variable, but their optical spectra did not show any indication of binarity. Blomme et al. (2005) found periodic behaviour in the archive radio data, which does suggest that it is a binary system. A subsequent X-ray Multi-Mirror Mission (XMM-Newton) observation by De Becker (2015) showed the X-ray emission to be overluminous, which can be attributed to heated material in the colliding-wind region. Sana et al. (2014) obtained optical interferometry observations with the Precision Integrated-Optics Near-infrared Imaging ExpeRiment (PIONIER) instrument at the Very Large Telescope Interferometer (VLTI) and found a companion at an angular distance of 3.33 ± 0.17 mas with a brightness difference in the H-band of 0.17 ± 0.19 mag.

Early attempts by Rauw et al. (2005) and Chini et al. (2012) to find spectroscopic evidence for binarity were not successful. However, Maíz Apellániz et al. (2019b) succeeded in assigning two spectral types to the components of HD 168112: O5 IV(f) + O6: IV: or O4.5 III(f) + O5.5 IV((f)) for the LiLiMaRlin data (Library of Libraries of Massive-Star High-Resolution Spectra, Maíz Apellániz et al. 2019a).

Finally, Putkuri et al. (2023) published a spectroscopic orbit determination, providing convincing evidence that HD 168112 is indeed a binary. They found it to be a highly eccentric system (with eccentricity e = 0.743 ± 0.005) and a period of P = 513.52 ± 0.01 days. The masses of the two components are very similar and the minimum value they derived is ~27.1 M_⊙ for the A component and ~24.6 M_⊙ for the B component.

Additionally, radio observations by De Becker et al. (2024) resolved the colliding-wind region of HD 168112. They used the European Very Long Baseline Interferometry Network (EVN) to observe at 1.6 GHz (18 cm). The colliding-wind region is clearly asymmetric and its high flux (integrated 1.7 ± 0.4 mJy) shows it to be non-thermal in origin, thereby providing further evidence for the colliding-wind nature of the HD 168112 binary.

Colliding-wind binaries are highly interesting because they provide the opportunity to study the acceleration of particles around shocks. This process occurs in a number of astrophysical environments, such as interplanetary shocks and supernova remnants. The colliding-wind binaries provide an environment that is quite different from the others regarding parameters such as density, magnetic field, ambient radiation field, and shock speed. They are also relevant for the mass-loss rate determinations in single stars, as they can provide an independent determination of the effect of clumping and porosity (e.g. Pittard 2007). This in turn will help determine the stellar mass-loss rates, which are very sensitive to clumping and porosity (e.g. Puls et al. 2008).

In this paper, we present independent spectroscopic observations of HD 168112, covering a time span of 13 yr. From these data, we derived an orbital solution. The combination of the spectroscopic data and the optical interferometry results of Sana et al. (2014) allowed us to further constrain the orbital parameters. We also determined the stellar parameters of both components. Furthermore, in this work, we analyse the radio observations we obtained near the 2013 periastron passage of HD 168112. These observations are part of a coordinated XMM-Newton and VLA project. The XMM-Newton data are analysed in a separate paper (Rauw et al. 2024).

In Sec. 2, we present both the spectroscopic data and radio data we obtained. Section 3 analyses the spectroscopic data, determining the orbital parameters and the stellar parameters of both components of the binary system. In Sec. 4 we analyse the new radio data. Section 5 combines our results with data from the literature to further determine the orbital and stellar parameters. Section 6 presents our conclusions.

2 Data

2.1 Optical spectroscopy

HD 168112 was monitored during 13 yr with the HERMES (High Efficiency and Resolution Mercator Echelle Spectrograph, Raskin et al. 2011) spectrograph on the 1.2 m Mercator telescope at the Roque de los Muchachos Observatory (La Palma, Canary Islands, Spain). We obtained 97 usable² HERMES spectra between June 2009 and September 2021. Our HERMES spectra are part of another observing programme than those used by Putkuri et al. (2023), and therefore provide an independent set of data.

The HERMES instrument has a resolving power R ≈ 85 000 and covers the wavelength range 3770–9000 Å. The data were reduced using the standard HERMES pipeline³. Spectral orders were extracted and flat-fielding and wavelength calibration (based on Th–Ar lamp spectra) was applied. Cosmic rays were removed, the orders were merged, and the barycentric velocity correction was applied.

An additional 16 spectra were obtained with the refurbished HEROS (Heidelberg Extended Range Optical Spectrograph) instrument on the 1.2 m Telescopio Internacional de Guanajuato Robótico Espectroscópico (TIGRE, Schmitt et al. 2014; González-Pérez et al. 2022), located at La Luz observatory, Mexico. The resolving power is 20 000, covering the wavelength range 3800–8800 Å. Data reduced with version 3.1 of the pipeline (Mittag et al. 2010) were downloaded from the TIGRE archive⁴. We applied the barycentric correction to the spectra. The observing log of both the HERMES and TIGRE observations is given in Table A.1.

2.2 Radio data

We obtained radio observations of HD 168112 near the periastron passage, using the Karl G. Jansky Very Large Array (VLA) of the National Radio Astronomy Observatory (NRAO). The VLA project SX222004 (PI: GR) is the VLA part of a joint XMM-Newton and NRAO proposal. The X-ray data are presented in a separate paper (Rauw et al. 2024). We observed on three dates: one near periastron (2023-Mar-18) and the others approximately two weeks before (Mar-03) and after periastron (Apr-05). For all three observations, the VLA was in B configuration.

Data were obtained at 3.6 cm (X-band) and 6 cm (C-band) On each date and for each band, three sequences of observations were done, consisting of a ~4.6 min observation on HD 168112, preceded and followed by a ~1.5 min observation on the phase calibrator J1832-1035 (which is 3.7° from HD 168112). The total observing time on HD 168112 is therefore ~14 min, for each date and each band. At the end of the run, the flux calibrator 3C286 was observed (~1 min for C-band, ~1.5 min for X-band). The data were taken in 32 spectral windows of 128 MHz each, with the X-band covering 8.0–12.0 GHz and the C-band 4.0–8.0 GHz. Each spectral window in turn consists of 64 channels of 2 MHz each.

3 Analysis of spectroscopic data

Consistent with its O5 III spectral classification (as a single star, Holgado et al. 2018), the spectra of HD 168112 show prominent hydrogen, He I and He II lines, as well as some metal lines. The hydrogen and helium lines do not show P Cygni profiles, indicating a relatively weak stellar wind. The C III 5696, 6730, and N III 4634,40,41 lines are in emission.

While most of the spectra are single-lined, some of them show clearly double-lined features, indicating the binary nature of HD 168112. Some of the apparently single-lined spectra show line broadening, which we attribute to radial velocity differences which are too small to induce a double-lined profile.

In our analysis, we applied a spectral disentangling technique to the data (details are given in Sec. 3.2). To start up the disentangling procedure, we obtained a first estimate of the radial velocities and a first estimate of the orbital parameters (Sec. 3.1). We then used the disentangled spectra to update the radial velocities and redetermined the orbit. Iterating this a number of times led to the final orbital parameters (Sec. 3.3). From the disentangled spectra, we then found the rotational and macroturbulent velocities (Sec. 3.5). By comparison to TLUSTY models, we then determined the stellar parameters (Sec. 3.6).

3.1 Preliminary orbit

We determined preliminary radial velocities by fitting a selected set of spectral lines with one or multiple Gaussians. While there are a good number of spectral lines that show the line-doubling, we limited this analysis to those with a high enough signal-to-noise ratio (S/N). We therefore first determined the S/N around a large number of candidate lines. The S/N at 5410 Å is listed in Table A.1, and the variation of S/N across the spectra is shown in Fig. A.1. On the basis of this, we chose to continue with the following spectral lines: He I 4471, He I 5876, He II 4200, He II 4541, He II 5411, C IV 5801,12, and O III 5592.

In the region around each of these lines, we flagged any remaining cosmic rays. We then normalised the spectra on each of these regions by fitting a second-degree polynomial through interactively selected continuum points. Where the line is clearly split, we obtained two radial velocities (RVs); otherwise only one value was measured. For each spectrum, we then took an average of the individual RVs, weighting them with the inverse of the square of the individual error bars. In this process, we excluded any outlier RVs.

To determine the orbital period from these RV data, we applied the string length method (Dworetsky 1983). We rescaled the radial velocities of each component separately to the range 0..1. This ensured that the phase differences and the radial velocity differences have equal weight. Figure 1 shows that the best result for the period is ~516 days. The second lowest value is at ~255 days, which is close to half of the best period. We checked the spectra taken near the additional periastron passages of an assumed ~255 days period. Only data taken around JD 2 457 200 fulfil this criterion, and they clearly show no splitting of the spectral lines. We therefore took ~516 days as the (approximate) value for the period.

We then determined the preliminary orbit based on these RVs, using the Markov Chain Monte Carlo (MCMC) code emcee (Foreman-Mackey et al. 2013). This code explores many combinations of the orbital parameters T₀ (time of periastron passage), P (period), e (eccentricity), ω (longitude of periastron), γ (systemic velocity⁵), and K₁, K₂ (semi-amplitude radial velocity of primary and secondary). For each combination of these parameters the code calculates the model radial velocities (RVModel_1,i, RVModel_2,i) for the times of the observed spectra. The code then finds the solution with the smallest χ², where: $${\chi }^2={\displaystyle \sum _i}{\left(\frac{{\text{ RVObs}}_{1,i}-{\text{ RVModel}}_{1,i}}{{\text{ RVObsError}}_{1,i}}\right)}^2+{\left(\frac{{\text{ RVObs}}_{2,i}-{\text{ RVModel}}_{2,i}}{{\text{ RVObsError}}_{2,i}}\right)}^2,$$

where RVObs_1,i and RVObs_2,i are the observed radial velocities of the first and second component, with their associated error bars RVObsError_1,i and RVObsError_2,i. For those spectra where we could not split the lines, we used RVObs_2,i = RVObs_1,i.

In the MCMC procedure, we used no significant priors, except requiring an eccentricity between 0 and 1. We used 32 walkers and started them with a range in value that covers the above derived approximate period, as well as estimates of the semi-amplitude radial velocities. We ran it for 20 000 iterations and discarded the first 5000. Based on estimates of the autocorrelation times, we thinned the resulting samples by taking only one in every 200.

Fig. 1 String length as a function of the period. The inset shows a zoomin on the best value found.

3.2 Disentangling

Using the orbital parameters derived above, we used a disentangling technique on each of the spectral lines. For the disentangling, we used the code by Shenar et al. (2020, 2022). This code is based on the shift and add technique introduced by Marchenko et al. (1998) and González & Levato (2006). At each step in this iterative procedure, an approximate spectrum is available for each of the two binary components. The next step consists of subtracting the Doppler-shifted intrinsic spectrum of the first component from the observed spectra, thus giving spectra that have only contributions of the second component. Doppler-shifting these spectra and adding them up gives the next approximation of the intrinsic spectrum of the second component. The procedure is then applied again, exchanging the roles of the first and second component. After a sufficient number of iterations, the disentangled spectra of the two binary components are obtained.

We applied the disentangling part of the Shenar et al. (2020, 2022) code to each of the spectral lines we used in the radial velocity determination (Sec. 3.1). We applied 1000 iterations; checks with a higher number of iterations gave the same result. For each of these lines we thus obtained two separate lines, one for each component of the binary. We then explored for each spectrum and each spectral line a two-dimensional grid of velocities. For every point in the grid we shifted the two disentangled spectral lines by the two chosen velocities and added them together. This was then compared to the observed spectral line at that epoch, and the χ² difference was determined. From the minimum χ², the RVs of the two components at that epoch were then derived. In this procedure, we assumed that the light ratio between the two components is 1.0; the correct value is determined in Sec. 3.6.

The combined RVs for each epoch were then determined as the equal-weight average of the RVs of the individual spectral lines. Outliers were removed in this procedure. The error bar on the combined RVs is given by the standard deviation on the individual RV values that contributed to the combined one.

3.3 Final orbital parameters

To determine the final orbital parameters, we iterated between the disentangling and the orbital parameter determination outlined in the previous sections. At each step in the iteration, we used the radial velocities as input to the emcee code, which determined the orbital parameters. From these orbital parameters, the disentangling code found the spectral lines of the two binary components. These were then used to redetermine the radial velocities.

During these iterations we checked the convergence of the orbital parameters. We stopped after 20 iterations, as the parameters had then sufficiently converged. To obtain the final values, we ran a longer emcee run with 200 000 iterations, discarding the first 50 000. The corner plot presenting the results of this emcee run is shown in Fig. 2. The fit of our orbital solution to the observed radial velocities, for each spectral line, and for their combination, is shown in Fig. 3. The orbital parameters are listed in Table 1. As noted in Table 1, the mass ratio is very close to one, making it difficult to decide which of the two components is the primary. To avoid introducing confusion in the literature, we have assigned the subscript “A” to the component that Putkuri et al. (2023) designate as the primary, and the subscript “B” to the secondary.

The MCMC approach is very good at determining the error bar on the orbital parameters due to the error bars in the radial velocities. But we also wanted to get an estimate of any systematic errors. We therefore ran the whole iterative procedure again with some changes to the input. We experimented with the wavelength range covered by each of the spectral lines, as well as systematically dropping one of the spectral lines from the procedure. We extended the range of the MCMC error bars to also include the values covered by these variant analyses. The final error bars are listed in Table 1.

Fig. 2 Corner plot showing the final orbital parameters and their correlations. The time of the periastron passage (T0) and the period are strongly correlated.

3.4 Comparison with Putkuri et al. (2023)

A comparison with the Putkuri et al. (2023) orbital determination shows that they have a slightly longer period (513.52 days) than our value of 512.17 days. Converting their T₀ to a value we can compare with our Fig. 2, we find 54377.39. This combination of T₀ and period is well outside the corresponding plot in Fig. 2. The Putkuri et al. (2023) period has an unexpectedly small error bar: 0.01 day. Accumulated over the ~11 cycles covered by the time span of their observations, this would result in a phase difference of 2 × 10⁻⁴, which would be hard to detect. Our ~0.4 day error bar over ~9 cycles of the time span covered by our observations results in a phase difference of 7 × 10⁻³, which is just about detectable around periastron.

When we plot our radial velocities in an orbital phase diagram using the Putkuri et al. (2023) orbital solution (Fig. 4), we first note that our measured radial velocities are somewhat smaller than the Putkuri et al. (2023) orbital fit suggests. This difference is due to the different way we measure the radial velocities. Putkuri et al. (2023) fitted Gaussian profiles to the spectral lines, while we applied an iterative procedure using profiles from the disentangled spectra, as detailed in the above section.

The orbital fit to our data using the Putkuri et al. (2023) parameters has a peak at periastron that is less sharp. It is important to note that around periastron we have a higher number of observations than Putkuri et al. (2023, their Fig. 2), which allows us to better estimate T0 and the period. It also allows us to derive a better constraint on the radial velocity semi-amplitudes, which turn out to be somewhat smaller than their results (which are K_A = 70.4 ± 1.1 km s⁻¹, K_B = 77.5 ± 1.2 km s⁻¹). The smaller semi-amplitudes also explain the differences in the derived parameters of minimum semi-major axis, minimum masses, and mass ratio.

Our value for the eccentricity is slightly higher than the Putkuri et al. (2023) one (e = 0.743 ± 0.005), but the error bars overlap. The values for the longitude of periastron are equal within the error bars.

Fig. 3 Comparison of our orbital solution (solid line) to the observed radial velocities (symbols). Each panel shows the observed values for a specific stellar line, with the bottom panel showing the radial velocities combined over all spectral lines. Blue indicates the A component, red the B component.

Fig. 4 Comparison of our orbital solution to the Putkuri et al. (2023) one for our data around periastron. The phases on the left panel are calculated with our ephemeris, the ones on the right panel with the Putkuri et al. (2023) one. The symbols show our observed radial velocities (combined data), on both panels. The colour coding is as in Fig. 3. The solid lines indicate the orbital solutions for the respective components.

3.5 Projected rotational velocity

To determine the projected rotational velocity (υ sin i), we used the IACOB-BROAD code (Simón-Díaz & Herrero 2014). This code determines both the υ sin i and the macroturbulent broadening (υ_mac). One part of the code uses the Fourier transform of the spectral line, and from the first zero in the transform finds υ sin i (Simón-Díaz & Herrero 2007). Another part of the code explores a number of υ sin i and υ_mac combinations, calculates the resulting line profile, and uses a goodness-of-fit technique (χ² minimization) to determine the best-fit values.

Stark-broadened lines are not suitable for the application of this technique. We therefore limited it to the two metallic lines that have a sufficiently high signal-to-noise in our disentangled spectra: O III 5592 and C IV 5812 (the C IV 5801 line is not well disentangled due to a nearby Diffuse Interstellar Band).

For the A component, the O III 5592 line gave υ sin i_fou = 81 km s⁻¹ from the Fourier technique and υ sin i_gof = ${82}_{-32}^{+17}$ km s⁻¹ and ${48}_{-47}^{+47}$ km s⁻¹ from the goodness-of-fit. The C IV 5812 line gave, respectively, υ sin i_fou = 94 km s⁻¹, ${86}_{-23}^{+14}$ km s⁻¹, and ${72}_{-28}^{+31}$ km s⁻¹. From the range covered by these determinations, we assigned υ sin i = 75 ± 25 km s⁻¹ to the primary component. We could not determine a reliable value for the macroturbulent broadening, only an upper limit of υ_mac of ~ 100 km s⁻¹.

For the B component, we found υ sin i_fou = 141 km s⁻¹, $v\sin i_{\text{gof }}={141}_{-23}^{+13}\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}$, and $v_{\mathrm{m}\mathrm{a}\mathrm{c}}={64}_{-44}^{+52}\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}$ from O III 5592, and υ sin ifou = 129 km s⁻¹, $v\sin i_{\text{gof }}={130}_{-17}^{+12}\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}$, and $v_{\text{mac }}={57}_{-46}^{+39}\mathrm{k}\mathrm{m}{\mathrm{s}}^{-1}$ from C IV 5812. From the range covered by these determinations, we combined this to υ sin i = 134 ± 21 km s⁻¹ and an upper limit for υ_mac of ~110 km s⁻¹ for the B component.

3.6 Stellar parameters

To determine the stellar parameters, we first created a grid of TLUSTY models. We started from the models with solar composition in the public OSTAR2002 grid (Lanz & Hubeny 2003). We used the TLUSTY code⁶ to refine this grid, covering the T_eff range of 30 000–40 000 K with a step of 1000 K and the log g = 3.5–4.0 range with a step of 0.1. We also covered the range T_eff = 40 000–46 000 K with a step of 1000 K, log g = 4.0–4.7 with a step of 0.1. We used a helium abundance of 0.10 and a microturbulent velocity of 10 km s⁻¹.

To calculate the model spectrum, we used SYNSPEC (Hubeny & Lanz 2017). These synthetic spectra were then convolved with a rotational profile and a macroturbulent profile, using a code similar to IACOB-BROAD. For the rotational broadening of the A component, we used υ sin i = 75 km s⁻¹. As shown in the previous section, we could derive only an upper limit for the macroturbulent velocity, and we therefore took the Putkuri et al. (2023) value of υ_mac = 75 km s⁻¹. For the B component, we used our value of υ sin i = 134 km s⁻¹ and the Putkuri et al. (2023) value of υ_mac = 100 km s⁻¹. Finally we applied the instrumental broadening, using a Gaussian with the resolving power of the instrument.

Contrary to the procedure of Putkuri et al. (2023), we did not fit the disentangled spectra separately, but we compared combined synthetic spectra to the observed data. This has the advantage that it avoids the potential artefacts associated with disentangling. To reduce the amount of data we had to process, we selected spectra for a number of representative orbital phases (Table 2). We gave a higher preference to phases around periastron. As the spectral lines are clearly split at these phases, they allow for a better constraint on the stellar parameters of the two components.

We used the following spectral lines for fitting: Hγ, He I 4026, He I 4471, He I 5876, He II 4200, He II 4541, and He II 5411. We looped over all combinations of a theoretical spectrum for the A component and one for the B component. We applied the required radial velocity shift as given by our orbital model. The two spectra were added up using a light ratio of 0.55/0.45 for the primary over the secondary. The combined theoretical spectra were then compared to the observed spectrum, and the χ² difference between the two was determined. For the weights in the χ² calculation we used the inverse square of the S/N for the spectral line. The total χ² for a given combination is the sum of the χ² over all observed spectra. From the minimum χ², we found the best fit solution. As a final verification, we checked all the fits by eye. The best-fit results for a selected number of phases are shown in Fig. B.1.

While the above procedure uses the 0.55/0.45 light ratio proposed by Putkuri et al. (2023), we also tried alternative values. We then judged by eye which gave the better fit. We took into account that a different light ratio can result in a different best fit combination of stellar parameters. Based on this, we found that the range 0.525–0.575/0.475–0.425 gives acceptable results. This conclusion is mainly based on the He I 4471 and 5876 lines which are most sensitive to the change in light ratio.

The final results for the effective temperature and gravity of the best-fit TLUSTY models are listed in Table 3. We then used the extinction-corrected value for the V magnitude of the combined system derived by Putkuri et al. (2023): m_ν = 5.283 ± 0.022. Inverting the Gaia parallax of 0.4985 ± 0.0204 mas (HD 168112 = Gaia DR3 4153666106124493696), we have a distance of 2.01 ± 0.08 kpc, which then gave an absolute M_ν = 6.229 ± 0.092. Using the 0.525–0.575/0.475–0.425 light ratio we split this into M_V,A = −5.580 ± 0.104 and M_v,b = −5.362 ± 0.110. From the TLUSTY models, we could derive the bolometric correction, and thence the bolometric luminosity. We could then derive the spectroscopic values for radius and mass. All stellar parameters are listed in Table 3.

Our results present various challenges. A first one is the difference of our stellar parameters (derived from TLUSTY models) compared to those of Putkuri et al. (2023, who use the FASTWIND code). For the A component, our error bars on T_eff and log g overlap with those of Putkuri et al. (2023, T_eff = 41 700 ± 1200 K, log g = 3.77 ± 0.12). However, this is not the case for the B component, where our values are lower than theirs (T_eff = 40 500 ± 800 K, log g = 3.6–3.8 ± 0.10). A possible cause for this is the difference between the plane-parallel hydrostatic TLUSTY code and the FASTWIND code that treats in a unified way both the atmosphere and the (hydrodynamical) stellar wind in spherical symmetry. Another contributing factor is that our procedure was not applied to the disentangled spectra, but worked directly on the observed data. The disentangled spectra can contain artefacts, especially in the wings of the hydrogen lines, which are well known to influence the gravity determination.

For the absolute magnitude M_V, we achieved sharper error bars, because we relied on the light ratio determined above, rather than the measured magnitude difference of 0.17 ± 0.19 in H-band (Sana et al. 2014). For the spectroscopic bolometric luminosity, we find results that are compatible with the Putkuri et al. (2023, log L/L_⊙A = 5.64 ± 0.12, log L/L_⊙,B = 5.53 ± 0.08) ones within the error bars. Our radii are slightly higher.

Another challenge is that the masses we derived are not nearly equal, in contradiction to the mass ratio of ${0.996}_{-0.183}^{+0.105}$ we found from the orbit determination. Comparing our values with the theoretically expected ones from Martins et al. (2005), we find that our A component mass is too high (by at least 15%), while the B component is at least 30% too low. A comparison with Putkuri et al. (2023) shows that our B component mass is compatible with their result (M_B < 26 M_⊙), but our A component mass is higher (although we have marginally overlapping error bars with their result of M_A = 35.4 ± 8.2 M_⊙). Within the error bars of the Putkuri et al. (2023) mass for the A component and our value for the B component, we can indeed arrive at a ratio close to one. Much better limits on the masses of both components can be obtained by combining the spectroscopic data with astrometric information, as is discussed in Sec. 5.

4 Analysis of radio data

Calibration of the radio data was done by NRAO with the VLA CASA (Common Astronomy Software Applications) Calibration Pipeline⁷ 6.5.4–9. This pipeline flagged data affected by Radio Frequency Interference, and used the calibrator observations to calibrate the HD 168112 visibilities. It also applied the flux and band-pass calibration using the known flux of 3C286 = J1331+3030. We checked by eye the resulting flagging and calibration, and applied some further (light) flagging to the data, as needed.

We then used CASA (CASA Team 2022) version 6.5.6.22 to convert the HD 168112 calibrated data into images and to deconvolve them, using the tclean command. During this imaging process, we combined all spectral channels, leading to a single image for each of the two bands. We did this to optimize the S/N. The image fully covers the primary beam, and uses a pixel size of 0″. 10 × 0″.10 for 3.6 cm and 0″.15 × 0″.15 for 6 cm, which oversamples the synthesized beam well. Towards the edge of the 6 cm image, the resolved image of a radio galaxy is seen; for this reason we limited the uv-range of the visibility data we used, and we applied multi-scale cleaning. We used the Briggs (1995) weighting scheme with a 0.5 robust parameter.

We next used the imfit command to fit a single elliptical Gaussian to a small area around the centre of the cleaned image. The integrated flux from this fit and its error bar are listed in Table 4. In the error bar we also included a 5% absolute calibration error⁸ which has been added in quadrature to the measured flux error. To judge the systematic errors, we repeated the reduction a number of times, dropping data from one of the antennas, or one of the three time sequences of the observations. In all cases, the results fall within the error bar listed in Table 4.

From the fluxes at two wavelengths, we derived the spectral index α, which is defined by F_ν ∝ ν^α. This is also listed in Table 4. The interpretation of these data is discussed in the next section.

5 Comparison with other data

Radio and optical interferometry data of HD 168112 have been published in the literature. In this section, we compare our results to those data.

Blomme et al. (2005) analysed the available archive radio data on HD 168112 from the Very Large Array (VLA) and from the Australia Telescope Compact Array (ATCA). They found considerable variability in the fluxes, but the radio data did not allow for a precise determination of the orbital period. A value between one and 2 yr was proposed, with a most likely value of ~1.4 yr. The latter is in good agreement with our value of 512.17 days.

Now that the orbital parameters are known, we can plot again the radio fluxes of Blomme et al. (2005) in the orbital phase diagram (Fig. 5), and add the new observations. The radio fluxes show a clear correlation with the orbital phase. The fluxes are at minimum near periastron, and the new observations presented here (in red on Fig. 5) are in excellent agreement with the older data, which were taken ~20–40 yr before.

The spectral index (Fig. 5, bottom plot) allows us to distinguish between thermal and non-thermal radiation. Thermal radiation is due to free-free emission in the stellar winds of both stars and has a spectral index of α ≈ +0.6 (Wright & Barlow 1975; Panagia & Felli 1975). Where the two stellar winds collide, shocks form on either side of the collision region. Around these shocks, a fraction of electrons can attain relativistic speeds due to the Fermi mechanism (Bell 1978; Reitberger et al. 2014; Pittard et al. 2021). By spiralling around in the magnetic field, these electrons then emit synchrotron radiation, which we detect as non-thermal radiation (Eichler & Usov 1993). Theoretical models show that the non-thermal spectral index can be almost as low as −1.0 (Pittard et al. 2006).

As had already been noted by Blomme et al. (2005), the radio fluxes of HD 168112 clearly show a non-thermal spectral index, at least at orbital phases away from periastron. Near periastron, however, the spectral index is more consistent with thermal emission. We can first estimate the expected thermal flux from both stellar winds. Using the stellar parameters from Table 3, we applied the Vink et al. (2001) mass-loss recipe to determine the mass-loss rate. We found 2.1 × 10⁻⁶ M_⊙ yr⁻¹ for component A and 1.2 × 10⁻⁶ M_⊙ yr⁻¹ for component B. For the terminal velocity, we used the value of 3250 km s⁻¹ from Leitherer (1988). We derived the distance from the Gaia parallax. The Wright & Barlow (1975) formula for the theoretical radio flux then gives 0.029 mJy (component A) + 0.014 mJy (component B) at 3.6 cm and 0.022 + 0.010 mJy at 6 cm. This is at least a factor 3–4 lower than the observed fluxes near periastron. Clumping in the stellar wind increases the radio flux for a given mass-loss rate. The Vink et al. (2001) mass-loss rates do not include clumping, but recent work by, for example, Hawcroft et al. (2021) on spectroscopic data shows that the inclusion of clumping requires a decrease in the mass-loss rate (compared to Vink et al. 2001) in order to explain the strengths of the spectral lines. To first order, the effect of the clumping and the reduction of the mass-loss rate compensate one another. We therefore cannot attribute the observed radio fluxes near periastron to clumping in the winds. It is therefore probable that the wind-wind collision region is also contributing thermal radiation, as the theoretical models of Pittard (2010) predict.

In this highly eccentric binary, one expects the stellar wind collision to be strongest at periastron and therefore the highest synchrotron flux to be generated. The question therefore remains why we do not detect this synchrotron radiation. One possible reason is that the shocks collapse at periastron. However, this is contradicted by the X-ray flux variations that are symmetrical around the periastron passage (Rauw et al. 2024). Another possibility is that the synchrotron photons are being absorbed by the free-free absorption in the two stellar winds. From Table 1, we find that at periastron, the two components have a projected separation of 220 R_⊙. This should be compared to the typical size of the radio emitting region, for which we use the radius where the optical depth τ_λ = 1. Using the Wright & Barlow (1975) equations and the mass-loss rate and terminal velocity determined above, we find a 3.6 cm radius of 420 R_⊙ for component A and 300 R_⊙ for component B. At 6 cm, the values are a factor (6/3.6)^2/3 ≈ 1.41 higher. The introduction of clumping has little influence on these values because, again, to first order, the effect of the clumping and the reduction of the mass-loss rate compensate one another. At periastron the synchrotron emission region is therefore well inside the two stellar winds. It is only when the stars are further apart in their orbit (for example, at apastron, the projected separation is 1570 R_⊙) that a much higher fraction of the synchrotron photons can escape and be detected. This effect has also been seen in other binaries, such as WR 140 (Williams et al. 1990, 1994) and Cyg OB2 #9 (van Loo et al. 2008).

De Becker et al. (2024) observed HD 168112 at 1.6 GHz (18 cm) with the European Very Long Baseline Interferometry Network (EVN). They obtain a high spatial resolution image that resolves the colliding-wind region. The region is clearly asymmetric and its high flux (integrated 1.7 ± 0.4 mJy) shows it to be non-thermal in origin. The date of their observation (2019-Nov-05) corresponds to phase 0.60 according to our orbital solution. A direct comparison with the results presented here is not possible, as the orbital phase and wavelength of the observation are quite different. Blomme et al. (2005) listed a much higher 20 cm flux of 11 ± 7 mJy at phase 0.68 (project AC116, date 1984-Nov-27) and 11 ± 2 mJy at phase 0.71 (C978, 2001-Oct-11). However, these fluxes were derived from VLA data, which do not resolve the colliding-wind region. As De Becker et al. (2024) show, the radio fluxes are dominated by the colliding-wind region and the contribution of the two stellar winds is negligible. The higher VLA fluxes compared to the EVN ones therefore show that the synchrotron emitting region must be geometrically more extended than is shown in the EVN image. In the outer parts of this extended region, the intensity is lower than in the inner parts, and so the outer parts are not detected on the EVN image. However, their large geometric extent makes them contribute substantially to the total flux.

Sana et al. (2014) found a visual companion to HD 168112 using the VLTI PIONEER instrument. The companion was detected at a distance of 3.33 ± 0.17 mas at a position angle of 303°27 ± 4° 12. The magnitude difference in the H band is 0.17 ± 0.19. The date of the observation is 2012.4446, which corresponds to JD 2 456 091.224. Using our orbital solution we found this corresponds to phase 0.32. At that phase the radius divided by the semi-major axis is r_A/a_A ≈ r_B/a_B = 1.6.

We can convert the angular separation of 3.33 mas to a projected linear distance of ~1436 R_⊙ using the Gaia parallax. To compare this with our spectroscopic orbit, we need to use the following astrometric conversion (Gallenne et al. 2023): $$\begin{array}{rl}\mathrm{\Delta }\alpha & =r[\sin \mathrm{\Omega }\cos (\omega +v)+\cos i\cos \mathrm{\Omega }\sin (\omega +v)],\\ \mathrm{\Delta }\delta & =r[\cos \mathrm{\Omega }\cos (\omega +v)-\cos i\sin \mathrm{\Omega }\sin (\omega +v)],\\ r& =\frac{a(1-e^2)}{1+e\cos v},\end{array}$$(1)

where Δα and Δδ are the offset in right ascension and declination, r is the (non-projected) separation, Ω is the longitude of the ascending node, ν is the true anomaly, and the other parameters have been defined in Sec. 3.1. The projected separation is then given by: $$\begin{array}{rl}\rho & =\sqrt{(\mathrm{\Delta }\alpha {)}^2+(\mathrm{\Delta }\delta {)}^2}\\ & =\frac{a(1-e^2)}{1+e\cos v}\sqrt{{\cos }^2(\omega +v)+{\cos }^{2}i{\sin }^2(\omega +v)}.\end{array}$$

From the spectroscopic value of a sin i and the astrometrically measured ρ, we then find the inclination to be ~63°.

Using this value for the inclination and the minimum masses from Table 1, we then find the masses to be ~26 M_⊙. This is compatible with the spectroscopic mass derived for the B component, but too small to be consistent with the A component (Table 3).

The position angle θ (measured from the north to the east) can be derived from: $$\theta =\frac{\pi }{2}-\arctan \left(\frac{\mathrm{\Delta }\delta }{\mathrm{\Delta }\alpha }\right),$$

where care must be taken to assign θ to the correct quadrant. From Eq. (1), we can also derive Ω, although there is an ambiguity as we do not know if the measured position angle is from component A to B or vice versa. Assuming the 303°27 position angle is from A to B, we find Ω ≈ 318°, otherwise Ω ≈ 138°.

The De Becker et al. (2024) resolved radio observation also provides an estimate of the position angle. On their Fig. 1, the colliding wind region is asymmetric with the west side wrapping around the star with the weakest wind (which we identify with the B component). From this we can estimate an A to B position angle⁹ of θ ≈ 298°. From this we find Ω ≈ 301°, which is compatible with one of the results from the optical interferometry.

Fig. 5 Radio fluxes at 3.6 cm (top panel) and 6 cm (middle panel), plotted as a function of the orbital phase in the 512.17 days period, and the corresponding spectral index (bottom panel). The black symbols show the data from Blomme et al. (2005) while the red ones give the data presented in this paper. The insets zoom in on the phase range around periastron.

6 Conclusions

We analysed 113 spectroscopic observations of HD 168112 taken during a 13 yr time span with the HERMES and TIGRE spectrographs. In this way, we obtained an orbital solution that is independent of the one determined by Putkuri et al. (2023). We find a period of ${512.17}_{-0.11}^{+0.41}$ days, an eccentricity of ${0.7533}_{-0.0124}^{+0.0053}$, and a mass ratio close to one. A comparison with the Putkuri et al. (2023) results shows that we have a slightly shorter period. To determine the radial velocities for each observation, we used an iterative approach with the disentangled spectra providing the shape of each spectral line. This approach differs from the Putkuri et al. (2023) one, who used Gaussian profiles. As a consequence, we find slightly smaller semi-amplitudes for both components, and therefore also smaller lower limits on the masses.

We also derived the stellar parameters for the two components. Contrary to Putkuri et al. (2023) we did not use the disentangled spectra for this, but instead tried various combinations of two model spectra and compared the combined spectrum (with the appropriate radial velocity shifts) to the observed spectra at various orbital phases. Within the (large) error bars we find similar T_eff and log g for the A component, but our B-component values are lower than the Putkuri et al. (2023) ones. The difference could be due to the different procedure we used (avoiding in our case the known disentangling problems), or to our use of plane-parallel TLUSTY models instead of the Putkuri et al. (2023) FASTWIND ones.

We also analysed radio observations that were taken at three epochs close to the periastron passage. We expect the synchrotron emission around the colliding wind region to result in a negative spectral index. But the spectral index of all three observations shows only thermal emission, suggesting that most of the synchrotron photons are absorbed in the two stellar winds. However, the total flux cannot be attributed to the stellar winds only. The colliding-wind region must therefore also be contributing thermal radio emission. Clearly, this effect can also occur in other colliding-wind binaries.

By combining our data with the optical interferometry result of Sana et al. (2014), we can constrain the inclination angle (i ~ 63°), and in this way estimate the masses of both components (~26 M_⊙). For the A component, this is considerably lower than the spectroscopic estimate, while the value is consistent with the B component. As the error bars on the orbital parameters are relatively small, we have a higher confidence in the masses derived in this way, rather than the spectroscopic masses. Finally, we showed that our orbital parameters are consistent with the resolved radio image obtained by De Becker et al. (2024).

Acknowledgements

We are grateful to the observers at the Mercator Telescope for collecting the data of our HERMES service mode programme. The Liège team acknowledges support from the Fonds National de la Recherche Scientifique (Belgium) and the Belgian Federal Science Policy Office (BELSPO) in the framework of the PRODEX Programme (contract HERMeS). This work has used the following software products: emcee (Foreman-Mackey et al. 2013), corner.py (Foreman-Mackey 2016), Matplotlib (Hunter 2007, https://matplotlib.org), NumPy (Harris et al. 2020, https://numpy.org), IDL (https://www.nv5geospatialsoftware.com/Products/IDL), Shenar et al. (2020, 2022) disentangling code, IACOB-BROAD code (Simón-Díaz & Herrero 2014), CASA (CASA Team 2022). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France, and NASA’s Astrophysics Data System. We thank the National Radio Astronomy Observatory (NRAO) for carrying out the Karl G. Jansky Very Large Array (VLA) observations and applying the CASA pipeline reduction. We thank the referee for comments that helped clarify the paper.

Appendix A Observing log

Fig. A.1 Relative signal-to-noise ratio (S/N) as a function of the wavelength for the HERMES (top panel) and the TIGRE (bottom panel) spectra. The S/N is relative to the S/N at 5410 Å, with the 5410 Å value being listed in Table A.1. At a given wavelength, the multiple crosses show the value of the individual spectra for the spectral line. The red line gives the median of these values, and the green lines the range that contains 68.3% of the data. The break in the TIGRE data around 5700 Å is due to the change from the blue arm to the red arm (Schmitt et al. 2014).

Appendix B Spectral line fitting

Fig. B.1 Best fit to the spectral lines for a selection of orbital phases, using the stellar parameters from Table 3. Each column shows an orbital phase (listed at the top), each row shows a spectral line (listed inside each plot). The observed spectrum is given by the black curve, the theoretical fits of each of the two components by the blue curves, and the combined fit by the red curve.

References

All Tables

All Figures

	Fig. 1 String length as a function of the period. The inset shows a zoomin on the best value found.
In the text

	Fig. 2 Corner plot showing the final orbital parameters and their correlations. The time of the periastron passage (T0) and the period are strongly correlated.
In the text

	Fig. 3 Comparison of our orbital solution (solid line) to the observed radial velocities (symbols). Each panel shows the observed values for a specific stellar line, with the bottom panel showing the radial velocities combined over all spectral lines. Blue indicates the A component, red the B component.
In the text

Fig. 4 Comparison of our orbital solution to the Putkuri et al. (2023) one for our data around periastron. The phases on the left panel are calculated with our ephemeris, the ones on the right panel with the Putkuri et al. (2023) one. The symbols show our observed radial velocities (combined data), on both panels. The colour coding is as in Fig. 3. The solid lines indicate the orbital solutions for the respective components.

In the text

	Fig. 5 Radio fluxes at 3.6 cm (top panel) and 6 cm (middle panel), plotted as a function of the orbital phase in the 512.17 days period, and the corresponding spectral index (bottom panel). The black symbols show the data from Blomme et al. (2005) while the red ones give the data presented in this paper. The insets zoom in on the phase range around periastron.
In the text

Fig. A.1 Relative signal-to-noise ratio (S/N) as a function of the wavelength for the HERMES (top panel) and the TIGRE (bottom panel) spectra. The S/N is relative to the S/N at 5410 Å, with the 5410 Å value being listed in Table A.1. At a given wavelength, the multiple crosses show the value of the individual spectra for the spectral line. The red line gives the median of these values, and the green lines the range that contains 68.3% of the data. The break in the TIGRE data around 5700 Å is due to the change from the blue arm to the red arm (Schmitt et al. 2014).

In the text

	Fig. B.1 Best fit to the spectral lines for a selection of orbital phases, using the stellar parameters from Table 3. Each column shows an orbital phase (listed at the top), each row shows a spectral line (listed inside each plot). The observed spectrum is given by the black curve, the theoretical fits of each of the two components by the blue curves, and the combined fit by the red curve.
In the text