© ESO, 2012

1. Introduction

NGC 6193 is a young open cluster surrounded by an H ii region in the Ara OB1 association. The core of the cluster hosts two bright and massive stars separated by about 10″: HD 150135 and HD 150136.

Several investigations in different wavelength domains (optical, radio, and X-rays) allowed the knowledge of HD 150136 to be improved. In the optical domain, Niemela & Gamen (2005) report this star as a binary system composed of an O3 primary and an O6 secondary. Indeed, spectral absorption lines of N iv λ5203 and N v λλ4603 − 19 (wavelengths are hereafter expressed in Å), as well as the emission line N iv λ4058, were observed in the composite spectra. These lines are considered to be typical features for O2–O3 stars (Walborn et al. 2002). This is a rather rare category in which only a few of these stars are reported as binary systems. Therefore, even though these objects are considered, with the WNLh stars, as possibly belonging to the most massive objects, their actual masses are scarcely known.

Located at a distance of 1.32    ±    0.12 kpc (Herbst & Havlen 1977), HD 150136 would thus harbour the nearest O3 star and would constitute a target of choice for investigating the fundamental parameters of such a star. Initially, the orbital period of the system was estimated to about 2.7 d whilst the mass ratio between the primary and the secondary stars was measured to be equal to 1.8 (unpublished, see Garmany et al. 1980). The analysis of Niemela & Gamen (2005) allowed the authors to refine these values to 2.662 d for the orbital period and 1.48 for the mass ratio. Moreover, the possibility of a third component bound to the system was envisaged (though no clear feature was detected) to explain that the He i lines in the spectra did not follow the orbital motion of any of the binary components.

In the radio domain, HD 150136 was reported as a non-thermal radio emitter (Benaglia et al. 2006; De Becker 2007), which suggests there are relativistic electrons accelerated in shocks issued from a wind-collision region. This non-thermal source could be of about 2″ (i.e., the resolution, Benaglia et al. 2006) and could be associated with a close visual component separated from HD 150136 by about 1.6″ (Mason et al. 1998), i.e., about 0.01 pc if we assume a distance of about 1.3 kpc. So far, there has however been no indication that this star is physically bound to the binary system.

The Chandra X-ray observations of the young open cluster NGC 6193 revealed that HD 150136 is one of the X-ray brightest O-type stars known (log L_X (erg s^-1) = 33.39 and L_X/L_bol = 10^-6.4, Skinner et al. 2005). Its emission probably results from a radiative colliding wind interaction and appears to be slightly variable on a timescale smaller than one day. According to Skinner et al. (2005), the discrepancy between the orbital period and the timescale of the variation in this X-ray emission is probably an occultation effect either linked to the stellar rotation or to a putative third star with a rather short period.

All these investigations thus suggest that HD 150136 is a multiple system rather than a binary one, even though evidence is still lacking. In this context, an intense spectroscopic monitoring at high resolution in the optical domain, first initiated to study the binary itself, was performed to bring new light on this system. We organize the present paper as follows. Section 2 describes the observations and the data reduction technique. Section 3 reveals the existence of a third star in the data of HD 150136 and focuses on the disentangling process in order to compute the individual spectra of the components. A first determination of the orbital solution of the known short-period system based on an SB2 configuration is presented in Sect. 4. We also demonstrate that these parameters can be optimized by assuming the external influence of a third component. The possible orbit of this third star is addressed in Sect. 5 and an SB3 partial solution derived. This discovery makes HD 150136 one of the first systems where it is demonstrated that the third component is physically bound to the inner system (see also Tr16-104, Rauw et al. 2001; and HD 167971, Blomme et al. 2007). However, it is the first time that the primary is an O 3 star. We thus in Sect. 6 use the CMFGEN atmosphere code (Hillier & Miller 1998) to derive the fundamental properties of the stars to improve our knowledge on the O 3 star parameters. Moreover, the configuration of the short-period binary seems to suggest a wind interaction zone possibly detectable through the He ii λ4686 and Hα lines. We devote Sect. 7 to analysis of the variability observed in these lines by applying a Doppler tomography technique. Finally, we discuss the implications of our results in Sect. 8 and we provide the conclusions in Sect. 9.

2. Observations and data reduction

A multi-epoch campaign was devoted to observing HD 150136. Between 1999 and 2009, we collected 79 high-resolution, high signal-to-noise spectra of this star with the échelle Fiber-fed Extended Range Optical Spectrograph (FEROS) successively mounted at the ESO-1.5 m (for spectra taken in 1999 and 2001) and the ESO/MPG-2.2 m (for observations between 2002 and 2009) telescopes at La Silla (Chile). This instrument has a resolving power of about 48 000 and the detector was a 2 k  ×  4 k EEV CCD with a pixel size of 15 μm  ×  15 μm. Exposure times were between 2 and 10 min to yield a signal-to-noise ratio over 200, measured in the continuum on the region [4800–4825] Å. FEROS is providing 39 orders allowing the entire optical wavelength domain, going from 3800 to 9200 Å, to be covered. The journal of the observations is listed in Table 1 (including the individual signal-to-noise ratios). For the data reduction process, we used an improved version of the FEROS pipeline working under the MIDAS environment (Sana et al. 2006). The data normalization was performed by fitting polynomials of degree 4–5 to carefully chosen continuum windows. We mainly worked on the individual orders, but the regions around the Si iv λλ4089, 4116, He ii λ4686, and Hα emission lines were normalized on the merged spectrum. We also note that the flatfields of the 2009 campaign presented some oscillations of unknown origin. These oscillations represent about 3% to 5% of the continuum flux. We applied a filter to smooth these flatfields and thus to limit the effects of the undesirable oscillations on the observed spectra.

Equivalent widths (EWs) and the Doppler shifts were estimated on each spectrum by fitting Gaussian profiles on the spectral lines. For the computation of the radial velocities, we adopted the rest wavelengths from Conti et al. (1977), for the lines with a rest wavelength shorter than 5000 Å, and Underhill (1994), for the other ones.

3. Third component and disentangling

On the basis of the high-resolution spectra analysed in the present paper, we report the unambiguous detection of a third star (T) signature for the system HD 150136 (Fig. 1). The data indeed exhibit a clear variation in the third component compared to the rest wavelength (He i λ5875.62) shown by the dashed line. Moreover, the last observation confirms the periodicity of the variation, suggesting that the third star is also a member of a binary system. Although several analyses have suspected the existence of this component, it has never been unveiled until now.

To separate the contribution of each star in the observed spectra, we apply a modified version of the disentangling program of González & Levato (2006) by adapting it to triple systems. This method consists of an iterative procedure that alternatively uses the spectrum of a first component (shifted by its radial velocity) and then the spectrum of a second component (also shifted by its radial velocity) to successively remove them from the observed spectra. The result is a mean spectrum of the remaining third star. This program also allows us to recompute the radial velocities of each star with a cross-correlation technique, even at phases for which the lines are heavily blended. The cross-correlation windows are built on a common basis that gathers the positions of the He ii λ4200, He i λ4471, He ii λ4542, He i λ4713, He i λ4921, O iii λ5592, C iv λλ5801−12, and He i λ5876 spectral lines. In addition, we add the N v λλ4603−19 lines for the primary star. The radial velocities obtained with this technique are listed in Table 1. The resulting spectra are disentangled in a [4000–5900] Å wavelength domain and at Hα. Though the He ii λ4686 and Hα lines present particular variations along the orbital phase, we include them in the process. However, the assumption of constant lines (on which the disentangling method relies) is no longer satisfied, thus making these line profiles unreliable.

To obtain the final individual spectra, we have to correct the disentangled spectra by the brightness ratio. For that purpose, we use the disentangled and the observed spectra to measure the EWs and constrain the spectral classification of the three components of HD 150136. This determination is based on the quantitative classification criteria presented by Conti & Alschuler (1971) and Conti (1973b), which uses the EW ratios of the He i λ4471 and He ii λ4542 lines to define the spectral type. This leads to an O3–O3.5 primary star, an O5.5–O6 secondary star and an O6.5–O7 third star. The weak N iii λλ4634 − 41 emission and the strong He ii λ4686 absorption lines of the secondary and the tertiary components suggest adding the ((f)) suffix to their spectral types. In addition, the criterion given by Walborn et al. (2002) confirms the classification of the primary. Indeed, its spectrum is reminiscent of that of HD 64568, i.e., close to an O3V((f^∗)). However, the N iv λ4058 line has an intermediate EW (~160 mÅ) between the values measured on the N iii λ4634 (~140 mÅ) and N iii λ4641 (~250 mÅ) lines, thereby indicating a spectral type between O3V((f^∗)) and O3.5V((f⁺)). The ((f^∗)) reports a spectrum with the N iv λ4058 emission line stronger than the N iii λλ4634 − 41 lines and a weak He ii λ4686 absorption line whilst the ((f⁺)) refers to a spectrum with medium N iii λλ4634 − 41 emission lines, the weak He ii λ4686 line, and the Si iv λλ4089 − 4116 lines in emission. To compute the different brightness ratios required to determine the relative luminosities of the stars, we compare the EWs measured on the observed spectra with the theoretical values issued from synthetic models whose spectral classification is the same as for our disentangled spectra. We stress that these models are computed from solar surface abundances (except for the nitrogen abundance). To avoid a possible bias due to differences between fits and models, we also focus on the tables of Conti & Alschuler (1971), Conti (1973a), and Conti & Leep (1974) for the secondary and the tertiary components (the spectral type of the primary is not available in these tables). This comparison between the “observed” EWs and the “canonical” EWs for isolated single stars is performed on the helium and oxygen lines, yielding brightness ratios of l_P/l_S = 2.54, l_S/l_T = 1.51, and l_P/l_T = 3.87. An error bar at 1-σ of about 0.14 is estimated between the “observed/canonical” ratios measured from the different spectral lines. The individual spectra, corrected from the brightness ratios, are shown as black lines in Figs. 2 to 4.

4. Inner short-period binary system

With the radial velocities refined by the disentangling program, we compute the orbital solution of the short-period system formed by the primary (P) and secondary (S) components. We first re-estimate the orbital period of the short-term binary by using the Heck-Manfroid-Mersch technique (hereafter HMM, Heck et al. 1985; as revised by Gosset et al. 2001) on the basis of RV_S − RV_P values. We obtain an orbital period of about 2.67456  ±  0.00001. We then use the Liège Orbital Solution Package (LOSP¹) to determine the SB2 orbital solution of the inner binary. Figure 5 displays the radial-velocity curves, and Table 2 lists the orbital parameters determined for an SB2 configuration. The system shows no significant deviation (test of Lucy & Sweeney 1971) from circularity, so we thus adopt e = 0. As a consequence, T₀ represents the phase of the conjunction with the primary in front. We see that the dispersion of some observational points around the theoretical radial-velocity curves is rather large, which is probably due to the influence of the third star.

Consequently, we compute the averages, for the primary and the secondary stars over each epoch, of the differences between the observations and the fitted radial-velocity curves. We also add to these values the constant systemic velocities quoted in Table 2. Moreover, we compare these results to the systemic velocities of the primary and the secondary components computed on each individual observing campaign. Both methods report comparable global trends of the velocities for the inner system. The adopted values are reported in Table 3. By plotting the mean systemic velocities of the inner system (see Table 3) and the mean radial velocities of the third component computed on each epoch as a function of time (Fig. 6), we detect a clear trend toward a periodical motion that is reminiscent of an additional orbit. Therefore, we suggest that HD 150136 is at least an SB3 system.

5. The outer system

5.1. The outer system as an SB2

We use the mean systemic velocities of the inner system determined from each epoch as a primary component, and we consider the third component as the secondary companion. These velocities do not allow us to determine the real orbital period of the wide outer system with accuracy. However, we can estimate this value between roughly 2950 and 5500 d, i.e., between 8 and 15 years. Consequently, computing an orbital solution is useless without further high-resolution monitoring of HD 150136. Nonetheless, the measured velocities can provide information about the orbit, as well as on the mass ratio between the inner system and the third component. To estimate the mass ratio, we calculate the average between the systemic velocities of the primary and the secondary stars and we take these values as radial velocities of the inner system (P + S). We focus on the radial velocities of both outer system components by fitting the following relation by a least-square method: $\mathit{R}{\mathit{V}}_{\mathrm{T}}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}{\mathit{c}}_{\mathrm{1}}\hspace{0.17em}\mathit{R}{\mathit{V}}_{\mathrm{P}\mathrm{+}\mathrm{S}}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{+}{\mathit{c}}_{\mathrm{2}}$where RV_P + S(φ), RV_T(φ) are the systemic radial velocities of the inner system and the radial velocities of the third component, respectively, whereas ${\mathit{c}}_{\mathrm{1}}\mathrm{=}\mathrm{-}\frac{{\mathit{M}}_{\mathrm{P}\mathrm{+}\mathrm{S}}}{{\mathit{M}}_{\mathrm{T}}}$ with M_P + S and M_T the masses of the inner system and the third star, respectively. The linear regression gives us values of about c₁ =  −2.60 and c₂ =  −79.58 km s^-1. From c₂, we deduce γ_T ≈  −20 km s^-1. To verify the impact of the period on the value of the mass ratio, we have also computed several orbital solutions by assuming different orbital periods sampled between 2950 and 5500 d. All these solutions led to a mass ratio of about 2.94. However, among the velocities exhibited in Fig. 6, the point corresponding to HJD = 2 452 037.8496 (i.e., to the observing campaign made in 2001) displays a clear deviation probably due to the shortage of observations during this campaign. If we remove this point, we obtain by the linear-regression method a mass ratio of about 3.1 and by the other method a mass ratio of about 3.2. We thus conclude on a mass ratio of about 2.9 ± 0.3 between the inner system and the third star. Similarly, from this sample of orbital periods, we derive an eccentricity higher than 0.3 for the outer system.

5.2. HD 150136 as a triple system

To better constrain the orbital parameters of the inner system, we take the influence of the third component into account by fitting all the components together. For that purpose, we determine with another method the orbital parameters of the inner system. This technique provides the approximate orbital solutions for the two SB2 systems in HD 150136, i.e., the inner system composed of the primary and the secondary stars and the outer one represented by the mean systemic velocities characterizing the centre of mass of the inner system and the radial velocities of the third component. Therefore, we perform on all the three sets of data (P, S, and T) a global least square fit with an assumed circular orbit for the inner binary. For both systems, we follow the same procedure. We first deduce the approached orbital solution using the Keplerian periodogram introduced by Zechmeister & Kürster (2009). This method consists in a fit of a simple sine function in the true anomaly domain. The approximate orbital solution is then globally refined using a Levenberg-Marquardt minimization. We search the period in [0.6T, 2.0T] for the outer orbit where T is the total time span of the observations. We use a constant step of frequency for both systems equal to 0.01/T. Figure 7 exhibits the orbital solution for the short-term binary system by considering the influence of the third star. We report in Table 2 the refined orbital parameters (SB3 solution). Among these values, we notably observe a decrease in the error bars. Following an F-test on the ratio of the two χ² of the fits (SB2 and SB3 solutions, see Table 2), we conclude that this solution is clearly improved in comparison with the one computed in Sect. 4 (see Fig. 5). The improvement is marked enough to remain significant even though the F-test is not fully valid due to the non-linearity.

6. Atmosphere analysis of the three stars

The three components of HD 150136 are then modelled separately by using the CMFGEN atmosphere code (Hillier & Miller 1998). Non-LTE models including winds and line-blanketing are computed. The computations are done in spherical geometry and in the co-moving frame. A more detailed description of this code is provided in Hillier & Miller (1998). For the velocity structure, we use TLUSTY hydrostatic structures in the photosphere associated with a β-velocity law v = v_∞(1 − R/r)^β for the wind part (v_∞ being the wind terminal velocity). The TLUSTY models are taken from the OSTAR2002 grid of models computed by Lanz & Hubeny (2003) and β = 0.8, typical for O dwarfs, is adopted (Repolust et al. 2004). The models include the following chemical species: H, He, C, N, O, Ne, Mg, Si, S, Ar, Fe, and Ni with the solar composition of Grevesse et al. (2007), except for N where the abundance is kept free. This is also the case for C and O abundances, but we see below that they are not well constrained and so are fixed to the solar values. We also use the superlevel approach to limit the required computation memory. When the model atmosphere has converged, the synthetic spectrum is generated through the formal solution of the radiative transfer equation. For that purpose, a microturbulent velocity varying linearly (with velocity) from 10 km s^-1, in the photosphere, to 0.1  ×  v_∞, in the outer boundary, is used. The emerging spectra are compared to the disentangled ones, corrected for the brightness ratios, in order to constrain the main parameters more accurately.

We thus determine the effective temperature (T_eff), the log g, the luminosity, the N content, the projected rotational velocity, and the macroturbulence of each component with the following diagnostics:

• the effective temperature: T_eff is computed on the basis of the ratio between the He i λ4471 and He ii λ4542 lines. However, other lines of helium such as He i λ4026, He ii λ4200, He i λ4388, He i λ4713, He i λ4921, and He i λ5876 are also considered as indicators. Generally, we estimate the error in the models on the determination of that parameter to be equal to 1000 K;

• the gravity: we use the wings of the Balmer lines (Hδ, Hγ, and Hβ) as main diagnostics to constrain the log g. We refrain from using the Hα line because it forms, at least partially, in the stellar wind. For the tertiary star, the determination of the gravity is more uncertain because of the disentangling process (see Sect. 8). Therefore, the uncertainty on the log g is equal to 0.1 dex for the primary and the secondary and to 0.25 dex for the third component;

• the luminosity: the luminosity is constrained by the V-band magnitude by adopting a distance of 1.32 kpc and a colour excess E(B − V) of 0.44 as explained in Sect. 3. The accuracy on the luminosity depends on the uncertainty on the distance of HD 150136;

• the surface abundances: the nitrogen abundance is determined on the N iii triplet in the 4500 − 4520 Å range but also on the N iv λ4058 and N iii λλ4634 − 41 line profiles. The typical uncertainty on the CNO surface abundances is of about 50%. To estimate these uncertainties, we run models with several different values of N and compare the corresponding spectra to the observed line profiles. When a clear discrepancy is seen on all lines of the same element, we adopt the corresponding abundance as the maximum/minimum value of the possible chemical composition. Carbon abundance is estimated to be solar for the three stars, whilst oxygen abundance is not determined because of the lack of diagnostic lines;

• the rotational and macroturbulent velocity: the synthetic spectra are successively convolved by a rotational profile and then by a Gaussian profile representing the macroturbulence. First, the projected rotational velocity is computed with the method of Simón-Díaz & Herrero (2007) based on the Fourier transform. Then, the macroturbulent velocity is determined so that the shape of the He i λ4713 line is reproduced best. When this line is not visible in the spectrum, we use the He i λ5876 and/or the C iv doublet λλ5801 − 12 as second indicator.

  Fig. 8 Left: the observed He ii λ4686 line profiles in the HD 150136 spectrum. Right: the restored He ii λ4686 line profiles in the HD 150136 spectrum when corrected for the absorption lines.

From V = 5.5 ± 0.12 (Niemela & Gamen 2005), a mean (B − V) = 0.16 ± 0.01, measured from Reed’s catalogue (Reed 2005), and a (B − V)₀ =  −0.28 (Martins & Plez 2006), we compute an absolute V magnitude for the entire system equal to ${\mathit{M}}_{\mathit{V}}\mathrm{=}\mathrm{-}\mathrm{6.4}{\mathrm{7}}_{-0.36}^{\mathrm{+}\mathrm{0.34}}$ by assuming a distance of about 1.32 ± 0.12 kpc. From the different brightness ratios, we then obtain ${\mathit{M}}_{{\mathit{V}}_{\mathrm{P}}}\mathrm{=}\mathrm{-}\mathrm{5.9}{\mathrm{1}}_{-0.36}^{\mathrm{+}\mathrm{0.35}}$, ${\mathit{M}}_{{\mathit{V}}_{\mathrm{S}}}\mathrm{=}\mathrm{-}\mathrm{4.9}{\mathrm{0}}_{-0.36}^{\mathrm{+}\mathrm{0.35}}$ and ${\mathit{M}}_{{\mathit{V}}_{\mathrm{T}}}\mathrm{=}\mathrm{-}\mathrm{4.4}{\mathrm{4}}_{-0.37}^{\mathrm{+}\mathrm{0.36}}$ for the primary, the secondary and the tertiary stars, respectively. The bolometric correction is computed on the basis of the expression given by Martins et al. (2005) from a typical effective temperature determined as a function of the spectral type of the three stars (see Table 4). These values coupled with the absolute V magnitude of the three stars yield the bolometric luminosities of about $\log \left({\mathit{L}}_{\mathrm{P}}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}\right)\mathrm{=}\mathrm{5.8}{\mathrm{6}}_{-0.11}^{\mathrm{+}\mathrm{0.12}}$, $\log \left({\mathit{L}}_{\mathrm{S}}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}\right)\mathrm{=}\mathrm{5.3}{\mathrm{2}}_{-0.10}^{\mathrm{+}\mathrm{0.12}}$, and $\log \left({\mathit{L}}_{\mathrm{T}}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}\right)\mathrm{=}\mathrm{5.0}{\mathrm{1}}_{-0.10}^{\mathrm{+}\mathrm{0.11}}$ for the three stars. These luminosities are in good agreement with the ones of main-sequence stars having similar spectral classifications.

The best-fit model spectra of the three components are displayed in red in Figs. 2 to 4, whilst the derived parameters are given in Table 4. However, the optical spectra do not give us enough information to constrain the wind properties of these components. We accordingly limit our investigation on the stellar parameters. Indeed, the Hα and He ii λ4686 lines are likely to be affected (to some extent) by the colliding-wind interaction zone, at least for the primary and the secondary components. Since the disentangling generates mean spectra on the entire dataset, the resulting Hα and He ii λ4686 line profiles could therefore be different from those of single stars with the same spectral classifications. We stress that, as we see in Fig. 2, both lines display P Cygni profiles for the primary, which means that the stellar wind of this star could be rather strong in comparison to its spectral classification. Therefore, we also investigate the UV domain, but the small number and the poor resolution of IUE data in the archives do not allow us to separate the individual contributions of the three components. Given the large ionizing flux generated by an O3 star, we thus assume that the flux collected in UV is mainly from the primary, whilst the He ii λ4686 and Hα lines in the optical range serve to constrain the wind parameters of the secondary and the third components. Therefore, we recommend considering the wind parameter values (Ṁ and v_∞) listed in Table 4 as preliminary ones useful for computing the model atmospheres.

7. Doppler tomography of the inner system

The O3V((f^∗)) – O3.5V((f⁺))  +  O5.5–6V((f)) inner system with its small separation between both components likely favours interactions between the stellar winds. If we look closer at the observed spectra, we clearly see that HD 150136 displays weak emissions in the He ii λ4686 and Hα lines. To better represent the emissivity of the inner system and notably the possible interactions between the wind of the primary and that of the secondary star, we use a Doppler tomography technique on those line profiles. This technique maps the formation region of these lines in velocity space. The absorption line profiles of the primary and the secondary, as well as the presence of the third component in the spectra, do not simplify the representation of the emissivity in the inner system (see left panel in Fig. 8). Therefore, a preliminary treatment to eliminate these contributions in the observed spectra must first be applied. To realize this “cleaning” process, the disentangled spectrum of the third star, shifted by its radial velocity, is subtracted from all the spectra in our dataset. We also take the He ii λ4686 line profiles of the synthetic spectra of the primary and the secondary obtained with the CMFGEN code (Sect. 6). We thus subtract them by shifting them with their respective radial velocities to remove the absorption parts of the lines. The resulting restored line profiles only take the emission part of the He ii λ4686 line into account (Fig. 8, right panel). However, such a process must be used with caution, especially for the Hα line. This line is indeed particularly dependent on the mass-loss rate and on the terminal velocity of the star. Therefore, this cleaning process requires fully reliable absorption line profiles, which is not the case for our three stars. The reconstructed emissions of the He ii λ4686 and Hα lines are nearly symmetric and follow the orbital motion of the primary component. To analyse their entire profiles, we apply a Doppler tomography technique to these two lines.

The Doppler tomography technique assumes that the radial velocity of any gas flow is stationary in the rotating frame of reference of the system. The velocity of a gas parcel, as seen by the observer, can be expressed by a so-called “S-wave” function, $\mathit{v}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}\mathrm{-}{\mathit{v}}_{\mathit{x}}\cos \mathrm{\left(}\mathrm{2}\mathit{\pi \phi }\mathrm{\right)}\mathrm{+}{\mathit{v}}_{\mathit{y}}\sin \mathrm{\left(}\mathrm{2}\mathit{\pi \phi }\mathrm{\right)}\mathrm{+}{\mathit{v}}_{\mathit{z}}\mathit{,}$where φ represents the orbital phase, (v_x,v_y) are the velocity coordinates of the gas flow and v_z is the systemic velocity of the emission line. This relation assumes an x-axis situated between the stars, from the primary to the secondary, and a y-axis pointing in the same direction as the orbital motion of the secondary star (see e.g., Rauw et al. 2002; Linder et al. 2008). Usually, the method of the Fourier-filtered back-projection based on Radon’s transform is used in astrophysics. However, this technique has the disadvantage of generating “spikes” in the tomogram if the observations do not densely sample the entire orbital cycle. To limit these artefacts, we thus apply another Doppler tomography program based on the resolution of the matricial system Af = g where A is the projection matrix, f is the vector of emissivity, and g is a vector concatenating the observed spectra. The elements of the projection matrix are equal to 1 if the projection of the S-wave onto the (v_x,v_y) plane is verified or to 0 otherwise. The solutions of the Af = g system are computed from an iterative process called SIRT (Simultaneous Iterative Reconstruction Technique, Kak & Slaney 2001). This technique starts from an initialisation of the solution vector and computes a correction at each iteration that improves the solution of the previous iteration. These corrections are given by ${\mathit{f}}^{\mathrm{(}\mathit{n}\mathrm{+}\mathrm{1}\mathrm{)}}\mathrm{=}{\mathit{f}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{+}\frac{\mathit{\mu }}{\mathrm{\left|}\mathrm{\right|}\mathit{A}\mathrm{|}{\mathrm{|}}^{\mathrm{2}}}{\mathit{A}}^{\mathit{t}}\mathrm{(}\mathit{g}\mathrm{-}\mathit{A}{\mathit{f}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{)}$where μ is a constant so that 0 < μ < 2 and  ||·||  is the matricial norm. This method converges to a least-square solution but requires more computational time than the Fourier-filtered back-projection (for more information on that technique, see Mahy 2011).

The tomograms are computed from multi-epoch spectra. Consequently, we have to correct the radial velocities to bring the spectra to the same systemic velocity on the entire campaign. The He ii λ4686 resulting Doppler map (Fig. 9) exhibits a single region of emission concentrated near the velocity of the centre of mass of the primary star. The maximum of emissivity is located at (v_x,v_y) ~ (50, −137) km s^-1. For the Hα line, we see a similar pattern as for the He ii λ4686 line. However, the analysis of the S-wave of the highest emission peak reveals a feature at (v_x,v_y) ~ (−24, −155) km s^-1. These resulting tomograms do not show any structures reminiscent of a wind-wind interaction or of a Roche lobe overflow. Although the emission peaks are not clearly centred on the centre of mass of the primary, we can assume given the size of the stars that the emissions of the He ii λ4686 and Hα lines appear to be formed in an unperturbed part of the primary wind. This result coupled with the P Cygni profiles obtained from the disentangling thus suggests that the wind of the primary star could be relatively strong for its spectral type. We notice that these Doppler maps are close to those computed by Rauw et al. (2002) for HDE 228766.

8. Discussion

8.1. The inner system

The projected rotational velocities and the radii quoted in Table 4 suggest that both components of the inner system of HD 150136 are in synchronous rotation, i.e., the angular velocity is similar for both stars. Moreover, from the minimum and evolutionary masses reported in Tables 2 and 4, we derive an inclination of about $\mathrm{47.2}{°}_{-2.5}^{\mathrm{+}\mathrm{2.9}}$ for the primary and of about $\mathrm{51.3}{°}_{-3.2}^{\mathrm{+}\mathrm{1.7}}$ for the secondary. Consequently, by assuming that our stellar parameters are close to the reality, we estimate the inclination of the inner system to be about 49° ± 5°. This value favours a non-eclipsing system and agrees with the previous inclination of about 50° given by Niemela & Gamen (2005). We then focus on the possible synchronous co-rotation of the inner system. From the orbital period and the radii, co-rotational velocities of about $\mathrm{27}{\mathrm{6}}_{-42}^{\mathrm{+}\mathrm{57}}$ km s^-1 and $\mathrm{18}{\mathrm{7}}_{-30}^{\mathrm{+}\mathrm{36}}$ km s^-1 for the primary and the secondary components, respectively, are found. Such velocities require an inclination of about $\mathrm{38}{°}_{-11}^{\mathrm{+}\mathrm{16}}$ and $\mathrm{47}{°}_{-15}^{\mathrm{+}\mathrm{37}}$ to agree with the measured V sini. These inclinations are close to the value estimated from the masses. Therefore, we may conclude that the inner system is probably in synchronous co-rotation.

By adopting a mass ratio of 1.6 and an asini = 28.9   R_⊙, we estimate the minimal radii of the Roche lobes of the primary and the secondary from the expression given by Eggleton (1983). We find (R_RLsini)_P = 12.1   R_⊙ and (R_RLsini)_S = 9.8   R_⊙. If this result confirms that the secondary does not fill its Roche lobe, the case of the primary depends on the inclination of the system. By assuming an inclination of 49° ± 5° and by taking the radius given in Table 4, this component would not fill its Roche lobe either.

8.2. The tertiary component

Previous investigations such as that of Niemela & Gamen (2005) or Benaglia et al. (2006), mentioned the existence of a third component in HD 150136. Sometimes, the close visual companion located at about 1.6″ (Mason et al. 1998) was quoted as this tertiary star. However, from the present analysis, we conclude that this close visual star cannot be this third component revealed by our analysis. First of all, the spectral classification of the close visual companion is estimated to be later than B3, whilst we clearly report here an O6.5–7V((f)) classification. Moreover, if we assume an orbital period in the range of 8 to 15 years for the outer orbit, the possible a   sini separation between the third component and the inner system is estimated between 3500 and 9100 R_⊙ i.e., between about 16 and 43 AU. With such a separation and assuming a distance of 1.3 kpc for HD 150136, we calculate a typical angular separation (acosi) between 14 and 38 mas, much smaller than 1.6″ as estimated for the close visual companion reported by Mason et al. (1998). On the basis of the spatial separation associated with a mass ratio between the inner system and the tertiary component of about 2.9, this third component should be detectable by interferometry according to the scheme given by Sana & Evans (2011). An observation made with SAM/NaCo (Nasmyth Adaptive Optics System, Paranal, Chile) in March 2011 allowed detecting a possible object at the limit of being significant (Sana et al., in prep.). This object would be located below 25 mas. Of course, this information must be taken with great care and needs to be confirmed by future observing campaigns.

Although HD 150136 is now established as a triple system, it could still have a higher multiplicity. Indeed, the possible close visual component located at about 1.6″ could be also bound to the system. In order to determine whether HD 150136 is a dynamically stable hierarchical triple system, we use the stability criterion quoted by Tokovinin (2004): ${\mathit{P}}_{\mathrm{out}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{out}}{\mathrm{)}}^{\mathrm{3}}\mathit{>}\mathrm{5}{\mathit{P}}_{\mathrm{in}}$where P_in and P_out are the orbital periods of the inner and the outer systems, respectively, and e_out is the eccentricity of the outer orbit. Even though we do not know precisely the exact values of the parameters of the external orbit, we can still report, in the worst case, that this criterion is satisfied for the configuration of HD 150136.

8.3. Evolutionary status

The stellar parameters determined with the CMFGEN atmosphere code allow us to give the positions of the three stars in the Hertzsprung-Russell (HR) diagram. As shown in Fig. 10, all these components are located relatively close to the ZAMS. The isochrones reveal an age of about 0–2 Myr for the primary and the secondary stars, whilst the tertiary would be slightly older with an age of about 1–3 Myr. Even though this discrepancy is barely significant, such a difference in age was already observed in several studies. Indeed, Martins et al. (2012, and references therein) reported that a possible reason to see this difference in the inferred age could be attributed to a different rotational rate between the components. By taking an inclination of 49° into account for the primary and the secondary components and by assuming that the tertiary is in the same orbital plane as the inner binary, we measure rotational velocities of about 227, 180, and 95 km s^-1 for the primary, the secondary, and the tertiary components, respectively. To represent them together, we use the evolutionary tracks computed with an initial rotational velocity of 300 km s^-1 (Fig. 10). However, evolutionary tracks based on the individual rotational rates are also investigated, but this does not change the conclusions: the discrepancy does not seem to come from a different rotational rate. It has to be stressed that the evolution of the rotational velocity, hence the efficiency of rotational mixing, is probably quite different in a close binary system where the stars corotate synchronously, on the one hand, and in a rather isolated third star, on the other.

This discrepancy in age is strengthened by determining the stellar parameters (see Sect. 6). Even though the parameters of the primary and the secondary stars agree with those of stars on the main sequence, we observe for the third component that its luminosity corresponds to a main-sequence star, whilst its gravity is more reminiscent of an evolved star. However, log g is determined on the basis of the wings of Balmer lines, and its estimate could be wrong if the wings of these lines are poorly reproduced by the disentangling process. In the case of HD 150136, we see that, due to the broadness of the Balmer lines, the spectral separation between the different components is too small to sample the full width of these lines. Therefore, the deduced value of the log g could be affected by the disentangling process.

Moreover, the parameters of the third star also show an evolutionary mass that is higher than the spectroscopic one, even though both values are lower than what the mass ratio between the inner system and the third star seems to indicate. We notice that the spectroscopic mass is derived from the radius and the log g of the star. Consequently, a wrong estimate of log g affects the determination of the spectroscopic mass. If we consider the evolutionary mass of the inner system (primary and secondary), we find a total mass of about 100 M_⊙. By adopting the mass ratio we have derived between the inner system and the third star, we obtain a mass of about 35 M_⊙ for the tertiary component. This value is close to the evolutionary mass, which would show that the spectroscopic log g would be wrong. This dynamical mass would also lead to a log g value of about 4.1, which is more reminiscent of a main-sequence star.

Finally, if we assume that the disentangling does not affect the true value of the log g, we have to look for the real phenomena that could explain such a value. From the different evolutionary schemes of Eggleton & Kiseleva (1996), it is however unlikely that the third star has a giant luminosity class. Indeed, according to these authors, a massive triple system such as HD 150136, i.e., with M_P > M_S > M_T, should first undergo a Roche lobe overflow in the inner binary. Since the primary star is more massive than the other ones, it should be the first to evolve, thus implying the mass transfer in the inner system. Therefore, the tertiary component could not be evolved, whilst the primary and the secondary still remain on the main-sequence band. Furthermore, the absence of modified surface chemical abundances and the low V sini derived for the third component also do not support an alternative scenario where the third component would have been captured by the inner system after having been ejected by a supernova kick or by dynamical interactions in the neighbourhood of HD 150136. All in all, it thus appears that the tertiary component and the two other stars belong to the main sequence and that the most probable cause of a wrong estimate of log g is the disentangling process.

8.4. Colliding wind and non-thermal radio emission

The presence of an O3V((f^∗))–O3.5V((f⁺)) star in the inner system seems to be mainly responsible for the variations observed in the He ii λ4686 and Hα lines. However, the proximity of the O5.5–6V((f)) secondary could imply a wind-wind interaction zone between the two stars, even though the disentangling and the Doppler tomography render that unlikely. According to the formalism of Stevens et al. (1992) and the parameters given in Table 4, we expect that the powerful wind of the primary star crashes onto the surface of the secondary star. This scenario is in line with the finding of Skinner et al. (2005). The collision of the primary wind with the secondary surface could generate an X-ray emission associated with the secondary inner surface. The short timescale X-ray variability could then result from occultation effects very much like those observed for CPD − 41° 7742 (Sana et al. 2005). Though this scenario seems to agree with the results of the Doppler tomography where we have seen that the emission is only due to the primary star, the position of the contact zone remains unknown because the wind parameters of the secondary star are uncertain. Therefore, we use the mass-loss rate of about Ṁ = 10^-6   M_⊙ yr^-1 derived by Howarth & Prinja (1996) as a function of stellar luminosity, and we assume an inclination of 49° for the system. These assumptions lead to a probable radiative primary wind, with a cooling parameter χ ≅ 0.87. On the basis of different values for the fundamental parameters, Skinner et al. (2005) derived χ = 0.2, also in agreement with a non-adiabatic regime for the X-ray emission from the inner wind-wind interaction region. On the other hand, any X-ray contribution coming from the colliding wind region located between the inner system and the third star should a priori be produced in the adiabatic regime, at least in a very large fraction of the long eccentric orbit.

The confirmation that HD 150136 includes a third star is worth discussing in the context of its non-thermal radio emission. The stellar wind material is quite opaque to radio photons, which are unlikely to escape if the absorbing column along the line-of-sight is too thick. The radius of the τ = 1 radio photosphere can be estimated on the basis of the formula given by Wright & Barlow (1975) and Leitherer et al. (1995), following the same approach as De Becker et al. (2004) in the case of the candidate binary system HD 168112. Using the stellar wind parameters given in Table 4, we calculate the radio photosphere radii, respectively for the primary, the secondary, and the tertiary components. These radii are plotted as a function of wavelength in Fig. 11. When discussing the radio emission from massive stars, one should keep in mind that stellar winds are radio emitters, but in the thermal regime. On the other hand, the wind-wind interaction in massive binaries could produce non-thermal radiation in the form of synchrotron radiation. Benaglia et al. (2006) reported on a clearly non-thermal spectral index between 3.6 and 6 cm, but the spectral index at longer wavelengths (13 and 20 cm) is more typical of a thermal spectrum. We therefore focus on the values derived at shorter wavelengths. While we are dealing with a triple system, let us separately consider two possible wind interaction regions: a first one within the short period sub-system (hereafter, SPC for short period collision) and a second one involving the third star (hereafter, LPC for long period collision). We can easily check whether any putative non-thermal radio emission from the SPC could contribute to the measured non-thermal radio flux. The curves in Fig. 11 are dominated by the primary contribution, which is easily explained by its high mass-loss rate. The radii at 3.6 and 6 cm are of the order of 260 and 370 R_⊙ for the primary. The projected stellar separation (a₁sini + a₂sini) is of about 28.9 R_⊙, which translates into a stellar separation of the order of 38 R_⊙ for an inclination of about 49° (see Sect. 8.1). The radio photospheres at 3.6 and 6 cm overwhelm thus completely the close primary + secondary system, and consequently the SPC is deeply buried in the opaque wind. This fact clearly rejects the possibility that the non-thermal radio photons come from the interaction between the primary and secondary components. Alternatively, they could come from the LPC. The radii of the radio photosphere at 6 cm constitute therefore minimum boundaries for the distance between the LPC and, respectively, the primary (solid curve in Fig. 11) and the tertiary (long-dashed curve). The sum of these two radii (i.e. about 390 R_⊙) should be considered as a lower limit for the stellar separation between these two stars. That value is indeed much lower than the stellar separation estimated in Sect. 8.2 on the basis of the probable period of the wide system. A more stringent lower limit could be derived using values at 20 cm, since the free-free absorption increases as a function of wavelength. However, the non-thermal nature of the radio photons at that wavelength is not guaranteed by the spectral index reported by Benaglia et al. (2006). We note that the radio photospheric radii estimated here intimately depend on the assumed stellar wind parameters, which should be considered with caution. In particular, the mass loss rate of the primary may be significantly underestimated here. However, our conclusion relies on conservative margins.

9. Conclusion

The analysis of high-resolution, high signal-to-noise data, collected over several epochs with FEROS, has allowed us to detect, for the first time, a third component physically bound to the inner system in the spectrum of HD 150136. In addition, we have revised the spectral classification of the three components to an O3V((f^∗))–3.5V((f⁺)) primary, an O5.5–6V((f)) secondary and an O6.5–7V((f)) tertiary. We also estimated their evolutionary masses as equal to 71, 36, and 27 M_⊙, respectively.

We refined the orbital solution of the inner binary system by fitting the radial velocities of the different components of the SB3 system together. The results have improved the accuracy of the majority of the parameters reported by Niemela & Gamen (2005). Moreover, we found that the inner system is probably in synchronous co-rotation. The agreement between the evolutionary and the Keplerian masses would suggest a probable inclination of about 49°, implying dynamical masses of about 64, 40, and 35 M_⊙.

From the systemic velocities measured on each epoch, we estimated that the orbit of the third star should have a period between 2950 and 5500 d and an eccentricity higher than 0.3. More data are necessary, however, to achieve an accurate orbital solution of the outer orbit. A possible alternative would be to observe HD 150136 by interferometry to constrain the missing orbital parameters.

The determination of the individual stellar parameters of the three stars has yielded a luminosity and a gravity corresponding to main-sequence stars for the primary and the secondary components, but these parameters are trickier for the third star. Indeed, we observed a clear discrepancy between its luminosity (comparable to main-sequence stars) and its log g (typical of giant stars). The origin of such a difference is thus unclear, but it seems that it could come from the disentangling procedure, notably because the spectral separation is not sufficient to sample the full width of the broadest lines, such as the Balmer lines, for the three components. Consequently, HD 150136 is probably composed of at least three O-type stars located on the main sequence and with inferred ages between 0–3 Myr. Finally, we showed that the non-thermal radio emission should come from the wind-wind interaction region between the third star and the inner short-period system, region that should be the site for particule acceleration in this multiple system.

Online material

---

Acknowledgments

This research was supported by the PRODEX XMM/ Integral contract (Belspo), the Fonds National de la Recherche Scientifique (FNRS) and the Communauté française de Belgique – Action de recherche concertée – ARC – Académie Wallonie-Europe. The authors thank Thomas Fauchez for a first detailed look to the data in the framework of a student work. We also thank Dr. F. Martins for his helpful comments and Pr. D. J. Hillier for making his code CMFGEN available. We are grateful to the staff of La Silla ESO Observatory for their technical support.

References

All Tables

All Figures