© ESO, 2017

1. Introduction

Massive stars constitute a small part of the stellar population but have a considerable influence on their environment. About 70% of them (Sana et al. 2012), and even 90% if we consider orbital periods longer than 3500 days (Sana et al. 2014), are moreover expected to be binaries or multiple systems. This multiplicity can affect the way that they evolve through tidally-induced-rotational mixing, exchange of matter and angular momentum through Roche lobe overflow, or even the coalescence of both components. In the case of young systems with large separations, the different components evolve independently, following single-star evolution (Pavlovski & Hensberge 2005) whilst, when the components interact with each other, their evolution can completely differ from that of single stars, affecting their size, mass, surface chemistry and/or even their rotational rate.

To better understand how massive stars evolve and improve evolutionary models, comparisons between the latter and spectroscopic observations are essential. Nevertheless, for single stars, the modeling by atmosphere codes does not provide all the information about their fundamental properties, especially on the real masses of the stars. Herrero et al. (1992) raised the attention to the mass discrepancy problem. This term refers to a disagreement between the mass determined by comparing the position of the star in the Hertzsprung-Russell (HR) diagram according to the predictions of stellar evolution calculations and the mass obtained from the spectroscopically derived values of the surface gravity and of the radius of the star. In this context, detached eclipsing binary systems are of crucial importance to provide direct measurements of the physical properties of stars such as their radius or their dynamical mass, allowing us to test the predictive abilities of evolutionary models for single stars.

HD 166734 (=V411 Ser) is particularly promising for such an analysis. Unfortunately, this system has neither a measured inclination nor a cluster membership. It was located by Conti et al. (1980) at a distance of 2.3 kpc from Earth. It is composed of two supergiant O-type stars in a wide eccentric orbit. This configuration, rather rare for such stars, suggests few interactions between each other, making them perfect objects for testing the theoretical schemes of stellar evolution. HD 166734 has not been investigated since the seventies when Conti et al. (1980) provided for the first time the parameters of its orbit. These authors described this system as an O7.5If primary and an O9I secondary moving on an eccentric orbit (e = 0.46) with a period of 34.54 days. These two objects had similar luminosities but the secondary presents a slightly larger minimum mass than the primary (31 M_⊙ vs. 29 M_⊙, respectively). On the basis of these minimum masses Conti et al. (1980) predicted eclipses in the light curve of HD 166734. Even though Otero & Wils (2005) mentioned HD 166734 as photometrically variable, the light curve has never been subject to a dedicated analysis.

The present paper aims at showing the advances that have been achieved in the analysis of this interesting system and is organized as follows. The spectroscopic and photometric data are described in Sect. 2, together with their reduction procedure. The radial velocity (RV) determination and the orbital solution of HD 166734 are presented in Sect. 3. Section 4 reports the analysis of the VRI light curves. In Sect. 5, we disentangle the spectra of each component and we model them with an atmosphere code to constrain the individual parameters and the surface abundances of each star. Section 6 is devoted to the winds of both components and to their interactions. Section 7 provides the evolutionary statuses of both components and discusses the distance to HD 166734. Finally our conclusions are given in Sect. 8.

2. Observations and data reduction

2.1. Spectroscopy

HD 166734 was observed through a multi-epoch campaign with two different instruments. We first acquired high-resolution optical spectra with FEROS (Fiber-fed Extended Range Optical Spectrograph) successively mounted on the ESO-1.52 m (for the observations taken before 2002) and on the ESO/MPG-2.2 m telescopes at La Silla (Chile). This instrument has a resolving power of about 48 000 and the detector was a 2k × 4k EEV CCD with a pixel size of 15 μm × 15 μm. It provides échelle spectra composed of 39 orders, allowing us to cover the [3800−9200] Å wavelength domain. The data reduction was performed with an improved version of the FEROS pipeline as described by Sana et al. (2006). The data normalization was performed by fitting polynomials of degree 4 or 5 to carefully chosen continuum windows.

We completed this dataset by acquiring medium-resolution échelle spectra with the ESPRESSO spectrograph. This instrument is mounted on the 2.12 m telescope at Observatorio Astronómico Nacional of San Pedro Mártir (SPM) in Mexico. The spectra are slitted over 27 orders to cover the [3800–6950] Å region. The spectral resolving power is of about 18 000. To increase the signal-to-noise ratio of our data, we took several consecutive observations and we combined them. The data were reduced by using the échelle package available within the ESO-MIDAS software.

The journal of the observations is provided in Table A.1.

2.2. Photometry

We obtained the light curve of HD 166734 in the VRI filters with the TAROT (Télescope à Action Rapide pour l’Observation des phénomènes Transitoires) robotic observatory located at La Silla. It is a fully autonomous 25 cm aperture telescope (F/D = 3.4) covering a 1.86° × 1.86° field of view thanks to an Andor CCD camera. Its spatial sampling is 3.3″/pixel. The reduction method of the TAROT observations is described in Damerdji et al. (2007). The object was observed between March 2011 and September 2016 in the V filter, between February 2010 and September 2016 in the R filter and between February 2011 and September 2016 in the I filter. These three light curves are composed of 5021, 5178, and 4928 points, respectively.

3. Radial velocity measurements and orbital solution

Figure 1 shows FEROS spectra as a function of the orbital phase (computed from orbital parameters determined below and listed in Table 1). It illustrates the variability of the line profiles for two particular regions featuring the He i 4471 absorption line as well as the C iii 5696 emission line.

Fig. 1 FEROS spectra in the region of He i 4471 (left) and C iii 5696 (right) shifted vertically as a function of phase (times 3 for clarity) from 0.0 at the bottom to 1.0 at the top.

To construct the orbital solution of HD 166734, we first measure the RVs of both components. For this purpose, we fit their spectral lines with Gaussian profiles. The computations of the RVs are done by adopting the rest wavelengths from Conti et al. (1977) for the lines with a rest wavelength shorter than 5000 Å, and those of Underhill (1995) for the others. Then we refine those values (notably for the spectra where the lines of both components are blended) by using a new method designed by one of the authors (YD). The latter is based on a minimum distance method considering a comparison between the observed spectrum and a composite (two-components) theoretical template. We perform this measurement for several lines individually, no matter their ionization stage. The final set of RVs (provided in Table A.1) is obtained by computing the mean values of the radial velocities measured on all the various spectral lines.

The determination of the orbital period is performed from the Heck-Manfroid-Mersch technique (hereafter HMM, Heck et al. 1985) that was revised by Gosset et al. (2001, see their appendix). This technique is basically a Fourier method generalized to address the unevenly sampled time series. The HMM power spectrum is determined independently of the mean of the data and thus does not suffer from the biases of the Scargle (1982) method. As input of this technique, we first use the RVs measured on the N iv 4058 line that is only visible in the primary’s spectrum. The periodogram exhibits two outstanding peaks at ν = 0.028453 d^-1 and at ν = 0.028890 d^-1 with a strong neighbor at ν = 0.028981 d^-1. However, given their large width, the two last-quoted frequencies overlap. We then use HMM on our final set of RVs (provided in Table A.1), on the basis of the difference between the secondary RVs and the primary RVs (RV_S–RV_P), and detect a frequency at ν = 0.028890 d^-1, giving the period of 34.614053 days and a neighbor at ν = 0.028957 d^-1, corresponding to a period of 34.534560 days. We determine an error of ϵ(f) = 1.1 × 10^-6 d^-1 (corresponding to an error on the period of 1.3 × 10^-3 days) calculated on our dataset from the expression given by Lucy & Sweeney (1971; see also Montgomery & O’Donoghue 1999; or Mahy et al. 2011, for further explanations). We compare these values to the result obtained by Conti et al. (1980). We stress that these authors determined an orbital frequency of ν = 0.028952 d^-1 from the difference between the RVs of the two stars, but that the computed periodogram also shows a neighbor at ν ~ 0.02896 d^-1 that overlaps with the other peak. To decrease the widths of the peaks, we thus need a dataset covering a longer timescale. We merge the two datasets and we determine in the periodogram once again two outstanding peaks at ν = 0.028879 d^-1, corresponding to a period of 34.627837 days and at ν = 0.028954 d^-1, corresponding to a period of 34.537542 days. By coupling the dataset of Conti et al. (1980) with our dataset, the error on the frequency becomes ϵ(f) = 2.4 10^-7 d^-1 which corresponds to an error on the period of 2.8 10^-4 days. The periodogram is exhibited in Fig. B.1. We then use the Liège Orbital Solution Package (LOSP¹) to determine the SB2 orbital solution of the system. The refined value of the orbital period, provided by this program (in order to improve the χ² of the fit) is 34.537723 ± 0.001330 days and is adopted in the following. We thus confirm the period computed by Conti et al. (1980) on the basis of the difference between their secondary RVs and their primary RVs and we emphasize that our refined value is similar to the period associated to the alias found in the periodogram computed from the two coupled datasets. The parameters are listed in Table 1 whilst the RV curves are shown in Fig. 2. From these orbital parameters, the inferior conjunction is expected at Φ ~ 0.018 and the superior conjunction is at Φ ~ 0.766.

Fig. 2 Radial velocity curves of HD 166734 computed with the RVs from the present paper (provided electronically). The primary is represented by filled symbols whilst the open ones indicate the RVs of the secondary. Circles represent the RVs measured on the FEROS spectra and diamonds indicate the RVs measured on the ESPRESSO data. Φ = 0.0 corresponds to the periastron passage.

Our orbital solution presents an orbit more eccentric (0.62 vs. 0.46, respectively) than that derived by Conti et al. (1980). Moreover, our monitoring (more intensive near the periastron passage) allows us to better constrain the semi-amplitudes, and thus the minimal mass of each component. In this new orbital solution, the primary is the most massive component of the system, which was not the case in the analysis of Conti et al. (1980). This result seems thus to suggest that the two components of HD 166734 show a “more classical” evolution than was first suspected (Conti et al. 1980).

4. Photometry

Before studying the light curve to constrain the parameters of both objects, we submit the three photometric datasets of HD 166734, obtained in the VRI filters, to a detailed Fourier analysis by HMM in order to search for the periods. Surprisingly, no outstanding peak is detected around 0.0289 d^-1 as would have been expected from the orbital frequency. Instead, we detect in each of the three datasets a series of peaks spaced with the same frequency of 0.0289 d^-1 (similar to the orbital frequency), starting at ν = 0.087 d^-1 until ν = 0.492 d^-1 (Fig. 3). The dominant member varies from one filter to another, going from ν = 0.174 d^-1 in the V and R filters to ν = 0.145 d^-1 in the I filter (see Figs. B.2, B.3). This series of frequencies corresponds in fact to harmonics of the orbital frequency, even though the fundamental and the first-order harmonics are not detected. This large number of harmonics is explained by the small size of the eclipse with respect to the orbital cycle. Moreover, these harmonics have aliases between 0.5 and 1.0 d^-1 and at higher frequencies. We also apply HMM on the publicly available ASAS data and still find most of these harmonics there.

Fig. 3 Top: power spectrum of the V light curve of HD 166734. The diamond indicates the location of the orbital frequency while the dots represent its harmonics. Bottom: spectral window computed from the V light curve.

Because of the inclination of the orbit and the orientation of its apses with respect to our line of sight, the light curve of HD 166734 exhibits only one eclipse during its orbital cycle, immediately after the periastron passage (Φ = 0.024) with a total duration in orbital phase of ΔΦ ~ 0.04. The phase of this eclipse is in agreement with the estimation of Otero & Wils (2005) but differs from the phase of the conjunction (see Sect. 3). This difference can come from the errors on the orbital parameters. We emphasize that if we fix the position of the eclipse to the conjunction it does not change the values of the orbital parameters within the uncertainties. This eclipse shows a difference of magnitude of about 0.22 mag and is preceeded by an increase of the total brightness of the system by about 0.03 mag at the periastron passage. We display in Fig. 4 the light curve of HD 166734 in the three different filters sampled over the orbital period. We stress that these light curves contain outliers (with an unknown origin) that do not change their general shape.

To limit the possible effects of the degeneracies that can be produced because of the presence of only one eclipse we use the PHOEBE (PHysics Of Eclipsing BinariEs, v0.31a, Prša & Zwitter 2005) software while keeping the parameters of the orbit fixed at their values obtained through the spectroscopic analysis (Sect. 3, Table 1). PHOEBE allows one to model the light curve and the RV curves of the object at the same time. This software is based on Wilson & Devinney’s code (Wilson & Devinney 1971) and uses Nelder & Mead’s Simplex fitting method to adjust all the input parameters and to find the best fit to the light curve. Furthermore, it allows one to model multiple bandpasses at once. We include the reflection effects between the two stars in the computation of the light curve and we also fix the effective temperatures to 32 000 K and 30 500 K, for the primary and the secondary components, respectively, as determined in Sect. 5. The photometric parameters derived from the modeling of the light curve are given in Table 2. The error bars are computed by exploring the parameter space. For this purpose, we fix one parameter and allow the others to vary in order to reach the minimum of the χ² corresponding to a 68.3% confidence level (1σ). Figure 4 displays the best-fit of the light curve obtained by PHOEBE in the three different bandpasses.

Fig. 4 Light curves of HD 166734 obtained with TAROT in the V (top panel), R (second panel) and I (third panel) filters. The red lines represent the best fit computed with PHOEBE. The magnitude outside the eclipse has been set to 0 in the three filters. The representation shown at the bottom is a zoom-in on the eclipse observed in the V light curve (plotted at the top of this figure).

The system is seen under an inclination of 63.0° ± 2.7°. The eclipse occurs when the primary is in front of the secondary (secondary eclipse) along the line of sight. The two components do not fill their Roche lobe over the entire orbit. They fill about 90%, within an error of 4%, at the periastron displaying a small deformation because of tidal interactions. With this inclination the masses obtained from the orbital solution are estimated to $\mathrm{39.}{\mathrm{5}}_{-4.4}^{\mathrm{+}\mathrm{5.4}}$ M_⊙ for the primary and $\mathrm{33.}{\mathrm{5}}_{-3.7}^{\mathrm{+}\mathrm{4.6}}$ M_⊙ for the secondary. Their radii are computed to 27.5 ± 2.3 R_⊙ and 26.8 ± 2.4 R_⊙, respectively.

The secondary component thus appears slightly more evolved (cooler and about the same size) than the primary. The modeling of the light curve also reports a primary component brighter than the secondary. From their radii and their effective temperatures, we determine stellar luminosities of $\log \frac{\mathit{L}}{{\mathit{L}}_{\mathrm{\odot }}}\mathrm{=}\mathrm{5.840}\mathrm{\pm }\mathrm{0.092}$ for the primary and of $\log \frac{\mathit{L}}{{\mathit{L}}_{\mathrm{\odot }}}\mathrm{=}\mathrm{5.732}\mathrm{\pm }\mathrm{0.104}$ for the secondary, infering a brightness ratio between the two stars of l₁/l₂ = 1.28 ± 0.14.

5. Individual parameters

5.1. Disentangling

We use the orbital parameters derived in Sect. 3 as input of the disentangling program. We apply the Fourier approach of Hadrava (1995) to separate the spectral contributions of both components. This method uses Nelder & Mead’s Downhill Simplex on the multidimensional parameter space to reach the best χ² fit between the recombined component spectra and the observed data. We emphasize that we only apply this technique on the FEROS spectra because they have a better signal-to-noise ratio and a higher resolution.

HD 166734 is a detached eclipsing binary. During the eclipse, the secondary component is however not totally eclipsed. There exists therefore an ambiguity in the determination of the continuum level. The method of Fourier presented by Hadrava (1995) induces oscillations in the continuum of the disentangled spectrum of each component when one of them is not completely eclipsed, as is the case for HD 166734. A renormalization of the disentangled spectra must thus be done with the light ratio derived from the light curve in Sect. 4 and provided in Table 2. A discussion of the spectral disentangling can be found in Pavlovski & Hensberge (2010) and in Pavlovski & Southworth (2012).

Fig. 5 Best CMFGEN fit (red) of the disentangled optical spectrum (black) of the primary component of HD 166734.

Fig. 6 Same as Fig. 5 but for the secondary component of HD 166734.

To better constrain the continuum level of both disentangled spectra, we decided, however, to compare them with those obtained from the method of Marchenko et al. (1998), revised by González & Levato (2006). As already stressed by Mahy et al. (2015), the technique of González & Levato (2006) gives reliable results outside broad lines such as the Balmer lines. In the wings of those lines, the disentangled spectra exhibit intrinsic excesses or deficiencies of the flux, making them difficult to normalize, and infering larger uncertainties on the shape of their wings. The two techniques have thus their respective advantages, and their comparison allows us to remove the artefacts linked to the disentangling method.

As mentioned above, the resulting spectra are normalized to remove any oscillation of the continuum linked to the Fourier disentangling technique and then corrected for the brightness ratio of l₁/l₂ = 1.28. The final disentangled primary and secondary spectra are displayed in Figs. 5 and 6, respectively.

5.2. Spectral-type classification

From the final disentangled spectra, we determine the equivalent widths of several lines to apply Conti’s criteria (Conti & Alschuler 1971; Conti 1973; Mathys 1988, 1989) and derive the spectral classifications of both components. We measure ratios of log (He i 4471/He ii 4542) = 0.06 and 0.31, which gives spectral classifications of O7.5 for the primary and of O9 for the secondary. We can then determine the luminosity class of each component using log (Si iv 4089/He i 4143). We find that this value is equal to 0.93 and to 0.65 for the primary and the secondary components, respectively, indicating that both stars are supergiant. We observe the lines of He ii 4686 and N iii 4634−41 in strong emission in the spectrum of the primary. According to Walborn (1971), we add the “f” suffix to the spectral type of the primary. For the secondary, we observe in its final disentangled spectrum that the absorption of the He ii 4686 line is compensated by the emission, whilst the N iii 4634–41 lines are moderately in emission, infering a “(f)” suffix to its spectral classification.

We also compare the component spectra of both stars of the system to the spectral atlas of Walborn & Fitzpatrick (1990), confirming their classification.

In summary, we attribute to the primary an O7.5 If spectral type whilst the secondary is classified as an O9 I(f) star, with an uncertainty of one subtype for each component. These classifications are in agreement with those obtained by Conti et al. (1980), within the uncertainties.

5.3. Modeling

To determine the physical parameters of both components, we use the CMFGEN atmosphere code (Hillier & Miller 1998) on the disentangled spectra corrected for the brightness ratio. This code provides non-LTE atmosphere models including winds and line-blanketing. CMFGEN requires input from an estimate of the hydro-dynamical structure connected to a β velocity law of the form v = v_∞(1−R/r)^β, where v_∞ is the terminal velocity of the wind. Our final models include the following chemical elements: H i, He i-ii, C ii-iv, N ii-v, O ii-v, Al iii, Ar iii-iv, Mg ii, Ne ii-iii, S iii-iv, Si ii-iv, Fe ii-vi, and Ni ii-v with the solar composition of Grevesse et al. (2007) unless otherwise stated. CMFGEN allows a super-level approach to reduce the memory requirements. Therefore, we include about 1771 super levels for a total of about 8075 levels. The formal solution of the radiative transfer equation allows us to compute an emergent spectrum. To do so, we use a depth-independent microturbulent velocity varying linearly from 20 km s^-1 at the photosphere to 0.1 × v_∞ at the top of the atmosphere. The choice of 20 km s^-1 for the microturbulence velocity at the photosphere is coherent with other studies made on supergiant massive star populations (see Bouret et al. 2012; or Crowther & Evans 2009). Once the emergent spectrum is calculated we convolve it twice, once with a rotational profile to include the projected rotational velocity (vsini) and once with a Gaussian profile to include macroturbulence (v_mac) in the line profiles. vsini is determined from the Fourier transform method given by Simón-Díaz & Herrero (2007) whilst v_mac is fixed by eye from the He i 4713 and O iii 5592 lines. We compute values of 95 ± 10 km s^-1 and 65 ± 10 km s^-1 for the vsini and the v_mac of the primary whilst values of 98 ± 10 km s^-1 and 65 ± 10 km s^-1 are obtained for the secondary, respectively.

The effective temperature of each component is estimated from the ratio between He i and He ii, mainly from the He i 4471 and He ii 4542 lines, even though the other helium lines are also considered. The typical uncertainty on our determination is 1000 K. The surface gravities of both components are determined from the wings of the Balmer lines Hγ and Hδ. The typical errors on the surface gravities of both components are 0.1 dex. Hα and Hβ are not taken into account because they are more affected by the stellar winds, or maybe even by the interaction between the winds of the two components. Furthermore, the shapes of these two lines can vary with respect to the orbital phase. In the final disentangled spectra, we obtain mean profiles computed over the entire orbital period that are not necessarily meaningful. The stellar luminosity is determined from the radii derived through the photometry (see Sect. 4) and from the effective temperature determined for each component through the CMFGEN best-fit (in an iterative way). The mass-loss rates and the clumping filling factors are adjusted to reproduce the Hα and Hβ lines as well as the He ii 4686 line whilst the shape of their emission is sensitive to the β exponent of the wind velocity law. The terminal velocity has been computed from the Newtonian escape velocity ${\mathit{v}}_{\mathrm{esc}}\mathrm{=}\sqrt{\mathrm{2}\mathit{G}\hspace{0.17em}\mathit{M}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{\Gamma }\mathrm{)}\mathit{/}\mathit{R}}$ (where G is the gravitational constant, Γ the Eddington factor computed from Vink et al. 2001; by assuming an electron scattering cross-section per unit mass of 0.32 cm² g^-1, Lamers & Leitherer 1993, M the mass of the star and R its radius) by assuming v_∞/v_esc = 2.6.

The best-fit model has been obtained by minimizing the calculated χ² by varying different parameters in the parameter space. To this end, we generate a non-uniform grid composed of several dozen models for each star. Once all the fundamental parameters are constrained, we run models with different surface abundances (for He, C, N, and O). Specific lines of each element are selected from the disentangled spectra and we quantitatively compare these lines to the synthetic spectra by means of a χ² analysis from which we derive the surface abundance and their uncertainties. This method is similar to that proposed by Martins et al. (2015), except that, rather than renormalizing the χ² function so that the minimum has a value of 1.0, we subtract an appropriate quantity to set the minimum to 1.0. The 1σ uncertainty is approximately set to the abundances for which ${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\mathrm{1}\mathrm{+}\mathrm{1}\mathit{/}{\mathit{\chi }}_{\min }^{\mathrm{2}}$. The lines used to constrain the abundances are C iii 4068–4070, and C iii 4153–56–63 for carbon, N iv 4058 (only for the primary), and N iii 4511–4515–4524 for nitrogen and O iii 5592 for oxygen. We exclude the C iii 4647–50–51, C iii 5696, and C iv 5801–5812 lines since their formation process depends on fine details of atomic physics and on the modeling (Martins & Hillier 2012). Moreover, special care is taken for the helium abundance. Some of the helium lines are also affected by the microturbulent velocity taken to compute the emerging spectrum (McErlean et al. 1998). A compromise is thus made to ensure the best CMFGEN fit. The individual parameters of both components are given in Table 3 and we display, in Figs. 5 and 6, the CMFGEN best-fit model for the primary and the secondary, respectively.

The individual parameters agree with two supergiant objects according to the tables of Martins et al. (2005), and Bouret et al. (2012). The best-fit is achieved with T_eff = 32 000 K and 30 500 K and log g = 3.15 and 3.10 for the primary and the secondary, respectively. We also derive an overabundance of helium and nitrogen as well as a depletion of carbon in both objects, as expected from their luminosity class. Furthermore, we emphasize that the oxygen is depleted for the primary whilst it is solar for the secondary.

The wind parameters are reminiscent of those of Galactic supergiant O-type stars, according to Bouret et al. (2012) even though stars as cool as the two components of HD 166734 were not analyzed by these authors. These parameters are thus in agreement with those of single massive supergiant objects.

6. Wind-wind interactions

The spectral disentangling has yielded the mean spectrum of each component. From these spectra, the emission observed in the He ii 4686 and the Hα lines seems to be linked to the primary component over the entire orbital cycle. Our CMFGEN analysis shows that this component has a larger mass-loss rate than the secondary, whilst its terminal velocity is assumed to be similar.

Fig. 7 Top: positions of the maximum of Hα emission as a function of the orbital phase. The red dots represent the positions of the maximum, whilst the blue line is the radial velocity curve of the primary component. Bottom: spectral variation of the Hα line as a function of the orbital phase. The emission vanishes at Φ = 0.00.

Upon inspection of the Hα profile in the observed spectra, it is clear that the maximum of emission of this line is associated to the primary on the whole orbit (see the upper panel of Fig. 7). A closer look at this profile at different phases (see the lower panel of Fig. 7) shows that at the periastron (Φ = 0.0) no emission is observed in the spectra, but this emission is present at Φ ~ 0.02 and even at Φ ~ 0.99. During this short time interval, when the two components are closest to each other, it appears that the emission generated by the two winds vanishes (see the X-ray analysis, Nazé et al. 2017).

The configuration of the system indicates that at the periastron the distance between the centers of both components is 71.2 R_⊙. Therefore, the distance that remains between the two photospheres (of about 16.9 R_⊙) is not large enough for the winds to strongly collide.

Fig. 8 Positions in the HR diagram of the two components of HD 166734. Top: evolutionary tracks (black lines) and isochrones (gray lines) from Ekström et al. (2012), computed with a rotational rate of 300 km s-1. Bottom: evolutionary tracks (black lines) and isochrones (gray lines) from Brott et al. (2011), computed with a rotational rate of 160 km s-1. The error bars represent the 1σ.

7. Discussion

7.1. Stellar evolution

The stellar luminosities determined from photometry as well as the effective temperatures derived from the CMFGEN modeling provide the exact locations of both stars in the Hertzsprung-Russell diagram. They are displayed in Fig. 8. Through our analysis, the orbital parameters of HD 166734 as well as the stellar parameters and the derived surface abundances of its two components suggest that they have evolved mainly as single stars (at least without Roche lobe overflow event). They are thus interesting objects to be compared with the current evolutionary tracks for massive stars.

The input physics and its implementation in six different evolutionary codes as well as their comparison with observed data were deeply discussed by Martins & Palacios (2013). In the present discussion, we compare the observational parameters of the two components of HD 166734 with the tracks of Geneva (Ekström et al. 2012) and with those of Brott et al. (2011; referred as STERN models by Martins & Palacios 2013, and hereafter in this paper).

7.1.1. Geneva models

The evolutionary tracks of Ekström et al. (2012) are computed with a rotational rate v_eq/v_crit of 0.4, a metallicity of Z = 0.014, an overshoot parameter of 0.1 and do not account for magnetism. The upper panel of Fig. 8 displays the evolutionary tracks and the isochrones of Geneva computed with rotation. Given the error bars on the luminosity and on the effective temperature, the two components are located, according to these tracks, close to the characteristic hook at the end of the main sequence. These positions therefore infer two different values for the initial mass of the two components. Indeed, within the error bars, the stars can be located either at the end of the main sequence, immediately before the stage where their temperature starts to increase, or before the yellow supergiant phase. Depending on their exact locations, the initial masses will be expected to be 44 or 47 M_⊙ for the primary and 37 or 39 M_⊙ for the secondary. According to the theoretical predictions, these stars would have ages between 5.2 and 5.4 Myr for the primary and between 5.8 and 6.0 Myr for the secondary (depending on the locations of both stars on the hook at the end of the main sequence). Furthermore, these tracks yield a radius of 27 R_⊙ for the primary, and of 26 R_⊙ for the secondary.

We also determine the present evolutionary masses of each component. For the primary we compute a mass of $\mathrm{34.}{\mathrm{5}}_{-9.5}^{\mathrm{+}\mathrm{10.2}}$ M_⊙ whilst for the secondary we obtain a mass of $\mathrm{31.}{\mathrm{0}}_{-9.5}^{\mathrm{+}\mathrm{10.5}}$ M_⊙. These masses can be compared to the spectroscopic masses computed from the luminosity, the effective temperature and the gravity of each star to be $\mathrm{37.}{\mathrm{7}}_{-16.1}^{\mathrm{+}\mathrm{29.2}}$ M_⊙ for the primary and $\mathrm{31.}{\mathrm{8}}_{-14.4}^{\mathrm{+}\mathrm{26.6}}$ M_⊙ for the secondary as well as to the dynamical masses of $\mathrm{39.}{\mathrm{5}}_{-4.4}^{\mathrm{+}\mathrm{5.4}}$ M_⊙ and $\mathrm{33.}{\mathrm{5}}_{-3.7}^{\mathrm{+}\mathrm{4.6}}$ M_⊙, for the primary and the secondary, respectively. The large rotational rate (v/v_crit = 0.4) is however not realistic for massive stars in the Milky Way (Simón-Díaz & Herrero 2014). This infers an overestimation of the ages of the stars, and an underestimation of their evolutionary masses. From the tracks of Geneva, the two components of HD 166734 are thus located close to the terminal age main sequence (TAMS) but given the fast evolution at the end of the main sequence, this part of the HR diagram is expected to be empty (McEvoy et al. 2015). The inclusion of extra mixing (either by tidal effects at periastron or by increasing the overshoot parameter) in the evolutionary models would allow to increase the lifetime of the stars on the main sequence. Even though the initial rotational rate is large for the Geneva models, the spin down is quick. Indeed, for our individual stars taken as single, the actual equatorial velocities are predicted to be 15 km s^-1 for the primary and 18 km s^-1 for the secondary, which is too small with respect to the rotational velocities of our stars. The spin down is notably ruled by the mass-loss rate included in the models, suggesting that the predictions concerning the mass-loss in the Geneva models seems to be overestimated or that another phenomenon is at work.

The initial rotational velocity of about 300 km s^-1 seems therefore too large for both components. We thus quickly compare the positions of the stars to the Geneva evolutionary tracks computed without rotation. From these models, the stars are on the main sequence, with estimated ages of about 3.0 Myr for the primary and 3.4 Myr for the secondary. The primary has an initial mass of 59 M_⊙ and an actual mass of 51.4 M_⊙ whilst the secondary has an initial mass of 50 M_⊙ and an actual mass of 45.3 M_⊙; thus larger than what we found for the two components from evolutionary tracks with rotation as expected from the theory.

7.1.2. STERN models

We decided to compare our observational values with the STERN evolutionary models (Brott et al. 2011). These models constitute a grid with a sufficiently fine mesh, attenuating the errors due to the interpolation. This contrasts with the Geneva models for which one rotational rate is taken into account. The STERN models are computed with a metallicity of Z = 0.0088, an overshoot parameter of 0.335 and the rotation is taken as a free parameter. The BONNSAI facility² (Schneider et al. 2014) was used as a Bayesian approach to derive the distribution of stellar parameters for each component of HD 166734. The best fits provide an initial rotational velocity of 150 km s^-1 for both objects, giving an equatorial rotational velocity of 110 km s^-1 at their current stage (in adequation with our observational parameters). The lower panel of Fig. 8 displays the locations of the stars on the HR diagram build from the STERN evolutionary tracks and the isochrones computed with an initial rotational velocity of 160 km s^-1. For the primary, the initial mass is expected to be 56.1$\begin{array}{c}\mathrm{+}\mathrm{4.7}\\ -5.4\end{array}$ M_⊙, giving an evolutionary mass of 47.8$\begin{array}{c}\mathrm{+}\mathrm{5.7}\\ -4.4\end{array}$ M_⊙. For the secondary, the initial mass is 47.4$\begin{array}{c}\mathrm{+}\mathrm{4.0}\\ -4.6\end{array}$ M_⊙ and the evolutionary mass is 41.2$\begin{array}{c}\mathrm{+}\mathrm{4.7}\\ -3.5\end{array}$ M_⊙. The STERN models provide an age of about 3.0 Myr for the primary and of about 3.5 Myr for the secondary; in agreement with the non-rotating tracks of Geneva. As emphasized by Martins & Palacios (2013), the characteristic hook at the end of the main sequence occurs at lower temperatures because of the large amount of overshooting. These authors also mentioned that an overshooting parameter of 0.335, as used for the computation of the STERN models, is too large for stars with masses above 20 M_⊙. Therefore, the ages estimated from these models are probably underestimated and their evolutionary masses overestimated. However, the STERN tracks better reproduce the rotational velocities of the supergiants.

7.1.3. Surface abundances

Whatever the evolutionary tracks we use, the secondary component appears more evolved than the primary. An enhanced mixing through tides at periastron could explain this difference in age (see Sect. 7.2). As we already mentioned in Sect. 5, the primary and the secondary exhibit an overabundance in helium and in nitrogen as well as a depletion in carbon. Furthermore, the primary is also depleted in oxygen. The abundances thus confirm the evolved status of each object. However, they do not differ from abundances of single evolved massive stars. When we compare the surface abundances of the two components of HD 166734 to the values of only the supergiant stars determined by Martins et al. (2015) in the framework of the MiMeS project and to the values derived by Bouret et al. (2012) in their analysis of a sample of Galactic supergiant stars, we see a clear correlation between N/O and N/C (see the upper panel of Fig. 9) that does not differ from the results of Bouret et al. (2012) or of Martins et al. (2015), represented in the upper panel of Fig. 9 by black dots. This agrees with the expectations of nucleosynthesis through the CNO cycle. The ratios N/C and N/O follow those of single stars according to the evolutionary models, within the uncertainties. The lower panel of Fig. 9 shows that the ratio N/C increases as surface gravity decreases, as was observed for single supergiant Galactic massive stars (Bouret et al. 2012; Martins et al. 2015). The ratios of N/C and N/O determined from the two components of HD 166734 are thus correlated and in the same range as those measured for presumably single supergiant objects, meaning that no or at least few interactions occur between these two components.

Fig. 9 Top:log  (N/C) (by number) as a function of log  (N/O). The expected trends for the case of CN and CNO equilibrium are shown by the dashed lines. Bottom: log  (N/C) (by number) as a function of log g. The two components of HD 166734 are represented by red dots, whilst the sample of Galactic supergiant stars analyzed by Bouret et al. (2012) and Martins et al. (2015) is represented by black dots. The evolutionary tracks are from Ekström et al. (2012) and are computed for a rotational rate (initial over critical) of 0.4. We use evolutionary tracks with initial masses larger than 10 M⊙.

7.2. Synchronization of the system

HD 166734 is a highly-eccentric binary system where the secondary appears to be more evolved than the primary. This problem linked to the ages can be due to tidal mixing at periastron. This tidal mixing would affect the secondary more than the primary because of its size. In addition, some matter could be transferred to the primary with the consequence of increasing its spin, inducing more mixing and therefore a higher temperature; this increase of the temperature would modify the apparent age of the primary by becoming younger. This effect should be however limited because of the configuration of the system.

We then may wonder whether these tidal interactions (that would be responsible for this tidal mixing) can affect the rotations of the components and can possibly synchronize the system (tidal interactions can indeed be responsible for a long timescale evolution of the binary orbital elements and the stellar rotation). To characterize these interactions in HD 166734, we focus on the theory of Hut (1981, 1982). If we compute the ratio between the rotational and the orbital angular momentum (Eq. (14) in Hut 1981), we find for the present stars of HD 166734 that a large part of the angular momentum is contained in the orbit. An important classification parameter (α) representing the ratio of the orbital angular momentum to the rotational angular momentum at equilibrium was introduced by Hut (1981). According to the value of this parameter, the circularization will take much more time than the synchronization. In the case of HD 166734, for which the orbital separation between the two components is relatively large, the α parameter is extremely large, meaning that the circularization time is expected to be longer than the synchronization time. Always according to Hut’s theory, the eccentricity of the orbit should decrease very slowly as well as the semimajor axis (de/ dt and da/ dt < 0) whilst the rotational angular velocities of both stars should increase as a function of time (dΩ/dt> 0) but all this evolution is expected on timescales longer than the current age of the system. In eccentric binaries, the term of synchronization between the rotation and the revolution is rather replaced by the pseudo-synchronization at the moment of the periastron. This pseudo-synchronization of the system is reached when the rotational angular velocity of a star is at 0.8 times the instantaneous orbital angular velocity at periastron (see Hut 1981). Indeed, assuming that the rotational axes of both components are perpendicular to the orbital plane, we use the combination of the inclination (63.0° ± 2.7°) of the system, the equatorial rotational velocities (107 and 110 km s^-1) and the radii (27.5 and 26.8 R_⊙) of both components, obtained in Sect. 4, to compute rotational periods of 13.0 days for the primary and of 12.3 days for the secondary. Both stars appear to have their rotations almost synchronized with each other but these values are smaller than the orbital period. When we compare the rotational angular velocities with the instantaneous orbital angular velocity at periastron, we compute ratios of about 0.6 for the primary and secondary components. We therefore conclude that pseudo-synchronization (Hut 1981) has not yet been achieved for HD 166734.

7.3. Distance to HD 166734

In Sect. 4, we determined the bolometric magnitudes of −9.85 ± 0.23 for the primary and −9.58 ± 0.26 for the secondary, giving stellar luminosities of $\log \frac{\mathit{L}}{{\mathit{L}}_{\mathrm{\odot }}}\mathrm{=}\mathrm{5.840}\mathrm{\pm }\mathrm{0.092}$ and $\log \frac{\mathit{L}}{{\mathit{L}}_{\mathrm{\odot }}}\mathrm{=}\mathrm{5.732}\mathrm{\pm }\mathrm{0.104}$, respectively. Given their spectral types and their effective temperatures, we can compute, from Martins & Plez (2006), bolometric corrections of −2.97 ± 0.09 and −2.83 ± 0.10 for the primary and the secondary, respectively. From these values, we compute absolute magnitudes of M_V = −6.88 for the primary and M_V = −6.75 for the secondary.

From the observed V magnitude of HD 166734 of 8.42 and its B−V of 1.07 (Ducati 2002), accounting for (B−V)₀ = −0.27 (Martins & Plez 2006), and a brightness ratio of 1.28 between the primary and the secondary components, we estimate, for HD 166734, a distance modulus of 11.85 and thus a distance of 2.34 ± 0.30 kpc, confirming the distance estimated by Conti et al. (1980).

Fig. 10 Section of the FEROS spectrum of HD 166734 showing the Na i interstellar lines, covering a velocity range of −55 to + 35 km s-1. The dotted lines indicate, from left to right, RVs of −36, −29, −13, + 4 km s-1 for both panels.

Next, we compare this distance with the kinematic distance estimated from the Galactic rotation curve. In Fig. 10, we show the Na i absorption line profiles. The Na i lines are relatively saturated on a width ranging from + 4 to −13 km s^-1 and with a shallower absorption profile between −29 and −36 km s^-1. From the Galactic rotation model of Fich et al. (1989), this velocity range gives a mean distance of 2.38 ± 0.17 kpc. This distance is in agreement with the distances estimated from the stellar luminosities or derived by Conti et al. (1980).

8. Conclusions

We present the results of the spectroscopic and photometric campaigns devoted to HD 166734. The spectra taken close to the periastron passage have allowed us to refine the orbital solution of the system, removing the existing doubts on the masses of both components. Our analysis has indeed shown that HD 166734 is a binary system with an orbital period of 34.537723 days and an eccentricity of 0.618. From the photometry, we determine an inclination of 63.0° ± 2.7° for the system and actual masses of $\mathrm{39.}{\mathrm{5}}_{-4.4}^{\mathrm{+}\mathrm{5.4}}$ M_⊙ and $\mathrm{33.}{\mathrm{5}}_{-3.7}^{\mathrm{+}\mathrm{4.6}}$ M_⊙ for the primary and the secondary.

The individual spectra of each component, obtained by spectral disentangling, indicate a primary with a spectral classification of O7If, an effective temperature of 32 000 ± 1000 K, a log g of 3.15 ± 0.10 and a luminosity of $\log \frac{\mathit{L}}{{\mathit{L}}_{\mathrm{\odot }}}\mathrm{=}\mathrm{5.840}\mathrm{\pm }\mathrm{0.092}$. We found that the secondary has the spectral type O9I(f) and its stellar parameters are T_eff = 30 500 ± 1000 K, log g = 3.10 ± 0.10, and $\log \frac{\mathit{L}}{{\mathit{L}}_{\mathrm{\odot }}}\mathrm{=}\mathrm{5.732}\mathrm{\pm }\mathrm{0.104}$. We thus conclude, unlike Conti et al. (1980), that the hottest and the most luminous star is also the most massive one. Both components show modified helium, carbon, nitrogen and oxygen surface abundances. They are both enriched in helium and nitrogen and are depleted in carbon. In addition, the primary is also depleted in oxygen. Their surface abundances do not however differ from the predictions of the “classical” evolutionary paths for single stars. This means that no or, at most, few interactions occur between the two stars, making them perfectly suitable as test objects for evolutionary models.

Systems containing a pair of massive evolved O-type stars are relatively rare. Despite their number, they represent excellent tests for atmospheric and evolutionary models. Analyses of the kind we did for HD 166734 are thus crucial to validate these models. In this context, we may quote other systems that can help with this validation. An interesting system is R139 (Taylor et al. 2011) composed of an O6Iafc primary and of an O6Iaf secondary on a wide eccentric orbit. This system is the most massive binary system known to host two evolved Of supergiants but the physical parameters of each component are not complete enough to allow a deep comparison with the evolutionary models. Finally, WR21a (Tramper et al. 2016; Gosset & Nazé 2016) is a very massive system composed of two younger O3/WN5 and O3 stars. Although this system is less evolved, it has orbital period (P = 31.67 days) and eccentricity (e = 0.64) comparable to HD 166734.

In the present paper, we compare the fundamental parameters of each component with evolutionary tracks of Geneva (Ekström et al. 2012) and of STERN (Brott et al. 2011). We conclude from our observations that the tracks of Geneva have an overly large initial rotational rate, that shortens the lifetime of the star on the main sequence and it seems that the mass-loss used to compute the tracks is also too large to reproduce the actual rotational velocities of the stars. At the same time, the STERN tracks have an overly large overshoot parameter, but the spin down of the stars is sufficient to agree with the actual rotational velocities of the stars. Finally, the ages and the evolutionary masses for both components are different as a function of the tracks used in the analysis. Even though other systems (such as those quoted above) must still be investigated, HD 166734 allowed us to show that the evolutionary models are not reliable enough to reproduce the exact characteristics of supergiant stars.

Acknowledgments

The authors thank the anonymous referee for his/her constructive remarks on the paper. This research was supported by the Fonds National de la Recherche Scientifique (F.R.S.-F.N.R.S.), by the PRODEX XMM contract (Belspo), and through the ARC grant for Concerted Research Actions, financed by the French Community of Belgium (Wallonia-Brussels Federation). L.M. thanks Pavel Mayer for his advice in the light curve modeling and John Hillier for making CMFGEN available. We are grateful to the staff of La Silla ESO Observatory for their technical support.

References

Appendix A: Journal of the spectroscopic observations of HD 166734

Appendix B: Periodograms

Fig. B.1 Top: power spectrum computed on the basis of RVS–RVP. Bottom: spectral window computed from the spectroscopic observations.

Fig. B.2 Top: power spectrum of the R light curve of HD 166734. The diamond indicates the location of the orbital frequency while the dots represent its harmonics. Bottom: spectral window computed from the R light curve.

Fig. B.3 Top: power spectrum of the I light curve of HD 166734. The diamond indicates the location of the orbital frequency while the dots represent its harmonics. Bottom: spectral window computed from the I light curve.

All Tables

All Figures

	Fig. 1 FEROS spectra in the region of He i 4471 (left) and C iii 5696 (right) shifted vertically as a function of phase (times 3 for clarity) from 0.0 at the bottom to 1.0 at the top.
In the text

	Fig. 2 Radial velocity curves of HD 166734 computed with the RVs from the present paper (provided electronically). The primary is represented by filled symbols whilst the open ones indicate the RVs of the secondary. Circles represent the RVs measured on the FEROS spectra and diamonds indicate the RVs measured on the ESPRESSO data. Φ = 0.0 corresponds to the periastron passage.
In the text

	Fig. 3 Top: power spectrum of the V light curve of HD 166734. The diamond indicates the location of the orbital frequency while the dots represent its harmonics. Bottom: spectral window computed from the V light curve.
In the text

	Fig. 4 Light curves of HD 166734 obtained with TAROT in the V (top panel), R (second panel) and I (third panel) filters. The red lines represent the best fit computed with PHOEBE. The magnitude outside the eclipse has been set to 0 in the three filters. The representation shown at the bottom is a zoom-in on the eclipse observed in the V light curve (plotted at the top of this figure).
In the text

	Fig. 5 Best CMFGEN fit (red) of the disentangled optical spectrum (black) of the primary component of HD 166734.
In the text

	Fig. 6 Same as Fig. 5 but for the secondary component of HD 166734.
In the text

	Fig. 7 Top: positions of the maximum of Hα emission as a function of the orbital phase. The red dots represent the positions of the maximum, whilst the blue line is the radial velocity curve of the primary component. Bottom: spectral variation of the Hα line as a function of the orbital phase. The emission vanishes at Φ = 0.00.
In the text

	Fig. 8 Positions in the HR diagram of the two components of HD 166734. Top: evolutionary tracks (black lines) and isochrones (gray lines) from Ekström et al. (2012), computed with a rotational rate of 300 km s-1. Bottom: evolutionary tracks (black lines) and isochrones (gray lines) from Brott et al. (2011), computed with a rotational rate of 160 km s-1. The error bars represent the 1σ.
In the text

Fig. 9 Top:log  (N/C) (by number) as a function of log  (N/O). The expected trends for the case of CN and CNO equilibrium are shown by the dashed lines. Bottom: log  (N/C) (by number) as a function of log g. The two components of HD 166734 are represented by red dots, whilst the sample of Galactic supergiant stars analyzed by Bouret et al. (2012) and Martins et al. (2015) is represented by black dots. The evolutionary tracks are from Ekström et al. (2012) and are computed for a rotational rate (initial over critical) of 0.4. We use evolutionary tracks with initial masses larger than 10 M⊙.

In the text

	Fig. 10 Section of the FEROS spectrum of HD 166734 showing the Na i interstellar lines, covering a velocity range of −55 to + 35 km s-1. The dotted lines indicate, from left to right, RVs of −36, −29, −13, + 4 km s-1 for both panels.
In the text

	Fig. B.1 Top: power spectrum computed on the basis of RVS–RVP. Bottom: spectral window computed from the spectroscopic observations.
In the text

	Fig. B.2 Top: power spectrum of the R light curve of HD 166734. The diamond indicates the location of the orbital frequency while the dots represent its harmonics. Bottom: spectral window computed from the R light curve.
In the text

	Fig. B.3 Top: power spectrum of the I light curve of HD 166734. The diamond indicates the location of the orbital frequency while the dots represent its harmonics. Bottom: spectral window computed from the I light curve.
In the text