© ESO, 2016

1. Introduction

In recent years, it has been found that the majority of the massive stars belong to binary or higher multiplicity systems (Sana et al. 2012, and references therein). This situation has considerable implications for our understanding of the evolution of massive stars, but at the same time, it also offers enormous possibilities to observationally constrain the properties of these stars. Particularly interesting are unevolved double-line spectroscopic eclipsing binary systems (so-called SB2Es). The joint analysis of the photometric eclipses and radial velocities inferred from spectroscopy allows us to determine absolute masses and radii that are essentially model independent. Moreover, for unevolved systems, the properties of these stars should be good proxies of those of single massive stars.

Of special interest are eccentric SB2Es that show significant apsidal motion (e.g. Bulut & Demircan 2007; Schmitt et al. 2016). The latter arises from the fact that the gravitational field of stars in a close binary system can no longer be approximated as the gravitational field of a point-like mass. As a result, the orbits of the stars can only to first order be described as closed ellipses. A better description is achieved when it is considered that the argument of periastron ω undergoes a secular perturbation (e.g. Schmitt et al. 2016, and references therein). The rate of apsidal motion is directly related to the internal structure of the components of the binary (e.g. Shakura 1985; Claret & Giménez 1992, 2010), thus allowing us to obtain information on the internal mass-distribution of the stars. Measuring the rate of apsidal motion can also provide an estimate of the masses of the components of non-eclipsing eccentric close binaries (Benvenuto et al. 2002; Ferrero et al. 2013), although this method is strongly model dependent. To date, there are only very few massive stars for which detailed studies of apsidal motion have been performed. Only four systems among the 128 eccentric eclipsing binaries listed by Bulut & Demircan (2007) host at least one star with a mass of more than 20 M_⊙. In this paper, we discuss a fifth system: HD 152218.

HD 152218 is a member of the young open cluster NGC 6231 in the Sco OB1 association. With a binary fraction of at least 0.63, the O-star population of this cluster appears to be very rich in binary systems, especially short-period systems (Sana et al. 2008b). HD 152218 is a binary system consisting of an O9 IV primary and an O9.7 V secondary (Sana et al. 2008a). The system has a rather high eccentricity (e ≃ 0.28) given its orbital period of about 5.604 days (Hill et al. 1974; Stickland et al. 1997; Sana et al. 2008a). Otero & Wils (2005) reported HD 152218 to display shallow eclipses. They classified the system as an Algol-type eclipsing binary, but did not perform a detailed analysis of its light curve. Comparing their radial-velocity curve obtained from data collected between 1997 and 2004 with previously published solutions from the literature (Struve 1944; Hill et al. 1974; Stickland et al. 1997), Sana et al. (2008a) noted the probable presence of a significant apsidal motion with between 1.4 and 3.3° yr^-1, depending on the data sets that were considered. XMM-Newton data furthermore revealed a modulation in the X-ray flux of HD 152218 with orbital phase that was interpreted as a result of a wind-wind interaction (Sana et al. 2008a; Rauw & Nazé 2016).

In the present paper, we re-address the issue of apsidal motion in the HD 152218 binary system. In Sect. 2 we introduce the observational data. The radial velocities and individual spectra of the binary components are analysed in Sect. 3, notably to establish the rate of apsidal motion. Section 4 presents our analysis of the photometric light curve, whilst Sect. 5 compares the observed rate of apsidal motion to predictions from theoretical models. Finally, in Sects. 6 and 7 we discuss our results and present our conclusions.

2. Observational data

2.1. Spectroscopic data

Our spectroscopic analysis is based on 28 high-resolution échelle spectra obtained with the FEROS spectrograph at ESO between April 1999 and May 2004. Details on these observations and the data reduction are given in Sana et al. (2008a). To derive the radial velocities (RVs) and to reconstruct the individual spectra of the primary and secondary components of HD 152218, we applied our spectral disentangling code, which is based on the method outlined by González & Levato (2006). In this procedure, the individual spectra are reconstructed in an iterative way by averaging the observed spectra shifted into the frame of reference of one binary component after subtracting the current best approximation of the companion spectrum shifted to its current estimated radial velocity. Improved estimates of the RVs of the stars are obtained by cross-correlating the residual spectra, obtained after subtracting the companion spectrum, with a mask that is equal to zero in the continuum and equal to one at the rest wavelengths of spectral lines. The revised RVs are then computed from the centroid of the cross-correlation peak. This approach is especially powerful near conjunction phases, where conventional Gaussian fits usually fail to properly deblend the lines. As an initial approximation of the RVs, we used the values derived by Sana et al. (2008a). Our newly derived RVs are used in the coming sections to constrain the physical parameters of the system.

2.2. Photometric data

We have extracted V-band photometry of HD 152218 taken from the All Sky Automated Survey (ASAS-3, Pojmański & Maciejewski 2004) carried out at the Las Campanas Observatory in Chile. The ASAS-3 system consisted of two wide-field (8.8° × 8.8°) telescopes, each equipped with a 200/2.8 Minolta telephoto lens and a 2048 × 2048 pixels AP-10 CCD camera. These instruments were complemented by a narrow-field (2.2° × 2.2°) 25 cm Cassegrain telescope equipped with the same type of CCD camera.

The ASAS-3 photometric catalogue provides magnitude measurements performed with five different apertures varying in diameter from 2 to 6 pixels. Substantial differences between large- and small-aperture magnitudes indicate either contamination by close neighbours, in which case large apertures should be avoided, or saturation. The ASAS-3 data are not uniform in terms of exposure time and thus in terms of saturation limit, as the exposure time was changed in the course of the project from 180 s (saturation near V ~ 7.5) to 60 s (saturation around V ~ 6). The saturation limit also depends on the quality of the focus and the zenith distance. The ASAS-3 catalogue lists a grade for each epoch at which the star was observed. However, this grade is calculated for the frame as a whole, not for a specific object inside the field of view. Hence, a frame can have grade A, but still contain stars that saturated.

We have filtered the ASAS-3 data of HD 152218, keeping only data with a quality flag A and a mean error on the photometry of less than 0.034 mag. Since HD 152218 has a V magnitude close to the saturation limit of the ASAS instruments, we requested in addition that the dispersion between the various aperture photometric measurements be less than the mean error on the photometry and that the V-band magnitudes should be in the range between 7.55 and 7.80. This procedure resulted in a light curve with a total of 126 photometric data points collected between HJD 2 451 963 and HJD 2 454 707.

3. Spectral analysis

3.1. RV analysis

In close eccentric binaries, the non-spherical shape of the stars induces a secular variation of the longitude of periastron (ω) known as apsidal motion. This in turn leads to slow variations in the shape of the RV curve. Before we analyse the RVs to derive the orbital parameters and constrain the apsidal motion of HD 152218, we briefly consider the observational consequences of apsidal motion.

As a result of the apsidal motion, the notion of the orbital period becomes ambiguous and different definitions of the period can be adopted (see Schmitt et al. 2016, for a discussion of this problem). A practical consequence is that the Fourier period search carried out on data taken over several years, such as performed by Sana et al. (2008a), is biased. Since the apses move in the same sense as the orbital motion, we expect the Fourier technique to underestimate the value of the period between consecutive periastron passages. Likewise, apsidal motion will also affect the value of the eccentricity determined from a fit of data spanning several years. The apsidal motion smears out the peaks of the RV curve, thus leading to an apparently lower value of e. The orbital period, eccentricity, longitude of periastron, and time of periastron passage derived from a set of RVs collected over several years are therefore most likely biased by apsidal motion, whilst the amplitudes K₁ and K₂ should be free of such biases.

To establish the amplitudes of the RV curves, we hence computed an orbital solution applying the Liège Orbital Solution Package (LOSP) code (Sana et al. 2006), which is an improved version of the code originally proposed by Wolfe et al. (1967), to the set of RV values obtained in Sect. 2.1. The resulting set of orbital parameters (see Col. 3 of Table 1) was used as a first guess of the parameters to determine the rate of apsidal motion.

To establish the orbital solution accounting for apsidal motion, we considered the radial velocities of the primary available in the literature. Our analysis includes the photographic measurements presented by Struve (1944) with the heliocentric Julian dates computed in the same way as indicated in Sana et al. (2008a), the photographic measurements of Hill et al. (1974), the IUE data of Stickland et al. (1997), the CAT+CES data from Sana et al. (2008a), and our newly determined radial velocities from the FEROS spectra. To establish the RVs of the components of HD 152218, Stickland et al. (1997) first brought the IUE spectra to an almost absolute frame of reference using the interstellar lines. The RVs were obtained through cross-correlation with comparison stars. For phases where the lines of the two stars were not fully deblended, the widths of the individual cross-correlation functions were fixed to the values obtained at phases when the lines were resolved. See Stickland et al. (1997) for further details. We did not consider the RV data from García & Mermilliod (2001) because they were found to be subject to clerical errors (see discussion in Sana et al. 2008a). All RV data were assigned an error that was set to the mean error about the RV solution obtained for the corresponding epoch.

For each time of observation t, we adjusted these RV data with the following relation (1)where (2)γ₁ is the apparent systemic velocity of the primary, which was set to the value found for the best orbital solution for that epoch. K₁ is the amplitude of the primary RV curve, which we set to 160.3 km s^-1 as this was the value consistently found for all of the epochs. φ(t) stands for the true anomaly computed by solving Kepler’s equation with a specific value of the orbital period P_orb, the eccentricity e, and the time of periastron passage T₀. ω is the argument of periastron, ω₀ is the value of ω at time T₀, and is the rate of secular apsidal motion.

This therefore is a problem with five free parameters: P_orb, e, T₀, ω₀, and . We scanned the parameter space in a systematic way, starting near the orbital solution obtained from the FEROS data. Some examples of projections of the five-dimensional parameter space onto 2D planes are shown in Fig. 1.

Fig. 1 Confidence contours for the best-fit parameters of the adjustment of the full set of 81 RV data of the primary star with Eqs. (1) and (2). The blue and cyan contours correspond to 1σ and 90% confidence levels, respectively. The best-fit solution is indicated by the filled symbol.

In this way, we found the best-fit parameters and their 1σ errors, which are listed in Col. 4 of Table 1.

As an example, Fig. 2 illustrates the fit of the radial velocity data at four different epochs. The change in morphology of the RV curve over 60 yr between the observations of Struve (1944) and Sana et al. (2008a) is evident.

Fig. 2 Examples of the comparison between the RV data of the primary star and relations 1 and 2 with the best-fit parameters given in Col. (4) of Table 1. The top panels correspond to data from Struve (1944, and Hill et al. (1974, . The bottom panels show one epoch of IUE data from Stickland et al. (1997, and one epoch of FEROS data presented by Sana et al. (2008a) and reanalysed here (right).

Of course, the photographic data of Struve (1944) and, to a lesser extent, of Hill et al. (1974) are of lower quality than the more recent IUE and optical echelle spectra. We therefore also applied the same technique to the 61 RV points consisting of the Stickland et al. (1997) data, the CES RVs from Sana et al. (2008a) and our newly determined RVs for the FEROS data. The results are in good agreement with those of the full set of RV data (see Col. 5 of Table 1), although, of course, with slightly larger error bars. This test demonstrates that our determination of the rate of apsidal motion is reliable.

3.2. Spectral types and stellar parameters

We used the separated spectra of the components of HD 152218 to reassess the spectral types and to determine more precise values of the effective temperature as well as estimates of the surface gravity of the stars.

For the spectral classification, we applied the criterion proposed by Conti (1973) and Mathys (1988) based on the ratio between the equivalent widths (EWs) of the He i λ 4471 and He ii λ 4542 lines. In this way, we found that log W′ =log  [EW(He i λ 4471)/EW(He ii λ 4542)] amounts to 0.34 ± 0.01 for the primary and 0.66 ± 0.02 for the secondary. These values correspond to spectral types O9 and O9.7 (with O9.5 within the error bar) for the primary and secondary, respectively. To assess the luminosity classes, we compared our separated spectra with the digital atlas of Walborn & Fitzpatrick (1990). The strengths of the Si iv λλ 4089, 4116 and He ii λ 4686 absorptions suggest main-sequence or, especially in the case of the primary, a slightly more evolved luminosity class. We therefore classify the primary as O9 IV, and the secondary as O9.7 V. These results confirm the classification previously proposed by Stickland et al. (1997) and Sana et al. (2008a).

Fig. 3 Top row: comparison between the primary (left) and secondary (right) spectra obtained with the disentangling process (black lines) and the best-fit CMFGEN models (red lines). Bottom row: χ2 map of the CMFGEN fits showing the errors on our best-fit temperatures and gravities (shown by the red crosses) for the primary (left) and secondary (right).

Fig. 4 He iλ 5016 line profiles as derived from the separated spectra and Fourier transform of the latter along with the best-match rotational profile for the primary (left panel) and secondary (right panel).

Using the individual spectra obtained through spectral disentangling, we derived the values of the projected rotational velocities of the primary and secondary star. For this purpose, we applied the Fourier transform method (Gray 2008; Simón-Díaz & Herrero 2007) to the profiles of the He i λλ 4026, 4713, 4922, and 5016 lines of both the primary and the secondary star, and to the O iii λ 5592 profile of the primary. These lines are relatively well isolated in the spectra and are therefore expected to be almost free of blends. The He i λ 4026 line might be weakly polluted by He ii λ 4026, although given the spectral types of the stars, such a pollution should be at a very low level. The results are presented in Table 2 and illustrated in Fig. 4.

We found that the primary v sini_rot values are in good agreement for the different lines, resulting in a mean v₁ sini_rot of (144.6 ± 4.4) km s^-1, in good agreement with the values of Howarth et al. (1997) and Stickland et al. (1997). For the secondary, we found a larger dispersion among the various lines, resulting in a mean v₂ sini_rot of (141 ± 19) km s^-1, which is in reasonable agreement with the results of Howarth et al. (1997) and Stickland et al. (1997).

To derive quantitative values of the effective temperature and the surface gravity, we used a grid of synthetic spectra generated with the CMFGEN model atmosphere code (Hillier & Miller 1998). The line profiles in the synthetic spectra were broadened according to the projected rotational velocities inferred above. The relative strengths of the spectral lines in the separated spectra indicate a brightness ratio l(O9) / (l(O9) + l(O9.7)) of 0.67 in good agreement with the value (0.67 ± 0.1) inferred by Sana et al. (2008a) and the corresponding UV flux ratio (0.67) inferred from IUE spectra (Penny et al. 1993). The temperatures are essentially determined by the ability to correctly reproduce the relative strengths of He i and He ii lines, whilst gravity is mainly determined from the wings of H i Balmer lines. In many early-type binaries, especially those with small or moderate RV amplitudes, spectral disentangling has difficulties to correctly reproduce the wings of the Balmer lines, which limits the possibilities of using these wings in constraining log g (e.g. Raucq et al. 2016). However, in the present case, the orbital motion is sufficiently important that the disentangling works well for these lines.

The comparison between the separated spectra and the synthetic spectra that provide the best match is shown in Fig. 3. Since we found no indication for non-solar composition, our fits were performed assuming a solar chemical composition (Asplund et al. 2009).

Our CMFGEN analysis yields temperatures for the O9 and O9.7 stars of 33 400 and 29 900 K, respectively. In the same way, we determine log g = 3.81 and 3.89 for the primary and secondary, respectively. The lower panels of Fig. 3 illustrate the χ² maps of the CMFGEN model fits in the (T_eff, log g) plane. We followed the approximation of Martins et al. (2015) in renormalizing the χ² so that its minimum has a value of 1.0. The 1σ confidence contour would then be given by the χ² = 2.0 contour (see Martins et al. 2015). Clearly, a correlation between the two parameters exists, which makes it difficult to derive meaningful uncertainties in this case. We therefore adopted a formal error estimate of 1000 K on T_eff and 0.15 on log g.

4. Light-curve analysis

As a result of the inclination of the orbit and because of the orientation of its apses with respect to our line of sight, the light curve of HD 152218 displays only a single eclipse per orbital cycle, shortly before periastron passage when the primary star passes behind the secondary (see Fig. 5). There is only a single rather shallow eclipse. Therefore it is easy to conceive that the inversion of this light curve can be subject to degeneracies or correlations between model parameters. To limit the effect of these problems, we used the results of our spectroscopic study to constrain as many model parameters as possible. The mass-ratio, the temperatures of both stars, the orbital eccentricity, and the longitude of periastron at the epoch of the photometric observations were all frozen in our study. The only free parameters were the Roche-lobe filling factors of both stars fill(O9 IV) and fill(O9.7 V) and the orbital inclination (i). The filling factors are defined here as the fraction of the polar radius of the star under consideration over the polar radius of its instantaneous Roche lobe at periastron passage.

We then ran the Nightfall code (version 1.86)¹ to search for the combination of model parameters that yields the best fit of the data. Nightfall handles eccentric eclipsing binaries assuming that the shape of the stars can be described by a Roche potential scaled with the instantaneous separation between the stars. To estimate the errors on the model parameters, we performed a systematic study of the 3D parameter space. For this purpose, we mapped the contours of χ² over three grids that sample the projections of the 3D parameter space onto the planes formed by the combinations of two out of the three parameters fill(O9 IV), fill(O9.7 V), i. The resulting contours for Δχ² = 2.30 and 4.61 (1σ and 90% confidence intervals for two degrees of freedom) are displayed in Fig. 6.

However, not all combinations of parameters that are acceptable as far as the fit of the light curve is concerned agree with the spectroscopic results. From spectroscopy, we have inferred an optical brightness ratio that leads to a ratio of bolometric luminosities of L_bol(O9 IV)/L_bol(O9.7 V) = 2.17 ± 0.75. With the temperatures of the stars, this ratio translates into Accounting for the ratio of the radii of the primary and secondary Roche lobes, we then inferred fill(O9 IV)/fill(O9.7 V) = 1.04 ± 0.18. The curves corresponding to this condition are also shown in Fig. 6. Their intersections with the χ² contours for 1σ confidence intervals were used to establish the error bars on the parameters (see Table 3). Whilst the orbital inclination is rather well constrained, the relative uncertainties on the filling factors and hence on the stellar radii are quite large. The latter situation mainly stems from a slight anti-correlation between the two parameters that leads to elongated and inclined confidence contours in the fill(O9.7 V) versus fill(O9 IV) plane (see Fig. 6).

Fig. 5 Photometric light curve of HD 152218 as derived from the ASAS-3 data along with the best-fit model computed with the Nightfall code.

From these considerations we can then estimate absolute radii and masses of the components of HD 152218 as given in Table 3. In this table we also compare the values of log g obtained from our CMFGEN analysis with those that we obtain using the absolute masses and radii as determined from the analysis of the light curve. Within the error bars, both estimates agree.

Assuming that the rotational axes of the stars are aligned with the orbital axis, we can use our best-fit value of the inclination and stellar radii and combine them with the projected rotational velocities from Table 2. In this way, we derive rotational periods of  days for the primary and  days for the secondary. Within the error bars, both stars thus appear to have their rotations nearly synchronized with each other. The rotational periods we derived are shorter by about a factor two than the orbital period. When we compare the rotational angular velocities with the instantaneous orbital angular velocity at periastron, we find ratios of and for the primary and secondary, respectively. We therefore conclude that pseudo-synchronization (Hut 1981) has been achieved in the system, or is on the verge of being achieved. These properties make HD 152218 an interesting target to search for tidally induced oscillations (Moreno et al. 2011; Palate et al. 2013, and references therein).

Fig. 6 Contours of 1σ (blue) and 90% confidence intervals (cyan) in the various parameter spaces of the light curve analysis. The red curve corresponds to a ratio of the filling factors of 1.04 (see text), whilst the magenta curves yield the corresponding 1σ interval, as derived from spectroscopy. The filled dot indicates the set of parameters that yields the best fit of the ASAS light curve.

5. Apsidal motion

The observationally derived rate of apsidal motion can be compared with the rates predicted from theoretical models. The rate of apsidal motion in a close binary system can be directly related to the internal structure constants k₂ of both stars (e.g. Shakura 1985). The internal structure constants are obtained from (3)where η₂(R_∗) is the solution of the Clairaut-Radau differential equation (4)evaluated at the stellar surface. In Eq. (4), r is the radial distance from the centre of the star, ρ(r) is the local density, and ⟨ ρ ⟩ (r) is the mean density inside a sphere of radius r (e.g. Sterne 1939; Schmitt et al. 2016). Solving the Clairaut-Radau equation applied to the density structure predicted by stellar evolution models thus yields the value of k₂ for a given stellar mass and age.

To this aim, we have designed a FORTRAN code using a fourth-order Runge-Kutta method with step doubling (Press et al. 1997). To validate our code, we adopted the same method as İnlek & Budding (2012), considering a series of polytropic models. For this purpose, we solved the Lane-Emden differential equation for different values of the polytropic index n and then applied our Clairaut-Radau integration routine to the resulting density profiles. We then compared our values of k₂ to those listed by Brooker & Olle (1955) and found a very good agreement with relative differences of less than 10^-5.

After the validation of our code, we applied it to a series of stellar structure models computed with the Code Liégeois d’Évolution Stellaire (CLES, Scuflaire et al. 2008) for stars of 15 and 20 M_⊙. Internal mixing was either restricted to the convective core or included a core overshooting of 0.2 times the pressure height scale (α_ov = 0.2). The chemical composition of the stars was assumed to be solar adopting the abundances of Asplund et al. (2009). Mass-loss was included according to the prescription of Vink et al. (2001). We tested both a solar metallicity (Z = 0.014) and a slightly increased metallicity of Z = 0.020.

Table 4 lists the model stellar radii and corresponding k₂ values as a function of stellar age. The decrease in k₂ values as the stars evolve is evident. At a given age, overshooting has a moderate effect on both the stellar radius, which decreases compared to the models without overshooting, and the value of k₂, which decreases for the 20 M_⊙ models, but increases for the 15 M_⊙ models. The effect of an enhanced metallicity is to increase the stellar radii and simultaneously decrease the values of k₂.

Fig. 7 Comparison between the observed value of (solid horizontal line) with its errors (dashed horizontal lines) and the values evaluated from stellar structure models as a function of stellar age. The solid curve yields the results obtained by interpolating the models of Claret & Giménez (1992). The filled squares correspond to the results obtained through the resolution of the Clairaut-Radau equation applied to the CLES models for a 20 M⊙ primary and a 15 M⊙ secondary. The open squares yield the results for CLES models assuming an overshooting with αov = 0.2. The filled diamonds stand for the CLES models with a slightly higher metallicity.

We then computed the theoretical rates of apsidal motion that are due to the linear part of the quadrupole moment for the binary components and the relativistic contribution using the following equation adapted from Shakura (1985): (5)where (6)and (7)stands for the rate of apsidal motion due to general relativity (Giménez 1985). Indices 1 and 2 refer to the primary and secondary star, respectively, and . As we are interested in determining how the theoretical rate of apsidal motion evolves with age, the stellar radii used in Eq. (5) were taken from the evolutionary models rather than from the photometric solution. The rotational periods were computed from the projected equatorial rotational velocities (see Table 2), assuming the rotation axes to be aligned with the normal to the orbital plane and adopting again the radii of the stars as predicted by the evolutionary models. As noted in Sect. 4, HD 152218 is currently most likely in a pseudo-synchronization state. As an alternative to our above approach to the calculation of the rotational periods, we might therefore assume that both stars have reached pseudo-synchronization early on. Using this assumption instead of the values of v sini_rot does not change the predicted values of in a significant way, at least not as far as the range of ages is concerned where the predicted matches the observed value.

The stellar radii play a key role in determining the value of because the f(e) term of the rate of apsidal motion depends on (R_∗/a)⁵, whilst the g(e) term still depends on (R/a)³ (see Eq. (5)). Therefore, the decrease of the stellar radii for the models with overshooting more than compensates for the changes of k₂, resulting in lower predicted values of . In the same way, the larger radii of the higher Z models, compared to models with Z = 0.014 of the same age, lead to higher theoretical values of .

Figure 7 illustrates the comparison between our observationally determined value of and the theoretical predictions based on the CLES models (Table 4) for primary and secondary masses of 20 and 15 M_⊙ and adopting the eccentricity and orbital period from Table 1. The solid curve in Fig. 7 furthermore illustrates the results of an interpolation of the Claret & Giménez (1992) models (again for primary and secondary masses of 20 and 15 M_⊙), which assume α_ov = 0.2 and Z = 0.020. In each case, we include a relativistic contribution of  yr^-1 evaluated according to Shakura (1985) and Giménez (1985).

The comparison between the model predictions and the observations yields an estimated age of the HD 152218 binary system in the range between 5.2 and 6.4 Myr, accounting for the uncertainties on and the uncertainties due to the different prescriptions for α_ov, and assuming Z = 0.014. When instead we were to assume Z = 0.020, our age estimate would be reduced by 0.3 Myr.

Finally, we have also tested models with an enhanced helium abundance corresponding to X = 0.66 instead of X = 0.70 adopted above. For a given age, these models have larger radii that lead to higher predicted values of and thus an age between 4.5 and 5.5 Myr. We caution, however, that the analysis of seven B-type stars by Mathys et al. (2002) resulted in log ϵ_He = 10.81 ± 0.20, whilst an earlier study of five B-stars by Keenan et al. (1984) yielded log ϵ_He = 10.92 ± 0.11. Both results are compatible with the solar He abundance (log ϵ_He = 10.93 ± 0.01, Asplund et al. 2009). Furthermore, our separated spectra of HD 152218 do not provide evidence of an enhanced helium abundance. This means that observational data support our X = 0.70 results.

Assuming the rotational axes to be aligned with the direction perpendicular to the plane of the orbit, Benvenuto et al. (2002) have shown that the rate of apsidal motion (Eq. (5)) can be expressed as a function of the primary mass as well as of several parameters that can be determined observationally (mass ratio, projected rotational velocities, eccentricity, etc.) or theoretically (internal structure constants and stellar radii). A simplified formalism based on the same idea was proposed by Jeffery (1984). We note that

and (10)which lead to (11)Benvenuto et al. (2002) used this formulation of apsidal motion to constrain the masses of the components of the non-eclipsing O3 V + O8 V binary HD 93205 by comparing the observed rate of apsidal motion against the predictions from a grid of theoretical models (see also Ferrero et al. 2013; Pablo et al. 2015, for other applications of this method). Whilst this method for constraining the masses is by construction highly model dependent and sensitive to the assumed age, it is interesting to apply it to HD 152218.

For this purpose, we have computed a grid of CLES evolutionary models with α = 0.0 and Z = 0.014 for masses between 11 and 24 M_⊙. Using our newly determined binary parameters (q, e, P_orb, a sini) and projected rotational velocities, we then determined the age of the system that allows us to reproduce the observed rate of apsidal motion for a given value of m₁. The results are shown in Fig. 8.

Fig. 8 Relation between the age of HD 152218 and the mass of the primary star. The blue solid curve corresponds to the combination of age and primary mass in our CLES models (assuming Z = 0.014 and α = 0.0) that reproduce exactly the observed rate of apsidal motion through Eq. (11). The blue dashed curves correspond to the 1σ errors on the observed . The red vertical lines correspond to the absolute mass of the primary determined from our orbital and photometric solution (Table 3) and its errors.

Using the absolute mass of the primary and its uncertainty (see Table 3), we can then estimate an age of  Myr, broadly confirming the results obtained above.

6. Discussion

The age of the early-type stellar population of the open cluster NGC 6231 has been estimated using various techniques, and an overview of the relevant literature can be found in Reipurth (2008). Age determinations of young open clusters mostly rely upon isochrone fitting to the population of early-type stars. They therefore obviously depend on the ingredients of the evolutionary models (e.g. opacities, metallicity, inclusion or not of rotation, etc.) and the accuracy of calibrations inferred from stellar atmosphere models (e.g. bolometric corrections, effective temperatures, etc.). These methods are also biased by the inclusion of binary systems that are not recognized as such. Finally, as revealed by studies of the populations of pre-main sequence stars, low-mass star formation in young open clusters proceeds over at least about 10 Myr, and it is rather unlikely that all massive cluster members formed simultaneously. For NGC 6231, Raboud et al. (1997) and Sung et al. (1998) reported on an age spread among the low-mass pre-main sequence stars of at least 7 Myr (Raboud et al. 1997) or 12 Myr (Sung et al. 1998), which was subsequently confirmed by Sana et al. (2007). We therefore expect a spread of ages and not a single value. With these caveats in mind, we briefly consider the existing age determinations of the massive stars in NGC 6231.

van Genderen et al. (1984) compared VBLUW photometry of 9 O-type members of the Sco OB1 association with isochrones and inferred an average age of 3.6 ± 0.6 Myr. However, only four stars in their sample are members of NGC 6231. Moreover, three of them (HD 152218, HD 152247, and HD 152248) are binaries, whilst the fourth star (HDE 326331) displays a variable radial velocity (Sana et al. 2008b).

Perry et al. (1991) compiled photometric data from the literature. Most of their stars had masses between 10 and 25 M_⊙ and were consistent with log  (age) = 6.9 ± 0.2 (i.e. an age of  Myr), although they noted that a few more massive stars were apparently younger.

Raboud et al. (1997) studied Geneva photometry of the NGC 6231 cluster and fitted isochrones to the upper main-sequence, inferring an age of 3.8 ± 0.6 Myr. These authors explicitly excluded known binary stars from their analysis. Still, the census of massive binaries in the cluster was not complete at that time, and it might well be that their sample still included some binaries.

Sung et al. (1998) performed isochrone fitting to UBVRIHα photometry and inferred an age of between 2.5 and 4.0 Myr for the massive star population.

Finally, Baume et al. (1999) studied UBVI photometry. These authors emphasized the difficulties that the large number of binaries introduces in the age determination. They estimated an age of between 3 and 5 Myr.

In summary, whilst there is no consensus on the age of NGC 6231, most determinations, with the exception of the study of Perry et al. (1991), favour a value near 3.8 ± 0.5 Myr. We stress that these literature results assumed Z = 0.020, although Baume et al. (1999) drew attention to the work of Kilian et al. (1994), who analysed the chemical composition for a sample of ten B-type stars in NGC 6231 and found their CNO mass fraction to be subsolar (between 0.7 and 1.3% compared to the then-accepted value of 1.8% for the Sun).

As our newly determined value of the age of HD 152218 is significantly higher than the estimates above, we need to ask ourselves whether our value could be biased in some way.

Since we used single-star evolution models in our analysis, our results rely upon the assumption that HD 152218 has not yet experienced any mass and angular momentum exchange through Roche-lobe overflow (RLOF). If the system were in a post-RLOF evolutionary stage, the exchange of mass between the stars could have led to an apparent rejuvenation of the mass gainer (Dray & Tout 2007) and would in any case affect our determination of the age. However, given the properties of the system, the assumption that the system has not yet experienced any RLOF interaction seems quite reasonable. In fact, contrary to well-established post-RLOF systems (e.g. Linder et al. 2008; Mahy et al. 2011; Raucq et al. 2016), the spectra of the components of HD 152218 do not show evidence of enhanced nitrogen abundance or asynchronous rotation. Furthermore, the orbital eccentricity is quite high for a short-period binary, and the Roche-lobe filling factors at periastron inferred from the photometric light curve (see Sect. 4) are well below unity. Therefore, it seems very unlikely that HD 152218 has experienced binary interaction in the past.

In Sect. 5 we investigated the effect of overshooting on the values of k₂ and . Whilst the assumptions on internal mixing clearly do have an effect on the age determination, this is rather moderate, as illustrated by the relatively good agreement between different models in Fig. 7. An additional effect that might influence the predictions of the theoretical models, but is not accounted for in our models, is rotational mixing. Rotational mixing in massive stars changes their position in the Hertzsprung-Russell diagram (e.g. Meynet & Maeder 2000). However, the rotational velocities that we have derived (see Table 2) are rather moderate for massive stars. Moreover, by comparing evolutionary tracks for a non-rotating 20 M_⊙ star and a star of the same mass rotating at ~132 km s^-1 (see Fig. 7 of Meynet & Maeder 2000), we note that there is almost no difference in luminosity (hence in radii) between the models when the stars have effective temperatures of about 33 400 K (log T_eff ≃ 4.52). It therefore appears unlikely that the neglect of rotation in our models could lead to a significant bias on the age determination. Claret (1999) investigated the effect of rotation on the internal structure constants. He found that for stars that are not too heavily distorted, rotation induces a reduction of the theoretical k₂ given by (12)where and refer to the internal structure constants with and without rotation, and Ω_s is the angular rotational velocity at the stellar surface. Applying this recipe to our case leads to a reduction of the internal structure constants by ~8 and ~7% for the primary and secondary star, respectively. Neglecting any effect of rotation on the stellar radii (see above), we then expect a reduction of the theoretical rate of apsidal motion by the same ~8%. Therefore, accounting for rotation leads to even higher values of the ages.

Finally, we made the assumption that the rotational axes are aligned with the normal to the orbital plane. Shakura (1985) presented a formalism (his Eq. (3)) that allows us to deal with situations where this assumption is not fulfilled. Assuming that the rotation periods of the two binary components are known, the angle-dependent part of the rotation term of takes the form (Shakura 1985) (13)where i denotes the inclination of the orbital plane with respect to the observer, i_rot is the angle of the rotation axis with respect to the observer, and γ is the misalignment angle between the normal to the orbit and the rotation axis. For alignment we have γ = 0 and i_rot = i, thus F(i,i,0) = 1. We adopt a Cartesian coordinate system with the z-axis pointing at the observer and the y-axis chosen such that the vector normal to the orbital plane is located in the (y,z) plane. To uniquely specify the rotation axis in space we need (in addition to the already introduced angle γ) an azimuthal angle θ that measures the precession angle of the stellar spin axis around the normal to the orbital plane. With the convention that θ = 0 yields a spin vector in the (y,z) plane, this yields the following relation between the angles θ, γ, i, and i_rot: (14)Substituting Eqs. (14) into (13), we obtain (15)from which we deduce (16)and (17)as the maximum and minimum (relative) contributions for the rotational terms to the apsidal motion rate.

Fig. 9 Maximum (solid curve) and minimum (dashed curve) values of the misalignment correction factor F (see Eq. (13)) as a function of misalignment angle γ.

In Fig. 9 we show these curves for the specific case of HD 152218 (with i = 66.3°); the area between the two curves shows the possible range of values of F (as in Eq. (13)) as a function of the misalignment angle. The highest possible rotational contribution is reached at a non-zero misalignment angle γ_max, which can be trivially computed from (18)However, this maximum contribution exceeds the value under perfect alignment only by 3% and for most parameters the rotational contribution to the apsidal motion rate is diminished compared to the situation of perfect alignment. Hence, if we were to assume a misalignment between the rotational and orbital axes, we would reduce the predicted as a function of age. As a result, such a misalignment would increase the age estimate based on the theoretical rate of apsidal motion.

Fig. 10 Position of the O-type stars of NGC 6231 in a Hertzsprung-Russell diagram, accounting for the binarity. The components of HD 152218 are shown by the filled magenta symbols. Evolutionary tracks (identified with their initial mass) for Z = 0.014 are taken from Ekström et al. (2012). The red lines yield isochrones for log age = 6.5 and 6.6 (dashed lines, i.e. ages of 3.2 and 4.0 Myr, from left to right), and log age = 6.71,6.76 and 6.81 (continuous lines, i.e. ages of 5.1, 5.8 and 6.5 Myr, from left to right). The latter isochrones bracket the range of ages inferred from the apsidal motion of HD 152218. The leftmost panel corresponds to evolutionary models without rotation, whilst the rightmost panel yields the tracks and isochrones for stars rotating at 40% of their critical velocity.

To conclude this discussion about a possible misalignment of the rotational axes, we briefly consider the precession of the rotation axes that would result for such a situation (Kopal 1968, 1969a,b; Alexander 1976). From an observational point of view, precession would manifest itself through slow variations of the v sini_rot values. Using the formalism of Kopal (1969b, his Eq. (5.46)) along with our theoretical values of k₂ and the moments of inertia derived from the CLES models, we estimate precession periods in the range 14−19 yr for the range of possible ages inferred in Sect. 5. The fact that the IUE and FEROS spectra, taken about eight years apart, indicate consistent values of v sini_rot and the lack of conspicuous variations of the line widths over the five years covered by the FEROS data suggest that any possible misalignment of the rotation axes probably is rather modest.

We therefore cannot identify any obvious bias in our result. Support for an age of about 6 Myr comes from the photometric light curve. Although the agreement between the stellar radii inferred from the photometric light curve and the theoretical radii is certainly not perfect, the stellar radii inferred from the solution of the light curve clearly favour an age significantly in excess of 3 Myr.

Finally, to check whether such an age is consistent with the general properties of the population of massive stars in NGC 6231, we reconsidered all the O-type cluster members (except for the companion of the Wolf-Rayet binary WR 79) listed by Sana et al. (2008b). For the known binary systems, we corrected the magnitudes for the contribution of the companion. Whenever a dedicated spectroscopic study was available, we adopted the brightness ratio established in this study. For all other cases, we estimated the brightness ratio based on the preliminary spectral types listed by Sana et al. (2008b). We also adopted a distance modulus of 11.07 ± 0.04 and R_V = 3.3 (Sana et al. 2008a). Effective temperatures and bolometric corrections were taken from the calibrations of Martins et al. (2005) and Martins & Plez (2006). Except for cases where the uncertainties on spectral types are large (weakly detected secondary components of long-period binaries), we assumed an uncertainty of half a spectral subtype both in temperature and bolometric correction. The location of the stars in a Hertzsprung-Russell diagram is shown in Fig. 10 along with single-star evolutionary tracks of Ekström et al. (2012) for Z = 0.014, both without rotation (left panel) and with a 40% critical rotation. This comparison shows that the O-type stars seem to span a range of ages, roughly between 3 and 6.3 Myr. Our values of the age of HD 152218 as determined in Sect. 5 are towards the higher end of this range, but are consistent with the range of ages inferred with the isochrones.

7. Conclusion

We have re-analysed the massive binary HD 152218 in the young open cluster NGC 6231. Spectral disentangling allowed us to reconstruct individual spectra and hence to derive stellar temperatures and gravities comparing the individual spectra with CMFGEN model atmospheres. Photometric data were used to constrain the orbital inclination and the Roche-lobe filling factors. Radial velocity data from the literature allowed us to establish a rate of apsidal motion of  yr^-1 , corresponding to a period of 176 yr. Comparison of this rate with the predictions of stellar structure models yields an age of 5.8 ± 0.6 Myr for HD 152218. This value is towards the higher end of the range of ages of the isochrones that encompass the location of the O-star population of NGC 6231 in a Hertzsprung-Russell diagram.

Acknowledgments

The Liège team acknowledges support through an ARC grant for Concerted Research Actions, financed by the Federation Wallonia-Brussels, from the Fonds de la Recherche Scientifique (FRS/FNRS), including an FRS/FNRS Research Project (T.0100.15), as well as through an XMM PRODEX contract (Belspo).

References

All Tables

All Figures

	Fig. 1 Confidence contours for the best-fit parameters of the adjustment of the full set of 81 RV data of the primary star with Eqs. (1) and (2). The blue and cyan contours correspond to 1σ and 90% confidence levels, respectively. The best-fit solution is indicated by the filled symbol.
In the text

	Fig. 2 Examples of the comparison between the RV data of the primary star and relations 1 and 2 with the best-fit parameters given in Col. (4) of Table 1. The top panels correspond to data from Struve (1944, and Hill et al. (1974, . The bottom panels show one epoch of IUE data from Stickland et al. (1997, and one epoch of FEROS data presented by Sana et al. (2008a) and reanalysed here (right).
In the text

	Fig. 3 Top row: comparison between the primary (left) and secondary (right) spectra obtained with the disentangling process (black lines) and the best-fit CMFGEN models (red lines). Bottom row: χ2 map of the CMFGEN fits showing the errors on our best-fit temperatures and gravities (shown by the red crosses) for the primary (left) and secondary (right).
In the text

	Fig. 4 He iλ 5016 line profiles as derived from the separated spectra and Fourier transform of the latter along with the best-match rotational profile for the primary (left panel) and secondary (right panel).
In the text

	Fig. 5 Photometric light curve of HD 152218 as derived from the ASAS-3 data along with the best-fit model computed with the Nightfall code.
In the text

	Fig. 6 Contours of 1σ (blue) and 90% confidence intervals (cyan) in the various parameter spaces of the light curve analysis. The red curve corresponds to a ratio of the filling factors of 1.04 (see text), whilst the magenta curves yield the corresponding 1σ interval, as derived from spectroscopy. The filled dot indicates the set of parameters that yields the best fit of the ASAS light curve.
In the text

Fig. 7 Comparison between the observed value of (solid horizontal line) with its errors (dashed horizontal lines) and the values evaluated from stellar structure models as a function of stellar age. The solid curve yields the results obtained by interpolating the models of Claret & Giménez (1992). The filled squares correspond to the results obtained through the resolution of the Clairaut-Radau equation applied to the CLES models for a 20 M⊙ primary and a 15 M⊙ secondary. The open squares yield the results for CLES models assuming an overshooting with αov = 0.2. The filled diamonds stand for the CLES models with a slightly higher metallicity.

In the text

Fig. 8 Relation between the age of HD 152218 and the mass of the primary star. The blue solid curve corresponds to the combination of age and primary mass in our CLES models (assuming Z = 0.014 and α = 0.0) that reproduce exactly the observed rate of apsidal motion through Eq. (11). The blue dashed curves correspond to the 1σ errors on the observed . The red vertical lines correspond to the absolute mass of the primary determined from our orbital and photometric solution (Table 3) and its errors.

In the text

	Fig. 9 Maximum (solid curve) and minimum (dashed curve) values of the misalignment correction factor F (see Eq. (13)) as a function of misalignment angle γ.
In the text

Fig. 10 Position of the O-type stars of NGC 6231 in a Hertzsprung-Russell diagram, accounting for the binarity. The components of HD 152218 are shown by the filled magenta symbols. Evolutionary tracks (identified with their initial mass) for Z = 0.014 are taken from Ekström et al. (2012). The red lines yield isochrones for log age = 6.5 and 6.6 (dashed lines, i.e. ages of 3.2 and 4.0 Myr, from left to right), and log age = 6.71,6.76 and 6.81 (continuous lines, i.e. ages of 5.1, 5.8 and 6.5 Myr, from left to right). The latter isochrones bracket the range of ages inferred from the apsidal motion of HD 152218. The leftmost panel corresponds to evolutionary models without rotation, whilst the rightmost panel yields the tracks and isochrones for stars rotating at 40% of their critical velocity.

In the text