© ESO 2022

1. Introduction

The young open cluster NGC 6231, at the core of the Sco OB1 association (Reipurth et al. 2008; Kuhn et al. 2017), hosts a substantial population of O-type stars, a large fraction of which were shown to be binary systems (Sana et al. 2008, and references therein). Based on the properties of the population of low-mass stars, it has been shown that the cluster has an age spread of 1–7 Myr with a slight peak near 3 Myr (Sana et al. 2007; Sung et al. 2013). Variations of the mean age across the cluster are, however, small and no obvious global trend (e.g., a radial age gradient) was found (Kuhn et al. 2017). The cluster core has a radius of 1.2 ± 0.1 pc and the cluster has probably considerably expanded from its initial size (Kuhn et al. 2017). Stars more massive than 8 M_⊙ appear more concentrated in the cluster core, whereas low and intermediate mass stars display no mass segregation (Raboud & Mermilliod 1998; Kuhn et al. 2017).

Since the rate of apsidal motion in eccentric short-period binaries depends on the evolutionary stage, and thus the age of the binary components, it can help us constrain the star formation history of the cluster. So far, this technique has been applied to two massive binaries of NGC 6231: Rauw et al. (2016) determined an age of 5.8 ± 0.6 Myr for HD 152218, while Rosu et al. (2020b) determined an age of 5.15 ± 0.13 Myr for HD 152248. Here, we report our analysis of a third system, HD 152219.

The spectroscopic binarity of HD 152219 was first reported by Hill et al. (1974) who presented an SB1 orbital solution with an orbital period of 4.16 days and an eccentricity near 0.1. Subsequent SB1 orbital solutions were published by Levato & Morrell (1983) and García & Mermilliod (2001). The systematic detection of the spectral signature of the secondary star was achieved by Sana et al. (2006) who proposed an O9.5 III + B1-2 V-III spectral classification and derived the first full SB2 orbital solution. Considering their own data combined with radial velocities (RVs) from the literature, Sana et al. (2006) obtained a dataset consisting of 79 observations spanning 13 186 days and derived an orbital period of 4.24028 days. The photometric eclipses of HD 152219 were discovered by Otero & Wils (2005) based on ASAS-3 data. These authors derived a photometric period of 4.24038 days.

Sana et al. (2006) reported apparent line profile variations of the primary star which they tentatively attributed to non-radial pulsations. However, in a subsequent study, Sana (2009) used intensive spectroscopic monitoring of the star to show that these line profile variations are restricted to a short phase interval around the primary eclipse and arise from the Rossiter-McLaughlin effect. Any additional line profile variability due to non-radial pulsations would have an amplitude below 0.5% of the continuum level (Sana 2009).

Over recent years, an increasing number of massive stars have been found to be part of multiple systems (Duchêne & Kraus 2013; Sana et al. 2012). In the case of HD 152219, some of the more distant neighbours might actually be cluster members without a gravitational connection to the binary system. Sana et al. (2014) detected six companions with the NACO instrument. The closest component was found at an angular separation of 83.6 mas, with a ΔK_s of 2.8 mag. Five additional companions were found at angular separations between 2.2 and 7.2 arcsec and with K_s magnitude differences (compared to the main binary) between 4.2 and 6.9 mag. The closest companion was however not detected during a speckle interferometric survey by Mason et al. (2009).

The present paper is organised as follows. The observational data are introduced in Sect. 2. In Sect. 3 we perform spectral disentangling, reassess the spectral classification of the stars, and analyse the reconstructed spectra by means of the CMFGEN model atmosphere code (Hillier & Miller 1998). The radial velocities (RVs) inferred from the spectral disentangling are combined with data from the literature in Sect. 4 to establish a value for the apsidal motion rate of the system. Photometric data are analysed in Sect. 5 by means of the Nightfall binary star code. The link between apsidal motion and the internal structure of the stars is recalled in Sect. 6 and a weighted-average mean of the internal structure constant is inferred. Stellar structure and evolution models are computed with the Clés code in Sect. 7 where theoretical results are confronted to observational constraints. Finally, in Sect. 8, we investigate the effects that could bias our interpretation of the apsidal motion in terms of the internal structure constant, notably the impact of a possible misalignment of the stellar rotation axis with respect to the normal to the orbital plane, and we give our conclusions in Sect. 9.

2. Observational data

2.1. Spectroscopy

To investigate the optical spectrum of HD 152219, we extracted 93 high-resolution échelle spectra, retrieved from the ESO archive. All were obtained with the Fiber-fed Extended Range Optical Spectrograph (FEROS) mounted on the European Southern Observatory (ESO) 1.5 m or 2.2 m telescopes in La Silla, Chile (Kaufer et al. 1999). These data were collected between May 1999 and June 2006 (Sana et al. 2006; Sana 2009). The FEROS instrument has a spectral resolving power of 48 000, and its detector is an EEV CCD with 2048 × 4096 pixels of 15 × 15 μm. Thirty-seven orders cover the wavelength domain from 3650 to 9200 Å. Exposure times range from 300 to 1200 s. The data were reduced using the FEROS pipeline of the MIDAS software. Residual cosmic rays were removed within MIDAS, and the telluric tool within IRAF was used along with the atlas of telluric lines of Hinkle et al. (2000) to remove the telluric absorptions near the Hα and He Iλ 5876 lines. The spectra were normalised with MIDAS by fitting low-order polynomials to the continuum. The journal of the spectroscopic observations is presented in Appendix A, Table A.1.

2.2. Photometry

Photometric data of HD 152219 exist in various public databases. Unfortunately, not all of them were useful for our analysis. Indeed the white light observations collected with the Optical Monitor Camera (Mas-Hesse et al. 2003) onboard the INTEGRAL spacecraft¹ are too scarce to perform a meaningful analysis of the light curve. Likewise, the HIPPARCOS/Tycho data (ESA 1997) have large uncertainties and HD 152219 (m_V ∼ 7.57) is heavily saturated in the All-Sky Automated Survey for Supernovae (ASAS-SN, Kochanek et al. 2017)² data. The remaining useful data are the V-band data from the All Sky Automated Survey (ASAS-3, Pojmański & Maciejewski 2004)³ that allowed Otero & Wils (2005) to discover the existence of photometric eclipses and the observations collected by the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015).

The ASAS-3 data were obtained at the Las Campanas Observatory with 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, complemented by a 25 cm Cassegrain telescope equipped with the same type of CCD camera but covering a narrower field (2.2° ×2.2°). The ASAS-3 archive provides aperture photometry for five different apertures varying in diameter from two to six pixels. We adopted the four pixel aperture data in our analysis. Following Mayer et al. (2008), we discarded the observations taken before HJD 2 452 455. To further clean the data from obvious outliers, we applied a 1.25σ clipping to the data about the mean light curve built from normal data points computed for phase bins of 0.01. With a few exceptions, this criterion allowed us to remove all points that were off by ≥0.05 mag from the mean magnitude in the phase bins. A few remaining discrepant points were associated with phase bins containing very few data points, but were not removed afterwards. The dispersion about the mean within 0.01 phase bins (average over all phase bins containing more than one data point) evolved from an initial value of 4.7 10⁻² mag to 3.0 10⁻² mag after the 1.25σ clipping. In this way, we ended up with a total of 469 data points obtained between HJD 2 452 456 and HJD 2 455 112. The phase-folded ASAS-3 light curve is illustrated in Fig. 1 adopting the orbital period given in Table 3.

Fig. 1. ASAS-3 light curve phase-folded using the orbital period given in Table 3. The top panel illustrates the initial ASAS-3 data set, whilst the bottom panel shows the data after the 1.25σ clipping.

TESS is an all-sky survey space telescope operated by NASA that provides high-precision photometry for a huge number of stars. TESS is equipped with four 10 cm aperture refractive cameras covering the sky in sectors of 24° ×96°. Each sector is observed during two consecutive spacecraft orbits. In the case of HD 152219, the data were taken with camera one during sector 12 (i.e. between 21 May and 19 June 2019, hereafter TESS-12) and sector 39 (i.e. between 27 May and 24 June 2021, hereafter TESS-39). The TESS bandpass ranges from 6000 Å to 1 μm (Ricker et al. 2015). The CCD detectors have 15 × 15 μm² pixels which correspond to (21″)² on the sky, and undersample the instrument PSF. Photometry is always extracted over several pixels. Bright objects saturate the central pixel, but their photometry can be recovered in the pipeline processing (Jenkins et al. 2016) from the excess charges that are spread into adjacent pixels via the blooming effect. For TESS-12, we extracted the 2 minutes high-cadence light curve of HD 152219 from the Mikulski Archive for Space Telescopes (MAST) portal⁴. These light curves provide background-corrected aperture photometry as well as so-called PDC photometry obtained after removal of trends that correlate with spacecraft or instrumental effects. We discarded all data points with a quality flag different from zero, as well as data taken immediately before the gap between the two orbits. For sector 39, we used the Lightkurve⁵ Python software to extract light curves with a 10 minutes cadence from the full frame images (FFI). The background level was evaluated in three different ways, either restricting ourselves to those pixels in a 50 × 50 pixels cutout that were below the median flux, or applying a principal component analysis including two or five components. We discarded the observations taken after 2 459 386 because of the large dispersion. This results in a total of 11 944 and 3411 data points for sectors 12 and 39, respectively, with a formal photometric accuracy ranging between 0.34 and 0.38 mmag. Since HD 152219 was located near the edge of one of the CCDs of the TESS instrument during the sector 12 observations, we also performed an extraction of the light curve using the FFI with a 30 minutes cadence. For this purpose, we used the Lightkurve software, restricting the source region to the central (21″)² pixel. The background level was evaluated as for sector 39. The results were identical within the errors, and in good agreement with those obtained from the 2 minutes cadence PDC photometry.

Since the region around HD 152219 is rather crowded, both the TESS and the ASAS-3 photometry of our source are contaminated by nearby sources. Within a radius of 1′, we found two stars with a magnitude difference in the I-band of less than 4 mag compared to HD 152219 (m_I = 7.27). The nearest contaminator is CPD−41° 7706 (V 945 Sco), a β Cep variable located at 38.5″ with m_I = 9.26 according to Sung et al. (1998). The second contaminator, CPD−41° 7712, is located at 54.0″ and has m_I = 8.84 (Sung et al. 1998). Assuming all the light from both sources contaminates the photometry of HD 152219, we expect a third light fractional contribution of ∼0.28. If instead, the contamination only comes from the nearest neighbour, the third light contribution would amount to ∼0.16. The masks adopted for the TESS light curves extraction (as specified in the CROWDSAP keyword of the PDC files) and the full frame images suggest that the likely third light contribution due to contamination by these two sources should be about 0.09. Whilst the PDC pipeline corrects the photometry for this level of third light, this is not the case for the FFI photometry, which might thus still contain a non-zero third light contribution.

Using the Fourier method of Heck et al. (1985) and Gosset et al. (2001), we analysed the periodogram of the TESS data to search for signals with frequencies up to 30 d⁻¹ (see Fig. 2). At frequencies up to 4 d⁻¹, the periodogram is largely dominated by HD 152219’s orbital frequency ν_orb = 0.2358 d⁻¹ and a number of its harmonics: frequencies up to 7 ν_orb are clearly seen, and even 9 ν_orb, 13 ν_orb, and 15 ν_orb are clearly detected. As outlined above, the TESS light curve could be contaminated by third light from two neighbouring sources, one of them is the β Cep variable V 945 Sco (Balona & Shobbrook 1983). Studies of this star consistently reported a signal near 14.91 d⁻¹ with an amplitude near 5.5 mmag in the V-band (Balona & Shobbrook 1983; Meingast et al. 2013). Given the magnitude difference between V 945 Sco and HD 152219, this signal should have an amplitude in the combined light of at most 0.8 mmag. In our periodogram of the TESS PDC photometry, the amplitude at this frequency is ∼0.5 mmag, which is consistent with the white noise level of the periodogram. Whilst we cannot totally rule out that the contaminating β Cep variable contributes to the non-orbital variability, its contribution must be very small, most probably comparable to or lower than the intrinsic variations of the OB-stars in HD 152219 and any residual instrumental effects. To reduce the amplitude of these non-orbital variations, we have constructed a mean light curve consisting of the normal points computed for phase bins of 0.01. Whilst this phase bin size does not completely remove the non-orbital variations, we note that adopting a larger bin size would smear out the details of the shape of the eclipses and thus impact the outcome of the analysis in Sect. 5. The TESS light curves were obtained over sufficiently short time intervals to ensure that we can assume a constant ω over the duration of each TESS sector.

Fig. 2. Fourier periodogram of the TESS data from sector 12. The insert zooms on the low frequency domain. The red lines yield the orbital frequency νorb = 0.2358 d−1 and its first six harmonics.

3. Spectral analysis

3.1. Spectral disentangling

The FEROS spectra were used as input to our disentangling code, which is based on the method described by González & Levato (2006), to establish the individual spectra of the binary components as well as their RVs at the times of the observations. We only considered those 93 spectra taken outside of the photometric eclipses to avoid biasing the reconstruction of the spectra of the binary components. The method as well as its limitations are described in details in Rosu et al. (2020a, Sect. 3.1).

For the synthetic TLUSTY spectra used as cross-correlation template for the RV determination, we assumed T_eff = 30 000 K, log g = 3.75, and v sin i = 160 km s⁻¹ for the primary star and T_eff = 21 000 K, log g = 4.00, and v sin i = 160 km s⁻¹ for the secondary star. We performed 60 iterations to disentangle the spectra and determine the RVs. The first 30 iterations were performed with the RVs fixed to the initial input values estimated from Gaussian fits of the He Iλ 4026 line. The last 30 iterations were then performed allowing the RVs to vary as we found that this number of iterations was necessary to achieve a relative error on the derived RVs of less than 1% between the last two iterations. To check the robustness of the results of the disentangling against any bias that the initial approximation of the spectrum of star B, as a featureless continuum, could introduce, we repeated the disentangling procedure by interchanging the roles of stars A and B. The agreement was very good and the final reconstructed spectra were taken to be the mean of the two approaches.

The spectral disentangling was performed separately over a number of wavelength domains: A0[3800:3920] Å, A1[3990:4400] Å, A2[4300:4570] Å, A3[4570:5040] Å, A4[5380:5860] Å, A5[5830:6000] Å, A6[6400:6750] Å, and A7[7000:7100] Å. As outlined in Rosu et al. (2020a), the presence of interstellar lines or diffuse interstellar bands close to spectral lines in some of these spectral domains (A0 and A4) affects the quality of the resulting reconstructed spectra. In these cases, the disentangling code erroneously associated some of the line flux of the non-moving interstellar medium (ISM) lines to the stars or vice versa. In addition, not enough spectral lines – especially lines of the secondary star – are present in the A5 and A7 domains. Both situations affect the determination of the stellar RVs. Therefore, we first processed the wavelength domains (A1, A2, A3, and A6) for which the code was able to reproduce the individual spectra and simultaneously estimate the RVs of the stars. We then computed the mean of the stellar RVs for each observation. The RVs from the individual wavelength domains were weighted according to the number of strong lines present in these domains (five lines for A1, three for A2, four for A3, and two for A6). We note that the Hγ and He Iλ 4388 lines are present in both A1 and A2 domains, meaning that their weight is artificially doubled in the computation of the mean RVs. Solving this issue by lowering the weight of A1 and/or A2 gives mean RVs which differ by less than 1%. Hence, we decided to keep the weight consistent with the number of lines in the domains. The resulting RVs of both stars are reported in Table A.1 together with their 1σ errors. We finally performed the disentangling on the four remaining domains (A0, A4, A5, and A7) with RVs fixed to these weighted averages, and, for the A0 and A4 domains only, using a version of the disentangling code designed to deal with non-moving ISM lines (Mahy et al. 2012).

Rosu et al. (2020a) outlined several limitations of the disentangling method for massive binaries. First of all, spectral disentangling can only reconstruct features that are well sampled by the Doppler excursions during the orbital cycle. If broad lines, such as H I Balmer lines, only partially deblend over the orbital cycle, this might lead to artefacts in the wings of these lines. Yet, in the case of HD 152219, the RV amplitudes are sufficient to ensure a reliable reconstruction of the wings of the Balmer absorption lines. Furthermore, there is no indication of a significant emission in the He IIλ 4686 and Hα lines arising from a wind interaction zone which could bias the reconstruction of these lines. We thus conclude that the reconstructed spectra of the components of HD 152219 should be free of such biases.

3.2. Spectral classification and absolute magnitudes

The reconstructed spectra of the binary components of HD 152219 allowed us to reassess the spectral classification of the binary components. For the primary star, we used the criterion of Conti & Alschuler (1971) complemented by Mathys (1988) to determine the spectral type of the star. We found that log W′ = log[EW(He Iλ 4471)/EW(He IIλ 4542)] amounts to 0.61 ± 0.01, which corresponds to spectral type O9.5. Hereafter, errors are given as ±1σ. In addition, we used the criteria of Sota et al. (2011) and Sota et al. (2014) based on the ratio between the strengths of several lines: the ratio He IIλ 4542/He Iλ 4388 amounts to 0.67 ± 0.01, the ratio He IIλ 4200/He Iλ 4144 amounts to 0.88 ± 0.01, and the ratio Si IIIλ 4552/He IIλ 4542 amounts to 0.48 ± 0.01. All three ratios being lower than unity, the spectral type O9.7 is clearly excluded and the spectral type O9.5 is confirmed. We further applied the criteria of Martins (2018) based on ratios of EWs. The ratio EW(He Iλ 4144)/EW(He IIλ 4200) amounts to 1.03 ± 0.01, suggesting a spectral type O9.2 with spectral types O9, O9.5, and O9.7 within the error bars; the ratios EW(He Iλ 4388)/EW(He IIλ 4542) and EW(Si IIIλ 4552)/EW(He IIλ 4542) amount to 1.81 ± 0.01 and 0.46 ± 0.01, respectively, both suggesting an O9.5 spectral type, excluding the O9.2 spectral type but including the O9.7 spectral type within the error bars. To assess the luminosity class of the primary star, we used the criterion from Conti & Alschuler (1971). We found that log W′ = log[EW(Si IVλ 4089)/EW(He Iλ 4143)] amounts to 0.21 ± 0.01, which corresponds to luminosity class III. The relative strength of the He Iλ 4026 and Si IVλ 4089 lines, as well as of the He IIλ 4686 and He Iλ 4713 lines, suggests a luminosity class between V and III according to the criteria and spectral atlas of Sota et al. (2011). Similar conclusions are drawn from the criteria of Martins (2018): The ratios EW(He IIλ 4686)/EW(He Iλ 4713) and EW(Si IVλ 4089)/EW(He Iλ 4026) amount to 1.77 ± 0.01 and 0.50 ± 0.01, respectively, both suggesting a luminosity class III, excluding luminosity classes IV and II. We determine that the primary star is of spectral type O9.5 III, which agrees with the determination of Sana et al. (2006). For the secondary star, we used the atlases and criteria of Gray (2009), Walborn & Fitzpatrick (1990) and Liu et al. (2019). We first observed that the He IIλ 4200 line is absent in the spectra, hence excluding spectral-type O. Qualitatively, we used the ratio between the strengths of the He Iλ 4471 and Mg IIλ 4481 lines, which suggests that the spectral type cannot be later than B2.5. This statement is reinforced by the fact that the Si IIλλ 4128-32 lines are only marginally detected, if at all, whilst they should be clearly visible for a spectral type later than B2.5. Furthermore, the presence of C IIIλ 4650, N IIIλ 4097, and Si IVλλ 4089, 4116 lines suggests a spectral type no later than B2. The fact that the Si IVλ 4089 line is clearly visible but weaker than the Si IIIλ 4552 line suggests a spectral type B1 with an uncertainty of one spectral type. The ratio of the strengths of the Si IIIλ 4552 and He Iλ 4387 lines suggests a main-sequence or giant luminosity class (Walborn & Fitzpatrick 1990). Furthermore, the O IIλ 4348 line is nearly absent, O IIλ 4416 is weak and the Balmer lines are broad, as expected for luminosity classes V or III. In conclusion, the classification criteria suggest that the secondary star should be a B1-2 V-III star, again agreeing with Sana et al. (2006).

The brightness ratio in the visible domain was estimated based on the ratio between the EWs of the spectral lines of the secondary⁶ star and TLUSTY synthetic spectra of similar effective temperatures. For this purpose, we used the Hβ, He Iλλ 4026, 4144, 4471, 4921, 5016, and 5876 lines and TLUSTY synthetic spectra of effective temperatures equal to 20 000, 21 000, and 22 000 K. The ratio EW_TLUSTY/EW_sec = (l₁ + l₂)/l₂ is equal to 10.95 ± 2.07, 10.77 ± 2.02, and 10.45 ± 1.89 for TLUSTY synthetic spectra of 20 000, 21 000, and 22 000 K, respectively. As the three results are very close, we adopt the value for a TLUSTY spectrum of 21 000 K and infer a value for l₂/l₁ of 0.101 ± 0.021.

The second Gaia data release (DR2, Gaia Collaboration 2018) quotes a parallax of ϖ = 0.644 ± 0.084 mas, corresponding to a distance of ${1507}_{-178}^{+231}$ pc (Bailer-Jones et al. 2018). From there we derived a distance modulus of $10.{89}_{-0.26}^{+0.33}$. Sung et al. (1998) reported mean V and B magnitudes of 7.552 ± 0.013 and 7.705 ± 0.035, respectively. Adopting a value of −0.26 ± 0.01 for the intrinsic colour index (B − V)₀ of an O9.5 V-III star (Martins & Plez 2006) and assuming the reddening factor in the V-band R_V equal to 3.3 ± 0.1 (Sung et al. 1998), we inferred an absolute magnitude of the binary system $M_V=-4.{70}_{-0.36}^{+0.29}$. The brightness ratio then yields individual absolute magnitudes $M_{V,1}=-4.{60}_{-0.36}^{+0.29}$ and $M_{V,2}=-2.{12}_{-0.41}^{+0.35}$ for the primary and secondary stars, respectively. Comparing with the typical magnitudes reported by Martins & Plez (2006), we find that M_V, 1 is between the values of an O9.5 V and an O9.5 III star.

These results, along with the radii determined in Sana et al. (2006) and in Sect. 5 below, seem to confirm that contrary to what the luminosity class criteria suggest, the primary star would be a main-sequence star rather than a giant or subgiant. The reason is likely the fact that the surface gravity log g is not only set by the ratio of mass over radius squared, but reflects the gradient of the binary potential evaluated at the surface of the star. Likewise, comparing the magnitude obtained for the secondary star to those reported by Humphreys & McElroy (1984) suggests that the luminosity classes IV and III can clearly be excluded. Instead, the secondary star’s absolute magnitude is coherent within the error bars with the magnitude of a B2 V star (Humphreys & McElroy 1984).

3.3. Projected rotational velocities

The projected rotational velocities of both stars were derived using the Fourier transform method (Simón-Díaz & Herrero 2007; Gray 2008), which has the advantage of being able to separate the effect of rotation from other broadening mechanisms such as macroturbulence. We applied this method to a set of well-isolated spectral lines, which are therefore expected to be free of blends. The results are presented in Table 1, and the Fourier transforms are illustrated for the C IIλ 4267 line in Fig. 3 for the primary and secondary stars. It is commonly recommended to use metal lines to compute the projected rotational velocity. The H I, He I, and He II line profiles are indeed more affected by non-rotational broadening mechanisms such as Stark broadening, which alter the position of the first zero in the Fourier transform (Simón-Díaz & Herrero 2007). However, when the projected rotational velocity of the star exceeds 100 km s⁻¹, the effect of Stark broadening on the He I lines is expected to be small and these lines can then be used as a complement to the metal lines. The results presented in Table 1 show that the mean v sin i_rot computed on metal lines alone or on all the lines agree very well. We adopted a mean v sin i_rot of 166 ± 10 km s⁻¹ for the primary star and of 95 ± 9 km s⁻¹ for the secondary star. For the primary star, our results agree with the values obtained by Levato & Morrell (1983) and Sana et al. (2006), but are lower than those obtained by Conti & Ebbets (1977) (see Table 1).

Fig. 3. Top row: C IIλ 4267 line profiles of the separated spectra obtained after application of the brightness ratio for the primary (left panel) and secondary (right panel) stars. Bottom row: Fourier transform of those lines (in black) and best-match rotational profile (in red) for the primary (left panel) and secondary (right panel) stars.

3.4. Wind terminal velocity

To estimate the wind terminal velocity v_∞ of the primary star of HD 152219, we extracted seven high-resolution spectra taken with the short wavelength primary (SWP) camera of the International Ultraviolet Explorer (IUE) satellite. The most prominent wind feature in these spectra is clearly the saturated C IVλλ 1548, 1551 P-Cygni profile. We measured the velocity of the violet edge of zero residual intensity in the absorption trough v_black (see Fig. 4). In this way, we derived v_black = (1930 ± 120) km s⁻¹. The dispersion most probably arises (at least partially) from the binarity of HD 152219. Indeed, the interaction of the primary wind with either the secondary star or the wind of the latter leads to a reduction of the primary wind velocity in some directions (and thus of v_black at some orbital phases). Since v_black was shown to be a good indicator of v_∞ for OB stars (Prinja et al. 1990), we thus adopted v_∞ = 1930 km s⁻¹.

Fig. 4. Determination of vblack on the absorption trough of C IVλλ 1548, 1551 as observed in the IUE spectra of HD 152219. The black and red curves correspond to observations SWP 56939 and SWP 56953, respectively.

3.5. Model atmosphere fitting

The reconstructed spectra of the binary components were analysed by means of the CMFGEN model atmosphere code (Hillier & Miller 1998) to constrain the fundamental properties of the stars. This code solves the radiative transfer equations, as well as the equations of statistical and radiative equilibrium in the comoving frame assuming the star and its wind are spherically symmetric. Line-blanketing and clumping are included in the code. For the wind, a standard β-velocity law is adopted through the relation

$$\begin{array}{c}\hfill v(r)=v_{\infty }{(1-\frac{R_{\ast }}{r})}^{\beta },\end{array}$$(1)

where v_∞ is the terminal wind velocity, R_* the radius of the star, and r the radial distance from the centre of the star. The wind clumping is introduced in the models using a two-parameter exponential law for the volume filling factor

$$\begin{array}{c}\hfill f(r)=f_1+(1-f_1)exp(-\frac{v(r)}{f_2}),\end{array}$$(2)

where the parameters f₁ and f₂ are the filling factor value when r → ∞ and the onset velocity of clumping, respectively. Including clumping in the models leads to a reduced estimate of the mass-loss rate compared to an unclumped model. To compare with values obtained with unclumped models, the mass-loss rate derived in this paper has to be multiplied by a factor of $1/\sqrt{f_1}$ (Martins 2011). Regarding the microturbulence velocity v_micro, CMFGEN assumes that it depends on the position r as

$$\begin{array}{c}\hfill v_{\text{micro}}=v_{\text{micro}}^{\text{min}}+(v_{}^{\text{micro}}-v_{\text{micro}}^{\text{min}})\frac{v(r)}{v_{\infty }},\end{array}$$(3)

where $v_{\text{micro}}^{\text{min}}$ and $v_{\text{micro}}^{\text{max}}$ are the minimum and maximum v_micro. The value of $v_{\text{micro}}^{\text{max}}$ is fixed to 0.1v_∞, whilst $v_{\text{micro}}^{\text{min}}$ is adjusted on metal lines.

The CMFGEN spectra were first broadened by the projected rotational velocities determined in Sect. 3.3. The stellar and wind parameters were then adjusted following the procedure outlined by Martins (2011). We proceeded slightly differently for the two stars as they have quite different spectral types.

3.5.1. Primary star

The macroturbulence velocity was adjusted on the wings of the O IIIλ 5592 and Balmer lines and we derived a value of v_macro = 120 ± 10 km s⁻¹. We adjusted the microturbulence velocity on the metal lines and inferred a value of 15 ± 3 km s⁻¹.

The effective temperature was determined by searching for the best overall fit of the He I and He II lines. This was clearly a compromise as we could not find a solution that perfectly fits all helium lines simultaneously. For instance, we found that the strength of the He Iλ 4471 line is underestimated while other He I lines are well-adjusted. Given the luminosity class of the primary, it seems unlikely that this discrepancy reflects the dilution effect discussed by Voels et al. (1989) and Herrero et al. (1992). Discarding some lines that are severely affected by blends (such as He Iλ 4121 blended with Si IVλ 4116 or He IIλ 5412 possibly affected by blends with diffuse interstellar bands), or turned out to be rather insensitive to a small variation in the effective temperature (such as He Iλ 4026 which reaches its maximum intensity at spectral type O9 and is also blended with the weak but non-zero He IIλ 4026 line), we finally adjusted the effective temperature on the He Iλ 4922 line and inferred a value of 30 900 ± 1000 K. With this effective temperature, the He Iλλ 4026, 4388, 4713, and 4922 lines as well as the He IIλ 4200 line are perfectly adjusted. A possible explanation for the difficulties encountered when attempting to determine the effective temperature could be the fact that we are applying a 1D model atmosphere code that assumes a spherical geometry to a spectrum that forms in the atmosphere of a distorted star in a binary system. Gravity darkening then leads to a highly non-uniform temperature distribution (Palate & Rauw 2012), which is not accounted for in the model atmosphere code.

The surface gravity was obtained by adjusting the wings of the Balmer lines Hβ, Hγ, Hδ, and H Iλλ 3835, 3890. We excluded Hε because this line is not correctly reproduced by the disentangling code owing to the blend with the interstellar Ca II line. We derived log g_spectro = 3.60 ± 0.10.

The surface chemical abundances of all elements, including helium, carbon, nitrogen, and oxygen, were set to solar as taken from Asplund et al. (2009). Indeed, HD 152219 is a typical case where it is impossible to certify that the chemical abundances of the elements differ from solar; the number of parameters we have to adjust exceeds the number of lines we can use for this purpose. For instance, to adjust the oxygen abundance, the O IIIλ 5592 line is the only one free of blends that could be used. To perfectly reproduce this line, an oxygen abundance nearly twice solar would be required, which is highly unlikely. Indeed, for such an oxygen abundance the synthetic CMFGEN spectrum predicts a series of O III absorption lines (λλ 4368, 4396, 4448, 4454, and 4458) that are not present in the observed spectra (neither before nor after disentangling). Moreover, this would also artificially reduce the carbon abundance, inferred from the C IIIλ 4070 line. This C III line is the only suitable carbon line which is not affected by subtle formation processes controlled by a number of unconstrained far-UV lines. However, it is heavily blended with an oxygen line. As a result, and since HD 152219 is a relatively unevolved system, we consider that adopting solar abundances for the chemical elements is the best option.

Regarding the wind parameters, the clumping parameters were fixed: The volume filling factor f₁ was set to 0.1, and the f₂ parameter that determines the wind velocity (and thus the position) where clumping starts was set to 100 km s⁻¹, as suggested by Hillier (2013). For the sake of completeness, we varied f₁ from 0.05 to 0.2 and f₂ from 80 to 120 km s⁻¹, and we did not observe any significant difference. Likewise, the β parameter of the velocity law was fixed to the value of 0.9 as suggested by Muijres et al. (2012) for an O9 V type star. Again for the sake of completeness, we tested a value of 1.0 for β but we did not observe any significant difference in the resulting spectrum.

In principle, the wind terminal velocity could be derived from the Hα line. However, because of the degeneracy between the wind parameters, we decided to fix the value of v_∞ to 1930 km s⁻¹ as determined in Sect. 3.4. Regarding the mass-loss rate, the main diagnostic lines in the optical domain are Hα and He IIλ 4686. We inferred a value of Ṁ = 6 × 10⁻⁸ M_⊙/yr based on these lines.

3.5.2. Secondary star

As the O IIIλ 5592 line is not present in the spectrum of the secondary star, we had to rely on He I and Balmer lines to estimate the macroturbulence velocity. Hence, we determined v_macro, the effective temperature, and the surface gravity simultaneously, in an iterative manner, in order to get the best-fit of the spectrum. We inferred a value of v_macro = 50 ± 10 km s⁻¹.

As for the primary star, we adjusted the effective temperature mainly based on the He Iλ 4922 line, and found a value of 23 500 ± 1000 K. With this effective temperature, He Iλλ 4026, 4471, 5016, and 5876 are overestimated whilst He Iλλ 4388, 4713, and 7065 are underestimated.

In principle, the surface gravity is obtained by adjusting the wings of the Balmer lines. However, in the present case, we found that the wings, as well as the depths, of the Balmer lines are systematically overestimated. At least for some of the Balmer lines (e.g., Hδ and Hγ), this discrepancy very likely reflects residual normalisation errors which affect the disentangled spectra and are amplified by the large brightness ratio that is applied to renormalise the secondary spectrum. Hence, we decided to adopt a log g_spectro = 3.92 ± 0.10 as inferred by Sana (2009). With this value, the wings of most hydrogen lines are overestimated.

We fixed the microturbulence velocity v_micro to 10 km s⁻¹, a reference value for such types of stars as suggested by Cazorla et al. (2017). For the same reasons as for the primary star, we set the surface chemical abundances to solar as taken from Asplund et al. (2009).

Regarding the wind parameters, the clumping parameters f₁ and f₂ were fixed as for the primary star. We fixed the β parameter of the velocity law to 1.1 as representative of early B-type stars (see, e.g., Lefever et al. 2010), and did not find any significant difference when varying the value of this parameter.

We adjusted the mass-loss rate based on the Hγ and Hα lines and inferred a value of Ṁ = 1.0 × 10⁻¹⁰ M_⊙/yr. The wind terminal velocity was then adjusted based on the strength of the Hα line. We found a value of v_∞ = 2250 ± 50 km s⁻¹.

The stellar and wind parameters of the best-fit CMFGEN model are summarised in Table 2. The normalised disentangled spectra of both components of HD 152219 are illustrated in Fig. 5 along with the best-fit CMFGEN adjustment.

Fig. 5. Normalised disentangled spectra (in black) of the primary and secondary stars of HD 152219 (we note that the spectrum of the secondary star is shifted by +0.3 in the y-axis for convenience) together with the respective best-fit CMFGEN model atmosphere (in red).

3.5.3. Spectroscopic radii and masses

The bolometric magnitudes of the stars $M_{\text{bol},1}=-7.{46}_{-0.37}^{+0.31}$ and $M_{\text{bol},2}=-4.{17}_{-0.43}^{+0.38}$ were computed assuming that the bolometric correction depends only on the effective temperature through the relation

$$\begin{array}{c}\hfill \mathrm{BC}=-6.89log(T_{\text{eff}})+28.07\end{array}$$(4)

(Martins & Plez 2006). We further obtained bolometric luminosities $L_{\text{bol},1}=7.{74}_{-2.18}^{+2.66}\times {10}^4{L}_{\odot }$ for the primary star and $L_{\text{bol},2}=3.{72}_{-1.29}^{+1.48}\times {10}^3{L}_{\odot }$ for the secondary star. From the relation between bolometric luminosity, effective temperature, and radius, we inferred spectroscopic radii $R_{\text{spectro},1}=9.7_{-1.5}^{+1.8}{R}_{\odot }$ for the primary star and $R_{\text{spectro},2}=3.7_{-0.7}^{+0.8}{R}_{\odot }$ for the secondary star.

We corrected the surface gravities determined with CMFGEN to account to first order for the impact of the centrifugal force and of the radiation pressure. In this way, we obtained log g_C of 3.70 ± 0.08 and 3.95 ± 0.09 for the primary and secondary stars, respectively. From there, we then inferred spectroscopic masses $M_{\text{spectro},1}=17.2_{-6.3}^{+7.1}{M}_{\odot }$ for the primary star and $M_{\text{spectro},2}=4.4_{-1.9}^{+2.1}{M}_{\odot }$ for the secondary star. In Sect. 5, we compare the spectroscopic radii and masses with model-independent determinations based on the RV curves and the photometric light curve. Beside the difficulties outlined above to determine the surface gravity of the secondary star, it is important to recall that CMFGEN assumes that the stars are spherical, static, and isolated. However, because of the binarity, these three assumptions are not verified. The effective temperature is thus not homogeneous all over the surface and the computed temperature has to be considered as a mean value over the star. In addition, because of the companion, the surface gravity |∇(Ω)| is not homogeneous over the stellar surface and the local surface gravity is lower than that of an isolated star of same mass and radius (Palate & Rauw 2012).

4. Apsidal motion

Our total dataset of RVs consists of sixteen primary RVs taken from Hill et al. (1974), one from Conti et al. (1977), four from Levato & Morrell (1983), three from Perry et al. (1990), eight from García & Mermilliod (2001), eight from Stickland & Lloyd (2001), and the six CAT+CES RV points from Sana et al. (2006). These data were complemented by our 93 RV points obtained as part of the disentangling process of the FEROS observations of HD 152219. The third datapoint from the Perry et al. (1990) dataset (HJD 2 439 958.773) was discarded as it was found to be offset by about 80 km s⁻¹ from the most-likely orbital solution (see also Sana et al. 2006). Likewise, we found that the four last RV points (HJD 2 450 593.600 – 2 450 598.792) given by García & Mermilliod (2001) were discrepant. These points were thus also discarded. We therefore ended up with a series of 134 primary RV data points spanning about 38 years.

Sana et al. (2006) and Sana (2009) reported significantly different RV amplitudes for different lines of the primary star. For instance, Sana (2009) found a difference of 15% between the smallest primary RV amplitude (107.7 km s⁻¹ for the He Iλ 4922 line) and the largest value (124.4 km s⁻¹ for the He IIλ 4686 line). In contrast, the RVs inferred from our disentangling method applied to different wavelength domains did not show any significant differences. We therefore suspect that the amplitude differences reported by Sana et al. (2006) and Sana (2009) arise from the double Gaussian fitting method used by these authors to establish the RVs.

For those literature references that explicitly quote errors on the RVs, we adopted the latter values. For the IUE data from Stickland & Lloyd (2001) and for the CAT/CES data of Sana et al. (2006), we adopted errors of 5 km s⁻¹ as representative of the measurement error on these values based on the spectral resolution and the method used to derive the RVs. Moreover these values are very close to the dispersions of the RV points over the best orbital solutions. For the RVs derived from the spectral disentangling, we estimated the errors from the dispersion of the RVs determined from the various spectral domains.

Since the secondary RVs are only available for the most recent data and are subject to larger errors than the primary RVs, we did not use them in the determination of the rate of apsidal motion. For each time of observation t, we adjusted the primary RV data with the following relation

$$\begin{array}{c}\hfill {\mathrm{RV}}_{\mathrm{P}}(t)={\gamma }_{\mathrm{P}}+K_{\mathrm{P}}[cos(\varphi (t)+\omega (t))+ecos\omega (t)],\end{array}$$(5)

where γ_P, K_P, e, and ω(t) are respectively the apparent systemic velocity, the semi-amplitude of the RV curve, the eccentricity, and the argument of periastron of the primary orbit. The true anomaly ϕ(t) is evaluated through Kepler’s equation which involves e and the anomalistic orbital period P_orb. We explicitly accounted for secular variations of the argument of periastron, by adopting

$$\begin{array}{c}\hfill \omega (t)={\omega }_0+\dot{\omega }(t-T_0),\end{array}$$(6)

where ω₀ is the value of ω at the epoch T₀ of periastron passage and $\dot{\omega }$ is the rate of secular apsidal motion. As the RVs of different spectral lines potentially yield slightly different apparent systemic velocities (see Table 3 of Sana et al. 2006), and since literature RVs were obtained from diverse sets of lines, we adjusted the systemic velocity of each subset of our total dataset so as to minimise the sum of the residuals of the data about the curve given by Eq. (5) for each combination of the six model parameters. For our best-fit solution, the systemic velocity amounts to −24.3 ± 1.6, −42.9 ± 3.8, −13.0 ± 2.9, −35.0 ± 3.5, −24.6 ± 2.1, −42.7 ± 1.8, −6.7 ± 2.0, and −24.5 ± 0.1 km s⁻¹ for the Hill et al. (1974), Conti et al. (1977), Levato & Morrell (1983), Perry et al. (1990), García & Mermilliod (2001), Stickland & Lloyd (2001), Sana et al. (2006), and our subset, respectively.

To find the values of the six free parameters (K₁, P_orb, e, T₀, ω₀, and $\dot{\omega }$) that provide the best-fit to the whole set of RV data, we scanned the parameter space in a systematic way. Figure 6 illustrates the projections of the 6D parameter space onto 2D planes and the best-fit values are given in Table 3. The 1σ confidence contours are computed in each two-dimensional plane adopting Δχ² = 2.30 as appropriate for two free parameters. Figure 7 illustrates the best fit of the RV data at all epochs.

Fig. 6. Confidence contours for the best-fit parameters obtained from the adjustment of the full set of 134 primary RV data of HD 152219 with Eqs. (5) and (6). The best-fit solution is shown in each panel by the black filled dot. The corresponding 1σ confidence level is shown by the blue contour.

Fig. 7. Comparison between the measured RVs of the primary (filled dots) and secondary (open dots, when available) with the RVs expected from relations (5) and (6) with the best-fit parameters given in Table 3.

The value of $\dot{\omega }$ that best fits our RV time series is (1.198 ± 0.300)° yr⁻¹. Mayer et al. (2008) derived a rate of apsidal motion of (1.64 ± 0.18)° yr⁻¹ based on the then available literature RVs and ASAS-3 photometry. Within the errorbars, the two values overlap.

The primary and secondary radial velocities of an SB2 binary are related to each other through a linear relation

$$\begin{array}{c}\hfill {\mathrm{RV}}_{\mathrm{P}}(t)=-q{\mathrm{RV}}_{\mathrm{S}}(t)+B,\end{array}$$(7)

where q is the mass ratio $(q=\frac{M_{\mathrm{S}}}{M_{\mathrm{P}}}$, M_P and M_S being the mass of the primary and secondary star, respectively) and B is a linear combination of the apparent systemic velocities of the two stars (B = γ_P + q γ_S). Applying this linear regression to the RVs obtained in the disentangling process, we found $q=\frac{M_{\mathrm{S}}}{M_{\mathrm{P}}}=0.413\pm 0.004$. We used this result to build an SB2 orbital solution for HD 152219. Except for two more discrepant points (at HJD 2 452 037.8 and 2 453 131.9), the secondary RVs were adjusted with ${\chi }_{\nu }^2=13.1$. Our best-fit parameters and their 1σ errors are listed in Table 3. We note that the mass ratio we found is higher than the value (q = 0.395 ± 0.003) inferred by Sana et al. (2006) from the He I lines. We further note that our best-fit eccentricity is lower than their value (e = 0.082 ± 0.011) but still compatible within the error bars, whilst their semi-amplitudes of the primary and secondary RV curves are respectively smaller (K_P = 110.7 ± 0.7 km s⁻¹) and larger (K_S = 279.9 ± 1.7 km s⁻¹), the latter being compatible within the error bars. Again these differences are likely due to the two-Gaussian fit used by Sana et al. (2006) to establish the RVs as well as to the fact that their solution did not account for the apsidal motion.

We further used Eq. (7) to convert the RVs of both stars into equivalent RVs of the primary star with

$$\begin{array}{c}\hfill {\mathrm{RV}}_{\text{eq}}(t)=({\mathrm{RV}}_{\text{P}}-q{\mathrm{RV}}_{\text{S}}(t)+B)/2,\end{array}$$(8)

whenever secondary RVs were available. For each time of observation t, we adjusted the equivalent RV data explicitly accounting for apsidal motion. The best fit is achieved for $\dot{\omega }=(1.69\pm 0.23)°$ yr⁻¹ and has ${\chi }_{\nu }^2=4.61$. This best-fit solution is not as good as the one based on the primary RVs only, as judged by the mean rms residual value, which is three times larger than the one obtained for the primary RVs fit. In addition, those results predict values for ω at the times of the TESS observations that are not compatible, within the error bars, with the values obtained from the fit of the TESS data. For those reasons, we kept the solution based on the adjustment of the primary RVs only as our best-fit solution.

5. Light curve analysis

We analysed the light curves of HD 152219 with the Nightfall binary star code (version 1.92) developed by R. Wichmann, M. Kuster, and P. Risse⁷ (Wichmann 2011). The shape of the stars is described by the Roche potential scaled with the instantaneous separation between the stars. In the absence of any surface spots, the eight model parameters are the mass ratio q, the effective temperatures T_eff,P and T_eff,S of the primary and secondary stars, the Roche lobe filling factors f_P and f_S of the primary and secondary stars, the eccentricity of the orbit e, the orbital inclination i, and the argument of periastron ω. The Roche lobe filling factor is defined as the ratio between the polar radius of the star and the associated polar radius of the Roche lobe at periastron. In the case of HD 152219, an additional parameter comes from the possible contribution of the third light discussed in Sect. 2.2. A quadratic limb-darkening law was adopted with the coefficients in the I band coming from Claret (2000). In view of the stellar effective temperatures, the gravity darkening exponent was set to 0.25 as appropriate for massive stars with radiative envelopes. The code accounts for reflection effects by taking into account the mutual irradiation of the stellar surface elements of both stars (Hendry & Mochnacki 1992).

Following the radial velocities analysis (see Table 3 and Sect. 4), we fixed q to 0.413 and e to 0.072. We fixed the effective temperature of the primary star to 30 900 K as derived from the CMFGEN analysis (see Table 2 and Sect. 3.5). In a first attempt to reproduce the TESS-12 light curve of the binary system, we also fixed the effective temperature of the secondary star to 23 500 K as found in the CMFGEN analysis. However, the solutions obtained in this way were unable to reproduce the depth of the primary eclipse. Considering both the difficulties encountered in the spectroscopic analysis to determine the effective temperature of the secondary star and the high uncertainty on the inferred value, we decided to leave T_eff,S free in the light curve analysis. For each adjustment, we therefore had five free parameters: T_eff,S, f_P, f_S, i, and ω. The results are given for the TESS-12 data in Table 4, together with the corresponding reduced ${\chi }_{\nu }^2$. We note that the secondary/primary light ratio in the V-band of 0.084, estimated by the Nightfall code, is compatible with the light ratio inferred from spectroscopy (see Sect. 3.2) within the error bars.

We scanned the parameter space of the five free parameters to estimate the error bars. For this purpose, we used the χ² mapping functionality implemented in Nightfall. This function scans two parameters at a time, fixing the other three to their best-fit value. We then computed the 1σ confidence contours in each two-dimensional plane adopting Δχ² = 2.30 as appropriate for two free parameters. The confidence contours obtained in this way are plotted in blue in Fig. 9. In this way, adopting the largest contours on the parameters to account for correlations, we obtained the values given in Table 4. For a near circular orbit seen under an inclination of 90°, the depths of the eclipses are mostly set by the ratio of the effective temperatures of the two stars. The primary star’s effective temperature, T_eff,P, determined in Sect. 3.5, is subject to uncertainties that hence propagate in the determination of T_eff,S. To quantify this dependence, we computed two light curves fixing all parameters to their best-fit value, setting T_eff,P to its value ±1σ (i.e. 31 900 and 29 900 K, respectively), and leaving T_eff,S as the only free parameter. We respectively obtained T_{eff, S} = 22 320 and 21 047 K, values which differ from the best-fit temperature by +623 and −650 K, respectively. Considering the results from the CMFGEN modelling (Sect. 3.5) and the light curve analysis, we hence adopted an error bar of ±1000 K for T_eff,S to be conservative.

We adopted the parameters of the best-fit to the TESS-12 data as our reference photometric solution. The best fit is illustrated in Fig. 8.

Fig. 8. Upper panel: best-fit Nightfall solution of the TESS-12 light curve of HD 152219. Lower panel: residuals over the best-fit solution.

Fig. 9. Confidence contours for the best-fit parameters obtained from the adjustment of the TESS-12 light curve of HD 152219 with the Nightfall code. The best-fit solution (i.e. with a third light contribution I3 = 0) is shown in each panel by the black dot, while the corresponding 1σ confidence level is shown by the blue contour.

The combination of the minimum semi-major axis and masses inferred from the radial velocity analysis (see Table 3) on the one hand and the best-fit inclination value on the other hand yields a semimajor axis of 32.81 ± 0.24 R_⊙ for the system and absolute masses for the primary and secondary stars M₁ = 18.64 ± 0.47 M_⊙ and M₂ = 7.70 ± 0.12 M_⊙. As expected from the discussion in Sect. 3.5.3, these masses are higher than the spectroscopic masses even though the primary masses are compatible within the error bars. From this solution, we inferred values for the primary and secondary stellar radii $R_1=9.{40}_{-0.15}^{+0.14}{R}_{\odot }$ and R₂ = 3.69 ± 0.06 R_⊙. Those radii are compatible within the error bars with the spectroscopic ones. This leads to photometric values of the surface gravities of the primary and secondary stars log g₁ = 3.76 ± 0.02 and log g₂ = 4.19 ± 0.01. As expected, these values are higher than those derived from the spectroscopic analysis (see discussion in Sect. 3.5). From the photometric stellar radii, the primary effective temperature derived in Sect. 3.5, and the secondary effective temperature derived here above, we inferred bolometric luminosities L_bol,1 = (7.26 ± 0.97)×10⁴ L_⊙ and L_bol,2 = (2.73 ± 0.51)×10³ L_⊙ for the primary and secondary stars, respectively. These values are compatible with those derived from the spectroscopic analysis (see Sect. 3.5.3). If we further assume the stellar rotational axes to be aligned with the normal to the orbital plane (but see Sect. 8), the combination of the orbital inclination, the stellar radii, and the projected rotational velocities of the stars derived in Sect. 3.3 yields rotational periods for the primary and secondary stars P_rot,1 = 2.86 ± 0.18 days and P_rot,2 = 1.97 ± 0.19 days. The ratio between rotational angular velocity and instantaneous orbital angular velocity at periastron amounts to 1.28 ± 0.08 and 1.86 ± 0.18 for the primary and secondary stars, respectively.

Mayer et al. (2008) inferred an orbital inclination i = 72 ± 3° from their analysis of the ASAS-3 photometric data. We note that our best-fits values are significantly larger, close to 90°. The total primary eclipse puts a stringent constraint on the inclination of the system. Hence, we checked that their inclination value of 72°, adopting for the other parameters our best-fit values, is not consistent with the total eclipse. To reproduce a total primary eclipse, an inclination of at least 80° is required.

The light curves of a binary system can in principle be used to determine the apsidal motion rate of the system (e.g., Zasche & Wolf 2019). However, in the case of HD 152219, the total time span of photometric observations is much smaller than the total time span of spectroscopic observations. Hence, we rather decided to use the photometric observations to perform a consistency check of the apsidal motion determined based on the spectroscopic data. To be consistent with the results obtained for the TESS-12 light curve, we adjusted the TESS-39 light curve keeping all parameters fixed to the values inferred for the TESS-12 except for ω. We found ω = 192.41° and ${\chi }_{\nu }^2$ = 0.500 for the best-fit adjustment. We observed that the depths of the eclipses were not perfectly reproduced: the theoretical depths were deeper than observed. This suggests that some third light contribution is present in the TESS-39 data and we hence adjusted the light curve keeping the third light contribution as a second free parameter. The best-fit adjustment was achieved for ω = 192.480° and I₃ = 0.032, has ${\chi }_{\nu }^2$ = 0.395, and is illustrated in Fig. 10. Adopting 1σ contours for the two-parameter space, we obtained $\omega =(192.5_{-6.1}^{+3.8})°$ and $I_3=0.{032}_{-0.014}^{+0.015}$.

Fig. 10. Upper panel: best-fit Nightfall solution of the TESS-39 light curve of HD 152219. Lower panel: residuals over the best-fit solution.

The resulting values of ω are provided in Table 5 together with the reduced ${\chi }_{\nu }^2$. Figure 11 illustrates the variations of ω with time. Both TESS-12 and TESS-39 data agree with the result inferred from the global fit of the RV data within the error bars. Unfortunately, the ASAS-3 data are too sparsely sampled to determine meaningful epoch-dependent values of ω.

Fig. 11. Values of ω as a function of time inferred from the photometric light curves and the RVs. The pink symbols correspond to the data of the fits of the TESS-12 and TESS-39 photometry. The blue dot indicates the ω0 value obtained from the global fit of all RV data. The solid blue line corresponds to our best-fit value of $\dot{\omega }$ inferred from the RVs, and the hatched cyan zone corresponds to the range of values according to the 1σ uncertainties on ω0 and $\dot{\omega }$.

We computed the phase difference between the primary and secondary eclipses. To do so, we adjusted a second-order polynomial to the eclipses for the TESS-12, TESS-39, and combined ASAS-3 data, and found the values of 0.451  ±  0.004, 0.454  ±  0.004, and 0.452  ±  0.010, respectively (see Fig. 12). These values were confirmed by computing the first order moment of the eclipses and from the times of minima determined by the Nightfall code. We adopted conservative error bars to take into account the difficulty to determine the minima of the secondary eclipse. We then computed the phase difference as a function of ω, defined as

$$\begin{array}{c}\hfill \mathrm{\Delta }\varphi =\frac{t_2-t_1}{P_{\text{orb}}}=\frac{\mathrm{\Psi }-sin\mathrm{\Psi }}{2\pi },\end{array}$$(9)

Fig. 12. Values of the phase difference Δϕ between the primary and secondary eclipses as a function of time and ω inferred from the photometric light curves and the RVs. The pink filled symbol and the pink open symbols correspond to the data of the fits of the ASAS-3 photometry, and the TESS-12 and TESS-39 photometry, respectively. The solid blue line corresponds to our best-fit value of Δϕ inferred from the RVs, and the hatched cyan zone corresponds to the range of values according to the 1σ uncertainties on e.

where t₁ and t₂ are the primary and secondary minima, respectively,

$$\begin{array}{c}\hfill \mathrm{\Psi }=\pi +2arctan\left(\frac{ecos\omega }{\sqrt{1-e^2}}\right),\end{array}$$(10)

and adopting e from the RV analysis. The resulting curve is plotted in Fig. 12. We observe that within the errorbars, the observational values of Δϕ agree with the curve showing the expected values from the RV analysis. Adopting the periastron length obtained from the RV analysis at the time of the TESS-12 observations, we find an eccentricity of $0.{070}_{-0.005}^{+0.006}$ that is compatible within the errorbars with an eccentricity of 0.072, therefore confirming the coherence between the RV and photometric analyses.

6. Internal structure constant k₂

The apsidal motion rate of a binary system is made of two contributions:

$$\begin{array}{c}\hfill \dot{\omega }={\dot{\omega }}_{\text{N}}+{\dot{\omega }}_{\text{GR}}.\end{array}$$(11)

The classical Newtonian contribution (N) is the dominant term in the majority of non-degenerate binary systems, except in a few cases where the general relativistic correction (GR) is the dominant term (see e.g., Baroch et al. 2021, and references therein).

The general relativity contribution to the rate of apsidal motion depends on the total mass of the binary, the orbital period, and the eccentricity of the system through the following expression

$$\begin{array}{c}\hfill {\dot{\omega }}_{\text{GR}}={\left(\frac{2\pi }{P_{\text{orb}}}\right)}^{5/3}\frac{3{(G(M_1+M_2))}^{2/3}}{c^2(1-e^2)},\end{array}$$(12)

when only the quadratic corrections are taken into account (Shakura 1985). In this expression, G and c stand for the gravitational constant and the speed of light, respectively. Using the eccentricity and the orbital period derived from the radial velocity analysis in Sect. 4 as well as the masses derived from the photometric analysis in Sect. 5, we inferred a value for ${\dot{\omega }}_{\text{GR}}$ of (0.159 ± 0.002)° yr⁻¹. As a consequence, the observational contribution of the Newtonian term ${\dot{\omega }}_{\text{N}}$ amounts to (1.039 ± 0.300)° yr⁻¹.

The Newtonian contribution takes the form adopted from Sterne (1939),

$$\begin{array}{cc}& {\dot{\omega }}_{\text{N}}=\frac{2\pi }{P_{\text{orb}}}[15f(e)\{k_{2,1}q{\left(\frac{R_1}{a}\right)}^5+\frac{k_{2,2}}{q}{\left(\frac{R_2}{a}\right)}^5\}\hfill \\ \hfill & +g(e)\{k_{2,1}(1+q){\left(\frac{R_1}{a}\right)}^5{\left(\frac{P_{\text{orb}}}{P_{\text{rot,1}}}\right)}^2\hfill \\ \hfill & +k_{2,2}\frac{1+q}{q}{\left(\frac{R_2}{a}\right)}^5{\left(\frac{P_{\text{orb}}}{P_{\text{rot,2}}}\right)}^2\}],\hfill \end{array}$$(13)

if the stellar rotation axes are orthogonal to the orbital plane, and where only the contributions arising from the second-order harmonic distortions of the potential are considered. In this expression, f and g are functions of the eccentricity of the system given by

$$\begin{array}{c}\hfill \{\begin{array}{cc}& f(e)=\frac{1+\frac{3e^2}{2}+\frac{e^4}{8}}{{(1-e^2)}^5},\hfill \\ \hfill & g(e)=\frac{1}{{(1-e^2)}^2}.\hfill \end{array}\end{array}$$(14)

The Newtonian term is itself the sum of the effects induced by stellar rotation and tidal deformation which in the present case amount to about 33% and 67% of the total Newtonian term, respectively. The main unknowns in Eq. (13) are k_2, 1 and k_2, 2, the internal structure constants of the primary and secondary stars. Also known as the apsidal motion constant, the internal structure constant k₂ of a star is expressed as

$$\begin{array}{c}\hfill k_2=\frac{3-{\eta }_2(R_{\ast })}{4+2{\eta }_2(R_{\ast })},\end{array}$$(15)

where

$$\begin{array}{c}\hfill {\eta }_2(R_{\ast })={\frac{dln{ϵ}_2}{dlnr}|}_{r=R_{\ast }}\end{array}$$(16)

is the logarithmic derivative of the surface harmonic of the distorted star evaluated at the stellar surface (r = R_*), expressed in terms of the ellipticity ϵ₂. η₂(R_*) is the solution of the Clairaut-Radau differential equation

$$\begin{array}{c}\hfill r\frac{d{\eta }_2(r)}{\mathit{dr}}+{\eta }_2^2(r)-{\eta }_2(r)+6\frac{\rho (r)}{\overline{\rho }(r)}({\eta }_2(r)+1)-6=0\end{array}$$(17)

with the boundary condition η₂(0) = 0 (Hejlesen 1987). The terms ρ(r) and $\overline{\rho }(r)$ are the density at distance r from the centre of the star and the mean density within the sphere of radius r, r being the current radius at which the equation is evaluated. Qualitatively, the internal structure constant is a measure of the internal mass distribution inside the star, that is to say the density contrast between the stellar core and external layers. An homogeneous sphere of constant density has its k₂ equal to 0.75 whilst a massive star having a very dense core and a diluted atmosphere can see its k₂ take a value as low as 10⁻⁴ (see e.g., Rosu et al. 2020b). During stellar evolution, the external layers become more and more diluted compared to the stellar core, hence k₂ decreases with time, making this quantity a good indicator of stellar evolution.

In the case of HD 152219, it is impossible to observationally constrain the individual values of the apsidal motion constants of both stars, as Eq. (13) is an underdetermined system. Nevertheless, we can define a weighted-average apsidal motion constant ${\overline{k}}_2$ for the binary system. For this purpose, we first rewrite Eq. (13) in the form

$$\begin{array}{c}\hfill {\dot{\omega }}_{\text{N}}=c_1{k}_{2,1}+c_2{k}_{2,2},\end{array}$$(18)

where c₁ and c₂ are two functions of known stellar parameters which expressions are straightforward from Eq. (13). Dividing Eq. (18) by (c₁ + c₂), we get

$$\begin{array}{c}\hfill {\overline{k}}_2=\frac{c_1{k}_{2,1}+c_2{k}_{2,2}}{c_1+c_2}=\frac{{\dot{\omega }}_{\text{N}}}{c_1+c_2}.\end{array}$$(19)

All terms appearing in the right-hand side of Eq. (19) are known from observations, hence we get ${\overline{k}}_2=0.00173\pm 0.00052$. As the secondary star is smaller and less massive than the primary star, its k₂-value exceeds that of the primary. Therefore, we have $k_{2,1}<{\overline{k}}_2<k_{2,2}$. Observationally, we have $\frac{c_1}{c_1+c_2}=0.95$ and $\frac{c_2}{c_1+c_2}=0.05$, meaning that the k_2, 1 has a much higher weight in the calculation of ${\overline{k}}_2$⁸. We further note that – and here we anticipate the analysis performed in Sect. 7.2 – if we make the hypothesis that $k_{2,1}={\overline{k}}_2$, and if k_2, 2 appears to be equal to n * k_2, 1, then the ensuing error beyond the initial hypothesis amounts to (n − 1)*5%.

7. Stellar structure and evolution models

We computed stellar evolution models with the Code Liégeois d’Évolution Stellaire⁹ (Clés, Scuflaire et al. 2008). An overview of the main features of Clés in the present context is presented in Rosu et al. (2020b, Sect. 2.1); we refer to this paper for further information and recall here only the most important points. For massive stars, Clés allows us to build stellar structure and evolution models from the Hayashi track to the asymptotic giant branch phase. The zero-age of a model is defined at the beginning of the pre-main sequence phase. The code adopts OPAL opacities from Iglesias & Rogers (1996) and solar chemical composition from Asplund et al. (2009). The FreeEOS equation of state and the rates of nuclear reactions are implemented following, respectively, Irwin (2012, a short description can be found in Cassisi et al. 2003) and Adelberger et al. (2011), whilst the mixing length theory is adopted to parameterise the mixing in convective regions (Cox & Giuli 1968). The stellar atmosphere is a radiative grey model atmospheres computed in the Eddington approximation, which is treated separately and added as a boundary condition. By default, internal mixing for massive stars is restricted to the convective core, but the user can introduce overshooting and/or turbulent mixing in the models (see, e.g., Rosu et al. 2020b). When overshooting is introduced in a model, the boundary of the central mixed region is displaced from its initial position given by the Ledoux criterion over a distance α_ovH_p. The temperature gradient is the radiative gradient in the overshooting region. Turbulent mixing can be introduced in the models through the inclusion of a pure diffusion term that reduces the abundance gradient of the considered element. It takes the following form

$$\begin{array}{c}\hfill -D_T\frac{\partial ln{X}_i}{\partial r},\end{array}$$(20)

where r is the radial coordinate, X_i is the mass fraction of element i at location r, and

$$\begin{array}{c}\hfill D_T=D_{\text{turb}}{\left(\frac{\rho }{{\rho }_0}\right)}^n+D_{\text{ct}}\end{array}$$(21)

is the turbulent diffusion coefficient (positive and expressed in cm² s⁻¹). In this expression, n, D_turb, and D_ct are three parameters chosen by the user. In the models, we did not include any microscopic diffusion. In the present analysis, we adopted a constant D_T throughout the star, that is to say we set n = 0.

The Clés code does not include any rotational mixing. Nonetheless, Rosu et al. (2020b) have shown that, to first order, the impact of rotational mixing on the internal structure can be simulated by means of the turbulent diffusion. Finally, the code is designed for single star evolution and does not account for the impact of binarity on the internal structure. However, Rosu et al. (2020b) have also shown that for an un-evolved binary system, the impact of binarity on the internal structures of the stars is quite modest and only slightly changes the age at which the star reaches a given stage of central condensation for a given external radius. Finally, the mass-loss rate induced by stellar winds is implemented through the Vink et al. (2001) recipe, to which the user can apply a multiplicative scaling factor ξ.

We built two sets of stellar evolution models for the primary and secondary stars of HD 152219. Whilst the primary star shows clear evidence of mass-loss through stellar wind, the mass-loss rate of the secondary star is much smaller. Hence, the initial stellar mass and the mass-loss prescription are only important for the primary star. In accordance with Asplund et al. (2009), we adopted as reference values the hydrogen mass fraction X = 0.715 and the metallicity Z = 0.015. These values are used to compute the models, unless otherwise stated (see Appendix B). The overshooting parameter is poorly constrained: Claret & Torres (2019) showed based on their MESA models that α_ov saturates around a value of 0.20 for stars more massive than 2 M_⊙ but their sample stopped at 4.5 M_⊙. More recently, Rosu et al. (2020b) explored values of α_ov up to 0.40 in their Clés models for the binary system HD 152248 located in the same cluster as our present target, and showed evidence for the need of high α_ov values. Tkachenko et al. (2020) reached the same conclusion from the analysis of their sample of stars based on MESA models. Martinet et al. (2021) suggested to use a value of α_ov = 0.20 for their GENEC rotating models for stars having a mass of 9 M_⊙ or higher. We emphasise here that different analyses using different stellar evolution codes with different prescriptions of the mixing processes reach the same qualitative conclusion that strong internal mixing is required.

For each model, the internal structure constant was computed solving the differential equation Eq. (17) using a Fortran routine in which a fourth-order Runge-Kutta method with step doubling is implemented. This code, validated against polytropic models, was already applied to the massive binaries HD 152218 (Rauw et al. 2016) and HD 152248 (Rosu et al. 2020b). This k₂ value is then corrected by the quantity

$$\begin{array}{c}\hfill \delta (log{k}_2)=-0.87{\lambda }_{\text{s}},\end{array}$$(22)

with

$$\begin{array}{c}\hfill {\lambda }_{\text{s}}=\frac{2{\mathrm{\Omega }}_{\text{s}}{2}^{}{R}_{\ast }^3}{3G{M}_{\ast }}\end{array}$$(23)

to account for the centrifugal deformation due to stellar rotation, following Claret (1999). In this expression, Ω_s stands for the observationally determined rotational angular velocity, computed from the v sin i_rot obtained in Sect. 3.3, and the stellar radius and the inclination determined in Sect. 5, while R_* and M_* stand for the radius and the mass of the corresponding model. In the following, unless explicitly specified, k_2, 1 and k_2, 2 are always corrected for the rotation of the star by Eq. (22).

We used the min-Clés Fortran routine¹⁰, in which a Levenberg-Marquardt minimisation method (Press et al. 1992) is implemented, to search for best-fit models. The routine determines the combination of the free parameters of the model allowing us to best reproduce the set of observationally determined present-day stellar properties. The initial masses and the age of the binary system are two major free parameters, also complemented by D_T in most cases. The set of currently observed properties consists of the masses, radii, and effective temperatures of the stars, as well as the weighted-average internal structure constant of the system: these parameters are summarised in Table 6. However, in min-Clés, only the first three parameters can directly be constrained in the models. Nonetheless, individual k₂ values can be constrained in the models (we return to this in Sect. 7.2). The min-Clés routine assumes symmetric uncertainties, hence we adopted error bars on each parameter given by the maximum value of the upper and lower bounds.

7.1. Preliminary analysis

In a first attempt to obtain stellar evolution models representative of both stars, we built two Clés models simultaneously, adopting the six corresponding constraints (M, R, and T_eff of both stars), and leaving the initial masses of the two stars and the age of the binary system as the only three free parameters. For both stars, we adopted ξ = 1, α_ov = 0.30, and no turbulent diffusion. The best-fit solution gives an age of 7.8 Myr, but none of the stellar parameters is correctly reproduced. Adding more free parameters would certainly help the convergence of the models but, due to degeneracies, might also lead to spurious solutions. Hence, we decided to proceed differently and to analyse both stars separately, keeping in mind that the ages of both stars have to match somehow.

As an illustration, we built two stellar evolutionary sequences, the first one with initial mass of 19.0 M_⊙ and ξ = 1 which should represent the primary star, and the second one with initial mass of 7.7 M_⊙ and ξ = 0.1 which should represent the secondary star. Both sequences assume α_ov = 0.30 and no turbulent diffusion. Whilst the stellar parameters of the primary star all cross the observational values, within the error bars, between 7 and 8 Myr, the stellar parameters of the secondary star cross the observational values, within the error bars, beyond 10 Myr. This preliminary analysis already shows the difficulty to build consistent stellar evolutionary models for both stars but also coherent in terms of age. Yet, the importance of these two evolutionary sequences rather lies in Fig. 13: The left panel shows the evolution as a function of the age of the internal structure constants of both evolutionary sequences together with the weighted-average mean of the k₂, while the right panel shows the evolution as a function of the age of the apsidal motion contributions of both stars taken separately as well as the total apsidal motion of the binary system. This second panel clearly shows that the primary star contributes for the most part to the binary apsidal motion, as expected. The evolution of ${\overline{k}}_2$ and k_2, 1 overlap nearly perfectly, suggesting that k_2, 1 dominates in ${\overline{k}}_2$. To quantify this assertion, we show in Fig. 14 the evolution as a function of the age of the coefficients $\frac{c_1}{c_1+c_2}$ and $\frac{c_2}{c_1+c_2}$ appearing in Eq. (19) for these two evolutionary sequences. The contribution of k_2, 1 to ${\overline{k}}_2$ increases with time and, assuming an age older than 5 Myr, amounts to minimum 90%. These arguments, together with the fact that the stellar parameters of the primary star are better constrained than those of the secondary star, motivate us to first build stellar evolutions models for the primary star.

Fig. 13. Evolution of k2 and $\dot{\omega }$. Left panel: evolution as a function of stellar age of the internal structure constant (to which the empirical correction Eq. (22) has been applied) for Clés models with initial mass of 19.0 M⊙ and ξ = 1 (k2, 1, green) and initial mass of 7.7 M⊙ and ξ = 0.1 (k2, 2, blue), both models have αov = 0.30 and no turbulent diffusion. The weighted-average mean of the k2-values, computed using Eq. (19), is also depicted (${\overline{k}}_2$, purple). Right panel: evolution as a function of stellar age of the total apsidal motion rate (purple) and of the contributions of the two stars (green and blue). The observational value of the corresponding parameter and its error bars are represented by the solid red line and the dashed red horizontal lines, respectively.

Fig. 14. Evolution as a function of stellar age of the contributions $\frac{c_1}{c_1+c_2}$ and $\frac{c_2}{c_1+c_2}$ to the weighted-average mean ${\overline{k}}_2$ of the internal structure constants k2, 1 and k2, 2 (to which the empirical correction Eq. (22) has been applied) for the models presented in Fig. 13.

7.2. Primary star

We built stellar evolution models for the primary star assuming the three constraints M, R, and T_eff, and adopting at least the initial mass of the star and its age as free parameters. The first model (Model I) is built fixing ξ = 1 and α_ov = 0.30, and assuming no turbulent diffusion. The only free parameters are thus the initial mass and the current age of the star. This model gives M_init = 19.05 M_⊙ and an age of 7.30 Myr, but has a slightly smaller effective temperature than observed, even though still within the errorbars (see Table 7). We note that the χ² is computed based on M, R, T_eff, and k_2, 1. Adopting α_ov = 0.0 gives even worse results in terms of M, R, and T_eff. Without any mixing, the global parameters of the star cannot be reproduced accurately.

Hence, we subsequently investigated the influence of ξ, α_ov, and D_T on the best-fit models¹¹. In three series of models, Series II to IV, we fixed ξ to 1, 0.5, and 0.25, respectively. For each series, we built five models with different values of α_ov, namely 0.20, 0.25, 0.30, 0.35, and 0.40, while leaving D_T as an additional free parameter. The results are summarised in Table 7. Series II, III, and IV models have an initial mass of 19.00, 18.82, and 18.73 M_⊙, respectively, and an age between 7.81 and 7.88 Myr, 7.92 and 7.99 Myr, and 7.98 and 8.05 Myr, respectively. The age is lower for larger values of α_ov and ξ. All these fifteen models correctly reproduce the three constraints M, R, and T_eff, meaning that models of equal quality, as far as the adjustment of the three stellar parameters is concerned, can be obtained for different triplets of D_T, ξ, and α_ov. This D_T − ξ − α_ov degeneracy (already highlighted for HD 152248 in Rosu et al. 2020b) is illustrated in Fig. 15. First, for a given value of ξ, the higher the value of α_ov, the lower the value of D_T. These two parameters affect the stellar structure and evolution the same way. They increase the mixing inside the stars, hence increasing the amount of hydrogen transported to the core. As a consequence, if one of these two parameters is increased, the other one decreases accordingly, and vice-versa. Second, for a given value of α_ov, the higher the value of ξ, the lower the value of D_T. This is expected as a higher value of ξ implies a higher mass-loss rate and thus a higher initial mass. When the star reaches the state where its mass equals the observational mass, its convective core is larger than for a star of lower ξ. The larger convective core in turn implies stronger mixing, thereby reducing the need for turbulent mixing.

Fig. 15. Degeneracy between the various processes in the stellar interior for the best-fit min-Clés models: turbulent diffusion coefficient DT as a function of the mass-loss rate scaling parameter ξ for different values of the overshooting parameter αov. The colours stand for αov = 0.20 (green), 0.25 (blue), 0.30 (purple), 0.35 (pink), and 0.40 (orange).

The biggest differences between Series II, III, and IV lie in the different initial masses and ages of the models. For the reasons aforementioned, the higher the ξ value, the higher the initial mass of the stars and, hence, the lower its age when the star reaches the state defined by the observational parameters. Finally, all these models have a k₂ value ∼ two times higher than the ${\overline{k}}_2$ value. There is no significant difference in the k₂ value for models having the same α_ov but a different ξ. However, for the models having the same ξ, k₂ decreases when α_ov increases notably due to a more pronounced density variation at a radius close to the junction between the convective core and the radiative envelope in the overshooting region.

The degeneracy between the best-fit models is illustrated in the Hertzsprung-Russell diagram in Fig. 16 where four evolutionary sequences are presented, together with Model I for comparison and the observational box defined by the observational radius and effective temperature and their respective errors. The evolutionary sequences are built based on the initial mass, ξ, α_ov, and D_T adopted from the Series II(1), II(5), III(1), and IV(1). The sequences corresponding to Series II(1) and II(5) are very similar, while the sequences corresponding to Series III(1) and IV(1) follow slightly different tracks. Nonetheless, all four sequences cross the observational box at the same point and follow very close paths near this point. The five dots over-plotted on each track represent the models having $k_{2,1}={\overline{k}}_2$, which is, as already stated in Sect. 6, the maximum value k_2, 1 can take. These five models lie well beyond the observational box, confirming that the best-fit models are too homogeneous. The filled triangles over-plotted on these tracks represent the models having k_2, 1 = 0.00139 (see Sect. 7.3 for the justification of this value). As expected, these models lie even further away from the observational box than those with $k_{2,1}={\overline{k}}_2$.

Fig. 16. Hertzsprung-Russell diagram: evolutionary tracks of Clés models corresponding to the following best-fit models: Model I (coral), Series II(1) (green), Series II(5) (light blue), Series III(1) (dark blue), and Series IV(1) (plum). The dots and triangles over-plotted on the corresponding tracks correspond to the models for which k2, 1 is equal to ${\overline{k}}_2$ and 0.00139, respectively, while the open circles correspond to the models for which k2, 1 is equal to ${\overline{k}}_2\pm \mathrm{\Delta }{\overline{k}}_2$. The observational value is represented by the red point, and its error bars are represented by the dark red parallelogram.

All series predict mass-loss rates lower than observed. Since the observational Ṁ is poorly constrained, one cannot totally rule out the possibility of a ξ-value different from 1. Yet, without any additional information, we decided, from now on, to restrict our analysis to models having ξ = 1.

We conclude that in order to find a model which reproduces M, R, and T_eff within their error bars, and which also satisfies the constraint that $k_{2,1}<{\overline{k}}_2$ within the error bars of ${\overline{k}}_2$, the value of D_T should be increased for a given value of α_ov. We computed five models (Models V to IX) with the value of α_ov fixed to 0.20, 0.25, 0.30, 0.35, and 0.40, respectively, leaving the age, initial mass, and D_T as free parameters. Compared to previously computed models, one additional constraint is enforced: $k_{2,1}={\overline{k}}_2$ within the error bars. The resulting models are reported in Table 7.

The masses and radii of all these models are comprised within the error bars on these parameters. However, none of these models has an effective temperature and k₂-value compatible with observations: both parameters are overestimated. This highlights the difficulty to reconcile the observational and theoretical k₂-values, especially since, as already discussed in Sect. 6, ${\overline{k}}_2$ is only the upper-bound value on k_2, 1. Requesting k_2, 1 to be much lower than ${\overline{k}}_2$ would lead to an even more pronounced discrepancy.

Even though the mass and radius can be correctly reproduced within the error bars, it seems impossible to find a model reproducing both the effective temperature and the internal structure constant simultaneously. This issue is highlighted in the Hertzsprung-Russell diagram in Fig. 17: Three evolutionary sequences are presented, together with the observational box. The evolutionary sequences are built based on the initial mass, ξ, α_ov, and D_T adopted from the Series II(1), Model V, and Model IX. As expected, the tracks corresponding to Models V and IX do not cross the observational box. The three dots over-plotted on each track represent the models having $k_{2,1}={\overline{k}}_2$, while the filled triangles over-plotted on these tracks represent the models having k_2, 1 = 0.00139 (see Sect. 7.3 for the justification of this value). These six models lie well beyond the observational box, confirming that the best-fit models still have a too small density contrast between the core and the external layers. The two stars over-plotted on the tracks of Models V and IX correspond to the best-fit models given in Table 7 (we note that the one corresponding to Series II(1) perfectly matches the observational value).

Fig. 17. Hertzsprung-Russell diagram: evolutionary tracks of Clés models corresponding to the following best-fit models: Series II(1) (green), Model V (purple), and Model IX (pink). The dots and triangles over-plotted on the corresponding tracks correspond to the models for which k2, 1 is equal to ${\overline{k}}_2$ and 0.00139, respectively, while the open circles correspond to the models for which k2, 1 is equal to ${\overline{k}}_2\pm \mathrm{\Delta }{\overline{k}}_2$. The stars on the purple and pink tracks correspond to the best-fit models Model V and Model IX, respectively. The observational value is represented by the red point, and its error bars are represented by the dark red parallelogram.

7.3. Secondary star

We built stellar evolution models for the secondary star assuming the three constraints M, R, and T_eff and adopting the initial mass, the age, and the turbulent diffusion as free parameters. We built seven models (Models I_S to VII_S) fixing ξ = 0.1 and adopting values for α_ov of 0.10, 0.15, 0.20, 0.25, 0.30, 0.35, and 0.40. The resulting best-fit models are summarised in Table 7. All models reproduce perfectly the mass and radius, but slightly overestimate the effective temperature even though still within the error bars. All these models have an initial mass of 7.70 M_⊙ and an age ranging between 10.9 and 12.5 Myr. We observe that D_T converges towards a value of 0 cm² s⁻¹ for all models. Assuming Models IX and I_S correctly reproduce the primary and secondary stars, respectively, we get an apsidal motion rate of 1.8793°  yr⁻¹ for the binary system. This value is 57% higher than the observational value and is therefore incompatible with the observational results. The second problem here lies in the fact that the age of these models is incompatible with the range of ages derived for the primary star in Sect. 7.2. The internal structure constant of these models, corrected for the rotation of the star through Eq. (22), amounts to ∼8 × 10⁻³. Assuming these models correctly reproduce the secondary star, then $k_{2,2}\sim 4.7{\overline{k}}_2$ and, making use of Eq. (19), we deduce that k_2, 1 should take the value of 0.00139. If it appeared that an enhanced mixing would be necessary to model the secondary star, then on the one hand k_2, 2 would be lower and k_2, 1 higher but, on the other hand, the age of the secondary star would be higher, thus reenforcing the age discrepancy between the primary and secondary stars.

7.4. Towards a common age of the stars

In order to examine the age discrepancy, we built series of evolutionary tracks with Clés for the two stars. We assumed α_ov = 0.30 for both stars. For the primary star, we assumed ξ = 1 and built 20 evolutionary tracks with initial mass ranging from 17 to 21 M_⊙ with a step of 1 M_⊙ and a D_T value of 0, 2 × 10⁶, 6 × 10⁶, and 10 × 10⁶ cm² s⁻¹. For the secondary star, we assumed ξ = 0.1 and built 9 evolutionary tracks with initial mass of 7.5, 7.7, and 8.0 M_⊙, and a D_T value of 0, 2 × 10⁵, and 5 × 10⁵ cm² s⁻¹.

Evolutionary tracks for the primary star are presented in the Hertzsprung-Russell diagram in Fig. 18. In this figure, only the evolutionary tracks crossing (or nearly crossing) the observational box are presented. The models with an initial mass of 17 and 18 M_⊙ have a too low mass compared to observational, even within the error bars, and therefore do not fulfil the constraint on the mass. On the contrary, the models with an initial mass of 20 and 21 M_⊙ have too high a value of the mass compared to the observational one and, at best, fulfil the constraint on the mass at the very end of the main-sequence phase. Only the evolutionary tracks having an initial mass of 19 M_⊙ have a mass compatible with the observational value during the major part of the main-sequence phase – or at least when crossing (or nearly crossing) the observational box. During the main-sequence phase, the turbulent diffusion impacts the evolutionary tracks as follows: the higher the turbulent diffusion in absolute value, the higher the bolometric luminosity of the star for a given effective temperature. Hence, for low-mass models (i.e. models with an initial mass of 17 and 18 M_⊙), only models with enhanced turbulent mixing cross the observational box, and only those having a D_T-value of 6 × 10⁶ and 10 × 10⁶ cm² s⁻¹ have a k₂-value (almost) compatible with ${\overline{k}}_2$ within the error bars when crossing the observational box. On the opposite, for the high-mass models (i.e. models with an initial mass of 20 and 21 M_⊙), only the models without or with a small D_T-value cross the observational box, but none of them has a k₂-value compatible with ${\overline{k}}_2$ within the error bars when crossing the observational box. Regarding the models with an initial mass of 19 M_⊙, only those with D_T equal to 0, 2 × 10⁶, and 6 × 10⁶ cm² s⁻¹ cross the observational box, but none of them has a k₂-value compatible with ${\overline{k}}_2$, even though the track with D_T = 6 × 10⁶ cm² s⁻¹ has a k₂-value very close to this value. These evolutionary sequences highlight the issue encountered in Sect. 7.2 about the difficulty to reconcile the mass, the radius, the effective temperature, and the internal structure constant of the star simultaneously.

Fig. 18. Hertzsprung-Russell diagram: evolutionary tracks of Clés models for the primary star of Minit = 17 M⊙ (green), 18 M⊙ (light blue), 19 M⊙ (dark blue), 20 M⊙ (pink), and 21 M⊙ (orange), and DT = 0 cm2 s−1 (solid line), 2 × 106 cm2 s−1 (dashed line), 6 × 106 cm2 s−1 (dotted-dashed line), and 10 × 106 cm2 s−1 (dotted line). All models have αov = 0.30 and ξ = 1. The dots over-plotted on the corresponding tracks correspond to the models for which k2, 1 is equal to ${\overline{k}}_2$, while the open circles correspond to the models for which k2, 1 is equal to ${\overline{k}}_2\pm \mathrm{\Delta }{\overline{k}}_2$. The filled and open squares over-plotted on the tracks having an initial mass of 19 M⊙ correspond to the models for which the mass equals the observational value within the error bars. The observational value is represented by the red point, and its error bars are represented by the dark red parallelogram.

Evolutionary tracks for the secondary star are presented in the Hertzsprung-Russell diagram in Fig. 19. All these tracks cross the observational box defined by the radius and effective temperature at some point of their evolution. However, only the models with an initial mass of 7.7 M_⊙ have a mass compatible with the observation when crossing the observational box. Given that the mass-loss rate is negligible, only evolutionary tracks having an initial mass ranging between 7.58 and 7.82 M_⊙ are acceptable. Regarding the evolutionary track of initial mass 7.7 M_⊙ and D_T = 0 cm² s⁻¹, the model closest to the observational value has a k_2, 2 value approximately equal to 4.7 ${\overline{k}}_2$, confirming the value we found in Sect. 7.3. Within the observational box, the age of the star ranges from 10 to 14 Myr, above the age of the best-fit models of the primary star of approximately 9 Myr. The filled 5-branch stars over-plotted on the tracks correspond to the models having an age of 9 Myr. These models are slightly below the observational box meaning that their stellar properties are not compatible with the observational values.

Fig. 19. Hertzsprung-Russell diagram: evolutionary tracks of Clés models for the secondary star of Minit = 7.5 M⊙ (water green), 7.7 M⊙ (dark blue), and 8.0 M⊙ (orange), and DT = 0 cm2 s−1 (solid line), 2 × 105 cm2 s−1 (dashed line), and 5 × 105 cm2 s−1 (dotted-dashed line). All models have αov = 0.30 and ξ = 0.1. The filled 5-branch stars over-plotted on the tracks correspond to the models for which the age is equal to the age of the best-fit models of the primary star (i.e. about 9 Myr). The observational value is represented by the red point, and its error bars are represented by the dark red parallelogram.

The influence of a potential non-solar metallicity and helium abundance is investigated in Appendix B. Besides being highly unlikely, we showed that a non-solar metallicity and helium abundance does not solve the discrepancies discussed above.

We inevitably come to the conclusion that we cannot find stellar evolution models for the two stars reproducing the stellar properties properly and having a compatible age at the same time. Whilst stellar evolution models have their intrinsic limitations, notably regarding the way the internal mixing is implemented, we might also have been biased by the observational data. There is no obvious issue arising from the spectroscopic and radial velocity analysis and for example the primary effective temperature and the mass ratio that we have derived should be robust. However, the situation is more complicated for photometry. The zero contribution of the third light to the photometric data of the binary system might have been underestimated. The Roche lobe filling factors of the two stars and the effective temperature of the secondary star have been determined based on the assumed value for the third light contribution. These parameters, and consequently also the radii of both stars, are not as reliable as their error bars might suggest. Yet, by taking I₃ = 0, we have actually the smallest possible radius of the secondary star, and larger radii would lead to even more severe age discrepancies than what we have found here.

8. Discussion: possible misalignment of the rotation axis

The highlighted difficulty to reconcile observational and theoretical internal structure constant values suggests there might be some bias in the interpretation of the apsidal motion in terms of the internal structure constant. One possibility could be that higher-order terms in the expression of the apsidal motion rate have been neglected. Yet, these higher-order terms have been shown to have a negligible contribution (Rosu et al. 2020b, and references therein). Another scenario could arise from the presence of a third body orbiting the binary system. However, to date, there is no signature of a third body in the spectrum of HD 152219 (see review in Sect. 1 as well as Sects. 3 and 5). The possibility that a ternary star would contribute to the apsidal motion rate without being detected though is a priori conceivable (see, e.g., Naoz et al. 2013 and Borkovits et al. 2015 for a review of the ensuing apsidal motion equations). At this stage, there is unfortunately no way to infirm or confirm the presence of a third body in the system. Though, it is highly unlikely that a non-detected third body would contribute for as much as 57% of the inner binary apsidal motion rate necessary to reconcile theoretical and observational k₂-values (see Sect. 7.3).

A third possibility could be a misalignment of the rotation axes of the stars with respect to the normal to the orbital plane. Indeed, the Newtonian contribution to the total apsidal motion rate given by Eq. (13) assumes the rotation axes of the stars to be parallel to the normal to the orbital plane. In case of a misalignment, the Newtonian contribution takes the general expression given by Eq. (3) of Shakura (1985):

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\dot{\omega }}_{\text{N}}=\frac{2\pi }{P_{\text{orb}}}[& 15f(e)\{k_{2,1}q{\left(\frac{R_1}{a}\right)}^5+\frac{k_{2,2}}{q}{\left(\frac{R_2}{a}\right)}^5\}\hfill \\ \hfill & \begin{array}{cc}\hfill -\frac{g(e)}{{(sini)}^2}\{& k_{2,1}(1+q){\left(\frac{R_1}{a}\right)}^5{\left(\frac{P_{\text{orb}}}{P_{\text{rot,1}}}\right)}^2\hfill \\ \hfill & [cos{\alpha }_1(cos{\alpha }_1-cos{\beta }_{1}cosi)\hfill \\ \hfill & +\frac{1}{2}{(sini)}^2(1-5{(cos{\alpha }_1)}^2)]\hfill \\ \hfill & +k_{2,2}\frac{1+q}{q}{\left(\frac{R_2}{a}\right)}^5{\left(\frac{P_{\text{orb}}}{P_{\text{rot,2}}}\right)}^2\hfill \\ \hfill & [cos{\alpha }_2(cos{\alpha }_2-cos{\beta }_{2}cosi)\hfill \\ \hfill & +\frac{1}{2}{(sini)}^2(1-5{(cos{\alpha }_2)}^2)]\}].\hfill \end{array}\hfill \end{array}\end{array}$$(24)

In this expression, α₁ (resp. α₂) stands for the angle between the primary (resp. secondary) star rotation axis and the normal to the orbital plane, and β₁ (resp. β₂) stands for the angle between the primary (resp. secondary) star rotation axis and the line joining the binary centre and the observer.

We define the function F_α, j (with j = 1, 2 for the primary and secondary stars, respectively) as

$$\begin{array}{c}\hfill F_{\alpha ,j}=cos{\alpha }_j(cos{\alpha }_j-cos{\beta }_{j}cosi)+\frac{1}{2}{sin}^{2}i(1-5{cos}^2{\alpha }_j).\end{array}$$(25)

We can express cos β_j in terms of the inclination and α_j using the cosine relationship in spherical trigonometry:

$$\begin{array}{c}\hfill cos{\beta }_j=cosicos{\alpha }_j+sinisin{\alpha }_{j}cos\theta ,\end{array}$$(26)

where θ is the azimutal angle of the rotation axes of the stars, such that θ = 0° and 180° correspond to the cases where the rotation axes, the line of sight, and the normal to the orbital plane lie in the same plane. Replacing in Eq. (25) cos β_j by its expression given by Eq. (26) gives us an expression of F_α, j which depends upon the inclination, the α_j, and θ only. We can further show that F_α, j takes its extremum values for θ = 0° and 180°. Assuming that θ = 180°, we get β_j = i + α_j.

We define the quantities

$$\begin{array}{c}\hfill C_{\text{mis},j}=-\frac{1}{{(sini)}^2}(1+\frac{M_j}{M_{3-j}}){\left(\frac{R_j}{a}\right)}^5{\left(\frac{P_{\text{orb}}}{P_{\text{rot},j}}\right)}^2{F}_{\alpha ,j},\end{array}$$(27)

and represent these as a function of the angle α_j in Fig. 20 for the specific case of i = 89.58°. C_mis, j takes its minimum and maximum values at α_j = 90° and 180°, respectively.

Fig. 20. Behaviour of Cmis, 1 and Cmis, 2 as a function of α1 and α2, the angle between the primary and secondary stellar rotation axis, respectively, and the normal to the orbital plane.

We computed ${\overline{k}}_{2,\text{mis}}$ considering values of α_j ranging from 0 to 180°. Figure 21 shows the ${\overline{k}}_{2,\text{mis}}$ value as a function of the two angles α₁ and α₂. The maximum value of ${\overline{k}}_{2,\text{mis}}$ is 3.39 × 10⁻³ and is reached for α₁ = α₂ = 90°. As expected, the ${\overline{k}}_{2,\text{mis}}$-value mostly changes with the primary stellar rotation axis inclination α₁ and is only slightly affected by a change in the secondary stellar rotation axis inclination. This reflects the predominant contribution of the primary star to the total rate of apsidal motion. For α₁ lower than about 40°, there is no significant impact on the ${\overline{k}}_{2,\text{mis}}$-value.

Fig. 21. Behaviour of ${\overline{k}}_{2,\text{mis}}$ as a function of α1 and α2, the angle between the primary and secondary stellar rotation axis, respectively, and the normal to the orbital plane.

In Sect. 7.2, the best-fit models we obtained for the primary star all had a k_2, 1-value between ∼2.4 and 2.5 × 10⁻³. For these values to be compatible with the observational value, the primary stellar rotation axis should be inclined by an angle of at least ∼50° with respect to the normal to the orbital plane. Such a high misalignment angle would be surprising in a close binary system, especially since the two stars are about to reach pseudo-synchronisation.

An important question is whether the condition of sub-critical stellar rotation leads to a restriction on the values of α₁ and α₂. The condition of sub-critical rotation velocity expresses the fact that the stellar rotation rate Ω_rot cannot exceed the critical value Ω_crit corresponding to the centrifugal force at the stellar equator exactly compensating the effective gravity (i.e. the gravity corrected for the effect of radiation pressure). This condition writes

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{rot}}\le {\mathrm{\Omega }}_{\mathrm{crit}}=\sqrt{\frac{G{M}_{\ast }}{R_{\ast }^3}(1-{\mathrm{\Gamma }}_e)},\end{array}$$(28)

where Γ_e is the Eddington factor:

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_e=\frac{{\kappa }_e{L}_{\ast }}{4\pi cG{M}_{\ast }},\end{array}$$(29)

with κ_e being the electron scattering opacity. For this parameter, we adopted the value κ_e ≃ 0.34 cm² g⁻¹ (e.g., Sanyal et al. 2015). Hence, we get ${\mathrm{\Gamma }}_e=2.6{10}^{-5}\frac{L_{\ast }{M}_{\odot }}{L_{\odot }{M}_{\ast }}$. In terms of the equatorial rotational velocity, this translates into

$$\begin{array}{c}\hfill v_{\mathrm{eq}}\le 436.7\sqrt{\frac{M_{\ast }{R}_{\odot }}{M_{\odot }{R}_{\ast }}(1-2.6{10}^{-5}\frac{L_{\ast }{M}_{\odot }}{L_{\odot }{M}_{\ast }})}{\mathrm{km}\mathrm{s}}^{-1}.\end{array}$$(30)

For the primary and secondary stars, we get v_eq, 1 ≤ 583 km s⁻¹ and v_eq, 2 ≤ 628 km s⁻¹, respectively. Using the projected rotational velocities derived in Sect. 3.3, namely v_eq, 1 sin β₁ = 166 km s⁻¹ and v_eq, 2 sin β₂ = 95 km s⁻¹, we get the following constraints on the angles β_j: β₁ ≥ 16.5° and β₂ ≥ 8.7°. In order to translate this condition on β_j into a condition on α_j, we need to assume a value for the azimuthal angle θ. As shown before, the largest impact of α_j on the apsidal motion rate occurs when θ equals 0° and 180°. For an azimuthal angle of 0°, we get the conditions that α₁ > 106.1° or α₁ < 73.0°, and α₂ > 98.3° or α₂ < 80.9°, whilst for an azimuthal angle of 180°, we get the conditions that α₁ > 107.0° or α₁ < 73.9°, and α₂ > 99.1° or α₂ < 81.7°. These conditions do not rule out the possibility to have inclination angles of ∼50° for the stellar rotation axis of the stars.

The possibility of a misalignment between the primary stellar rotation axis and the axis of the binary system can in principle be tested via the Rossiter-McLaughlin effect (RM, e.g., Hosokawa 1953; Ohta et al. 2005; Giménez 2006; Albrecht et al. 2007) during the eclipse of the primary star by the secondary. As the secondary passes in front of the primary, assuming a prograde orbit, it first eclipses those parts of the rotating primary that produce the blue wings of the spectral lines, thereby shifting the line centroids to the red. Likewise, near the end of the eclipse, the secondary occults those parts of the primary that produce the red wings, resulting in a shift of the line centroids to the blue. The signature of a RM effect in HD 152219 was reported by Sana (2009) from measurements of the He IIλ 4686 line which, in this system, is only formed in the spectrum of the primary star.

We have designed a Fortran code to simulate the RM effect for different values of the α₁ and β₁ angles, adopting the photometric value of the orbital inclination and the primary’s v_eq, 1 sin β₁. Synthetic line profiles were computed adopting the quadratic limb darkening law of Claret (2000) and assuming the stars to have spherical shapes with radii as determined from our photometric solution. The centroid of the synthetic profiles was computed by means of the first order moment of the line profile.

Sana (2009) reported maximum RV deviations of about 26 km s⁻¹, with the negative and positive RV deviations having nearly identical amplitudes. We found the amplitudes predicted by our code difficult to compare with those reported by Sana (2009). Indeed, Sana (2009) used Gaussian profile fitting to determine the line RVs, whilst we used the first order moment instead. Our method resulted in an amplitude of the RM effect of at most 19 km s⁻¹. Because of the highly non-gaussian shape of the line profiles during the eclipse, the Gaussian fitting results in a different amplitude of the RM effect, which is difficult to reproduce as it also depends on the wavelength range over which the Gaussian fitting is performed. Nevertheless, the fact that the negative and positive RV deviations have nearly identical amplitudes yields a constraint on the value of α: α ≤ 15° and α ≥ 165°.

Moreover, given the proximity of the system’s orbital inclination to 90°, the largest amplitude of the RM effect is expected for an alignment of the rotation and binary axes. In view of the amplitudes reported by Sana (2009), and keeping in mind the above caveat about these amplitudes, it seems thus unlikely that there exists a very large misalignment between the rotation and binary axes. Whilst this conclusion applies of course only to the primary star, we stress that it is actually the primary that has the largest rotational contribution to the rate of apsidal motion. Hence, a misalignment of the axes seems unlikely to be responsible for the discrepancy between the observational and theoretical values of $\dot{\omega }$.

9. Conclusion

The eccentric massive binary HD 152219, belonging to the open cluster NGC 6231, has been analysed both from an observational and a stellar evolution point of views.

We made use of the huge set of spectroscopic observations of the system to reconstruct the individual spectra of the components by means of a disentangling code. These spectra were subsequently analysed by means of the CMFGEN model atmosphere code to asses the properties of the binary components, notably an effective temperature of 30 900 ± 1000 K for the primary star. The radial velocities obtained via disentangling were combined with those coming from the literature to establish the rate of apsidal motion along with the orbital period and eccentricity. For these parameters, we respectively found $\dot{\omega }=(1.198\pm 0.300)°$ yr⁻¹, $P_{\text{orb}}=4.{24046}_{-0.00004}^{+0.00005}$ d, and $e=0.{072}_{-0.005}^{+0.004}$. The TESS-12 light curve was analysed by means of the Nightfall code and we notably inferred an effective temperature of 21 697 ± 1000 K for the secondary star and an inclination of $(89.{58}_{-2.28}^{+0.42})°$ for the system. The Roche lobe filling factors allowed us to infer stellar radii of $9.{40}_{-0.15}^{+0.14}$ R_⊙ and 3.69 ± 0.06 R_⊙ for the primary and secondary stars, respectively. Combining the results from the radial velocity analysis and light curve analysis, we derived masses of 18.64 ± 0.47 M_⊙ and 7.70 ± 0.12 M_⊙ for the primary and secondary stars, respectively. We inferred a value of 0.00173 ± 0.00052 for the weighted-average mean of the internal structure constant of the binary.

We built Clés stellar evolution models using the min-Clés routine to search for best-fit models of the two stars adopting the mass, radius, and effective temperature as constraints, sometimes complemented by k₂ for the primary star. A simultaneous analysis of the two stars gave no coherent results. We considered different prescriptions for the internal mixing occurring inside the stars. We notably performed several tests with different values of the overshooting parameter α_ov, the mass-loss rate scaling factor ξ, and the turbulent diffusion D_T. The primary star models were shown to have a smaller density contrast than suggested by the observations, that is to say, the theoretical internal structure constant was two times greater than the observational weighted-mean ${\overline{k}}_2$-value. To simultaneously satisfy the three previously mentioned constraints together with the constraint that $k_{2,1}<{\overline{k}}_2$ in the models was impossible: The models were shown to have a strongly enhanced turbulent diffusion but both the effective temperature and k₂-value were overestimated. The best-fit models, that is to say the models having the best compromise in terms of M, R, T_eff, and k_2, 1, all have an initial mass of 18.3 ± 0.5 M_⊙ and an age of 9.5 ± 0.6 Myr. Regarding the secondary star, given that there is no strong constraint on k_2, 2, we only enforced the stellar models to reproduce M, R, and T_eff. The best-fit models obtained in this way all have an initial mass of 7.70 ± 0.12 M_⊙, a k_2, 2-value approximately equal to 4.7 ${\overline{k}}_2$, and an age estimate of 11.7 ± 0.6 Myr. Assuming these stellar models are representative of the secondary star, the k_2, 2-value implies that k_2, 1 = 0.00139, thereby increasing even more the discrepancy with the stellar evolution models of the primary star. All things considered, the age estimated for the secondary star is incompatible with the age of the primary star and, furthermore, with the age estimate of the cluster. We also tested the impact of non-solar metallicity and helium abundance, and found that both an enhanced metallicity and depleted helium abundance go towards adjustments of the primary star that better reproduce M, R, T_eff, and k_2, 1. However, the corresponding models of the secondary star yield even higher ages, further increasing the inconsistency between primary and secondary ages.

Finally, we investigated the impact of a misalignment of the stellar rotation axis with respect to the normal to the orbital plane. Such a misalignment reduces the contribution of the rotational term to the Newtonian apsidal motion rate, thereby implying an increased value of ${\overline{k}}_2$ compared to the aligned configuration. In the present case, we found that a severe misalignment angle is required to get a ${\overline{k}}_2$-value high enough for the stellar models of the primary star to satisfy the constraint on k_2, 1. Such a high misalignment is however ruled out by the observational constraint obtained via the analysis of the Rossiter-McLaughlin effect which sets an upper limit of 15° on the misalignment angle.

The difficulty to reconcile the k₂-value of the primary stellar models with the observational constraint on the ${\overline{k}}_2$-value as well as the internal inconsistency between the ages of the two stars are indications that some physics of the stellar interior are still not completely understood. Obviously, all the tools used in our observational or theoretical analysis have their own limitations. Whilst we are quite confident in our determination of the masses of the stars as well as in the effective temperature of the primary star, we are less confident in the inferred values of the secondary effective temperature and the radii of both stars.

Acknowledgments

S.R. acknowledges support from the Fonds de la Recherche Scientifique (F.R.S.- FNRS, Belgium). M.F. acknowledges support from STFC consolidated grant ST/T000252/1. We thank Dr John Hillier for making his code CMFGEN publicly available. We thank Dr Rainer Wichmann for making his code Nightfall publicly available as well as for discussions regarding his code. We also thank Dr Eric Gosset for our numerous discussions about the χ². We finally thank Dr Yaël Nazé for her contribution in the extraction of the TESS light curves. This paper makes use of data collected by the TESS mission, whose funding is provided by the NASA Explorer Program.

References

Appendix A: Journal of the spectroscopic observations

This appendix provides the journal of the spectroscopic observations of HD 152219 including the radial velocities obtained by means of the disentangling code (Table A.1).

Appendix B: Influence of the metallicity and helium abundance in stellar evolution models

As a last attempt to reconcile the ages of the primary and secondary stars, we investigated the influence of a potential non-solar metallicity and helium abundance. Both the metallicity and the initial helium abundance affect the internal structure essentially through a change in the opacity. The uncertainties on the opacity hence justify the influence analysis we performed here.

Regarding the metallicity, we built eight evolutionary sequences for the primary star, all having an initial mass of 19 M_⊙, α_ov = 0.30, and ξ = 1, but different values of Z, namely 0.010 and 0.020, and different values of D_T, namely 0, 2 × 10⁶, 6 × 10⁶, and 10 × 10⁶ cm² s⁻¹. These evolutionary tracks are presented in the Hertzsprung-Russell diagram in Fig. B.1 together with the evolutionary tracks for Z = 0.015 presented in Sect. 7.4. For more clarity, only the tracks crossing at some point of their evolution the observational box are shown. For Z = 0.010, only the two models with the lowest values of D_T cross the observational box, while for Z = 0.020, all models cross the observational box. In addition, all tracks have a mass compatible with the observational value when crossing the observational box. Compared to Fig. 18, we observe that a higher metallicity has the same impact on the evolutionary tracks as a lower mass and/or lower D_T: The evolutionary track with Z = 0.020 and the highest value of D_T has now a value of k_2, 1 compatible with ${\overline{k}}_2$ when it crosses the observational box.

Fig. B.1. Hertzsprung-Russell diagram: evolutionary tracks of Clés models for the primary star having Minit = 19 M⊙, αov = 0.30, ξ = 1, X = 0.715, Z = 0.010 (green), 0.015 (dark blue), and 0.020 (orange), and DT = 0 cm2 s−1 (solid line), 2 × 106 cm2 s−1 (dashed line), 6 × 106 cm2 s−1 (dotted-dashed line), and 10 × 106 cm2 s−1 (dotted line). The dots over-plotted on the corresponding tracks correspond to the models for which k2, 1 is equal to ${\overline{k}}_2$, while the open circles correspond to the models for which k2, 1 is equal to ${\overline{k}}_2\pm \mathrm{\Delta }{\overline{k}}_2$. The filled and open squares over-plotted on the tracks correspond to the models for which the mass is equal to the observed value within the error bars. The observational value is represented by the red point, and its error bars are represented by the dark red parallelogram.

Regarding the helium abundance, we built eight evolutionary sequences for the primary star, all having an initial mass of 19 M_⊙, α_ov = 0.30, ξ = 1, and Z = 0.015, but different values of X, namely 0.700 and 0.730, and different values of D_T, namely 0, 2 × 10⁶, 6 × 10⁶, and 10 × 10⁶ cm² s⁻¹. These evolutionary tracks are presented in the Hertzsprung-Russell diagram in Fig. B.2 together with the evolutionary tracks having X = 0.715 presented in Sect. 7.4. For more clarity, only the tracks crossing at some point of their evolution the observational box defined by the radius and effective temperature are plotted. For all values of X, only the three models with the lowest values of D_T cross the observational box at some point of their evolution. In addition, all tracks have a mass compatible with the observed value when crossing the observational box. Compared to Fig. 18, we observe that a higher X-value, hence a lower helium abundance, has the same impact on the evolutionary tracks as a lower mass and/or lower D_T. The impact of a change in X is however smaller than that of a change in Z and none of the models has a k_2, 1 compatible with ${\overline{k}}_2$ when crossing the observational box.

Fig. B.2. Hertzsprung-Russell diagram: evolutionary tracks of Clés models for the primary star having Minit = 19 M⊙, αov = 0.30, ξ = 1, Z = 0.015, X = 0.700 (water green), 0.715 (dark blue), and 0.730 (pink), and DT = 0 cm2 s−1 (solid line), 2 × 106 cm2 s−1 (dashed line), 6 × 106 cm2 s−1 (dotted-dashed line), and 10 × 106 cm2 s−1 (dotted line). The dots over-plotted on the corresponding tracks correspond to the models for which k2, 1 is equal to ${\overline{k}}_2$, while the open circles correspond to the models for which k2, 1 is equal to ${\overline{k}}_2\pm \mathrm{\Delta }{\overline{k}}_2$. The filled and open squares over-plotted on the tracks correspond to the models for which the mass is equal to the observed value within the error bars. The observational value is represented by the red point, and its error bars are represented by the dark red parallelogram.

We then built three best-fit models for the primary star and the secondary star adopting the constraints on M, R, and T_eff, fixing α_ov = 0.30 for both stars, and ξ = 1 for the primary star and 0.1 for the secondary star, and leaving the age, initial mass, and D_T as free parameters, for the three following couples of values for (X, Z): (0.715, 0.020), (0.730, 0.015), and (0.730, 0.020). The results are given in Table B.1 under the names of Model Z+, X+, and Z+X+, respectively. All six models correctly reproduce the observational mass, radius, and effective temperature of the corresponding star. Regarding the three primary star models, compared to Series II(3), all models have a higher age, a slightly higher initial mass, and a higher turbulent diffusion coefficient D_T. Their lower k_2, 1-value compared to Series II(3) make these three models better suited to reproduce the primary star. Regarding the secondary star, compared to Model III_S, the three models have a higher age, the same initial mass, and a higher D_T-value, as well as a lower k_2, 2-value. The evolutionary tracks corresponding to these best-fit models are presented in Fig. B.3 together with Series II(3) for the primary and Model III_S for the secondary. Compared to the models having a solar metallicity and helium abundance, the age discrepancy between the primary and secondary best-fit models is enhanced when higher metallicity or/and lower helium abundance is/are considered.

Fig. B.3. Hertzsprung-Russell diagram: evolutionary tracks of Clés models for the primary star (upper panel) and secondary star (lower panel) corresponding to the best-fit models listed in Table B.1 together with Series II(3) for the primary and Model IIIS for the secondary. All models have αov = 0.30, primary models have ξ = 1 while secondary models have ξ = 0.1. The dots over-plotted on the evolutionary tracks of the primary correspond to the models for which k2, 1 is equal to ${\overline{k}}_2$, while the open circles correspond to the models for which k2, 1 is equal to ${\overline{k}}_2\pm \mathrm{\Delta }{\overline{k}}_2$. The filled 5-branch stars over-plotted on the secondary tracks correspond to the models for which the age is equal to the age of the best-fit model of the primary star of same X and Z. The observational value is represented by the red point, and its error bars are represented by the dark red parallelogram.

The increase in D_T with metallicity can be explained as follows (see e.g. Rosu et al. 2020b). An increase in metallicity, at a given mass, increases the opacity which in turn lowers the effectiveness of the radiative transport, hence decreasing the effective temperature and luminosity. The same holds true for the increase in D_T with the decrease in helium abundance which is the direct consequence of an increase in the opacity with X. Since an increase in D_T leads to an increase in the luminosity, it compensates the decrease in T_eff and L induced by the enhanced metallicity. Given that the turbulent diffusion increases the mass of the convective core, it directly follows that k₂ decreases accordingly. In this search for best-fit models with enhanced metallicity, the modification of the various parameters is not a direct consequence of a change in the metallicity but rather an indirect consequence resulting from the corresponding change in the turbulent diffusion.

There is no observational evidence that HD 152219 would have a non-solar metallicity or helium abundance. Even though a higher metallicity and a subsolar helium abundance would improve the best-fit models of the primary star by allowing us to simultaneously satisfy the constraints on M, R, T_eff, and k_2, 1, such a high metallicity seems rather unlikely, especially since Kilian et al. (1994) found a subsolar CNO mass fraction from the analysis of ten B-type stars in NGC 6231. At this stage, insufficient information is available to conclude that the metallicity and helium abundance of the star deviate from solar.

All Tables

All Figures

	Fig. 1. ASAS-3 light curve phase-folded using the orbital period given in Table 3. The top panel illustrates the initial ASAS-3 data set, whilst the bottom panel shows the data after the 1.25σ clipping.
In the text

	Fig. 2. Fourier periodogram of the TESS data from sector 12. The insert zooms on the low frequency domain. The red lines yield the orbital frequency νorb = 0.2358 d−1 and its first six harmonics.
In the text

	Fig. 3. Top row: C IIλ 4267 line profiles of the separated spectra obtained after application of the brightness ratio for the primary (left panel) and secondary (right panel) stars. Bottom row: Fourier transform of those lines (in black) and best-match rotational profile (in red) for the primary (left panel) and secondary (right panel) stars.
In the text

	Fig. 4. Determination of vblack on the absorption trough of C IVλλ 1548, 1551 as observed in the IUE spectra of HD 152219. The black and red curves correspond to observations SWP 56939 and SWP 56953, respectively.
In the text

	Fig. 5. Normalised disentangled spectra (in black) of the primary and secondary stars of HD 152219 (we note that the spectrum of the secondary star is shifted by +0.3 in the y-axis for convenience) together with the respective best-fit CMFGEN model atmosphere (in red).
In the text

	Fig. 6. Confidence contours for the best-fit parameters obtained from the adjustment of the full set of 134 primary RV data of HD 152219 with Eqs. (5) and (6). The best-fit solution is shown in each panel by the black filled dot. The corresponding 1σ confidence level is shown by the blue contour.
In the text

	Fig. 7. Comparison between the measured RVs of the primary (filled dots) and secondary (open dots, when available) with the RVs expected from relations (5) and (6) with the best-fit parameters given in Table 3.
In the text

	Fig. 8. Upper panel: best-fit Nightfall solution of the TESS-12 light curve of HD 152219. Lower panel: residuals over the best-fit solution.
In the text

	Fig. 9. Confidence contours for the best-fit parameters obtained from the adjustment of the TESS-12 light curve of HD 152219 with the Nightfall code. The best-fit solution (i.e. with a third light contribution I3 = 0) is shown in each panel by the black dot, while the corresponding 1σ confidence level is shown by the blue contour.
In the text

	Fig. 10. Upper panel: best-fit Nightfall solution of the TESS-39 light curve of HD 152219. Lower panel: residuals over the best-fit solution.
In the text

Fig. 11. Values of ω as a function of time inferred from the photometric light curves and the RVs. The pink symbols correspond to the data of the fits of the TESS-12 and TESS-39 photometry. The blue dot indicates the ω0 value obtained from the global fit of all RV data. The solid blue line corresponds to our best-fit value of $\dot{\omega }$ inferred from the RVs, and the hatched cyan zone corresponds to the range of values according to the 1σ uncertainties on ω0 and $\dot{\omega }$.

In the text

Fig. 12. Values of the phase difference Δϕ between the primary and secondary eclipses as a function of time and ω inferred from the photometric light curves and the RVs. The pink filled symbol and the pink open symbols correspond to the data of the fits of the ASAS-3 photometry, and the TESS-12 and TESS-39 photometry, respectively. The solid blue line corresponds to our best-fit value of Δϕ inferred from the RVs, and the hatched cyan zone corresponds to the range of values according to the 1σ uncertainties on e.

In the text

Fig. 13. Evolution of k2 and $\dot{\omega }$. Left panel: evolution as a function of stellar age of the internal structure constant (to which the empirical correction Eq. (22) has been applied) for Clés models with initial mass of 19.0 M⊙ and ξ = 1 (k2, 1, green) and initial mass of 7.7 M⊙ and ξ = 0.1 (k2, 2, blue), both models have αov = 0.30 and no turbulent diffusion. The weighted-average mean of the k2-values, computed using Eq. (19), is also depicted (${\overline{k}}_2$, purple). Right panel: evolution as a function of stellar age of the total apsidal motion rate (purple) and of the contributions of the two stars (green and blue). The observational value of the corresponding parameter and its error bars are represented by the solid red line and the dashed red horizontal lines, respectively.

In the text

	Fig. 14. Evolution as a function of stellar age of the contributions $\frac{c_1}{c_1+c_2}$ and $\frac{c_2}{c_1+c_2}$ to the weighted-average mean ${\overline{k}}_2$ of the internal structure constants k2, 1 and k2, 2 (to which the empirical correction Eq. (22) has been applied) for the models presented in Fig. 13.
In the text

	Fig. 15. Degeneracy between the various processes in the stellar interior for the best-fit min-Clés models: turbulent diffusion coefficient DT as a function of the mass-loss rate scaling parameter ξ for different values of the overshooting parameter αov. The colours stand for αov = 0.20 (green), 0.25 (blue), 0.30 (purple), 0.35 (pink), and 0.40 (orange).
In the text

Fig. 16. Hertzsprung-Russell diagram: evolutionary tracks of Clés models corresponding to the following best-fit models: Model I (coral), Series II(1) (green), Series II(5) (light blue), Series III(1) (dark blue), and Series IV(1) (plum). The dots and triangles over-plotted on the corresponding tracks correspond to the models for which k2, 1 is equal to ${\overline{k}}_2$ and 0.00139, respectively, while the open circles correspond to the models for which k2, 1 is equal to ${\overline{k}}_2\pm \mathrm{\Delta }{\overline{k}}_2$. The observational value is represented by the red point, and its error bars are represented by the dark red parallelogram.

In the text

Fig. 17. Hertzsprung-Russell diagram: evolutionary tracks of Clés models corresponding to the following best-fit models: Series II(1) (green), Model V (purple), and Model IX (pink). The dots and triangles over-plotted on the corresponding tracks correspond to the models for which k2, 1 is equal to ${\overline{k}}_2$ and 0.00139, respectively, while the open circles correspond to the models for which k2, 1 is equal to ${\overline{k}}_2\pm \mathrm{\Delta }{\overline{k}}_2$. The stars on the purple and pink tracks correspond to the best-fit models Model V and Model IX, respectively. The observational value is represented by the red point, and its error bars are represented by the dark red parallelogram.

In the text

Fig. 18. Hertzsprung-Russell diagram: evolutionary tracks of Clés models for the primary star of Minit = 17 M⊙ (green), 18 M⊙ (light blue), 19 M⊙ (dark blue), 20 M⊙ (pink), and 21 M⊙ (orange), and DT = 0 cm2 s−1 (solid line), 2 × 106 cm2 s−1 (dashed line), 6 × 106 cm2 s−1 (dotted-dashed line), and 10 × 106 cm2 s−1 (dotted line). All models have αov = 0.30 and ξ = 1. The dots over-plotted on the corresponding tracks correspond to the models for which k2, 1 is equal to ${\overline{k}}_2$, while the open circles correspond to the models for which k2, 1 is equal to ${\overline{k}}_2\pm \mathrm{\Delta }{\overline{k}}_2$. The filled and open squares over-plotted on the tracks having an initial mass of 19 M⊙ correspond to the models for which the mass equals the observational value within the error bars. The observational value is represented by the red point, and its error bars are represented by the dark red parallelogram.

In the text

Fig. 19. Hertzsprung-Russell diagram: evolutionary tracks of Clés models for the secondary star of Minit = 7.5 M⊙ (water green), 7.7 M⊙ (dark blue), and 8.0 M⊙ (orange), and DT = 0 cm2 s−1 (solid line), 2 × 105 cm2 s−1 (dashed line), and 5 × 105 cm2 s−1 (dotted-dashed line). All models have αov = 0.30 and ξ = 0.1. The filled 5-branch stars over-plotted on the tracks correspond to the models for which the age is equal to the age of the best-fit models of the primary star (i.e. about 9 Myr). The observational value is represented by the red point, and its error bars are represented by the dark red parallelogram.

In the text

	Fig. 20. Behaviour of Cmis, 1 and Cmis, 2 as a function of α1 and α2, the angle between the primary and secondary stellar rotation axis, respectively, and the normal to the orbital plane.
In the text

	Fig. 21. Behaviour of ${\overline{k}}_{2,\text{mis}}$ as a function of α1 and α2, the angle between the primary and secondary stellar rotation axis, respectively, and the normal to the orbital plane.
In the text

Fig. B.1. Hertzsprung-Russell diagram: evolutionary tracks of Clés models for the primary star having Minit = 19 M⊙, αov = 0.30, ξ = 1, X = 0.715, Z = 0.010 (green), 0.015 (dark blue), and 0.020 (orange), and DT = 0 cm2 s−1 (solid line), 2 × 106 cm2 s−1 (dashed line), 6 × 106 cm2 s−1 (dotted-dashed line), and 10 × 106 cm2 s−1 (dotted line). The dots over-plotted on the corresponding tracks correspond to the models for which k2, 1 is equal to ${\overline{k}}_2$, while the open circles correspond to the models for which k2, 1 is equal to ${\overline{k}}_2\pm \mathrm{\Delta }{\overline{k}}_2$. The filled and open squares over-plotted on the tracks correspond to the models for which the mass is equal to the observed value within the error bars. The observational value is represented by the red point, and its error bars are represented by the dark red parallelogram.

In the text

Fig. B.2. Hertzsprung-Russell diagram: evolutionary tracks of Clés models for the primary star having Minit = 19 M⊙, αov = 0.30, ξ = 1, Z = 0.015, X = 0.700 (water green), 0.715 (dark blue), and 0.730 (pink), and DT = 0 cm2 s−1 (solid line), 2 × 106 cm2 s−1 (dashed line), 6 × 106 cm2 s−1 (dotted-dashed line), and 10 × 106 cm2 s−1 (dotted line). The dots over-plotted on the corresponding tracks correspond to the models for which k2, 1 is equal to ${\overline{k}}_2$, while the open circles correspond to the models for which k2, 1 is equal to ${\overline{k}}_2\pm \mathrm{\Delta }{\overline{k}}_2$. The filled and open squares over-plotted on the tracks correspond to the models for which the mass is equal to the observed value within the error bars. The observational value is represented by the red point, and its error bars are represented by the dark red parallelogram.

In the text

Fig. B.3. Hertzsprung-Russell diagram: evolutionary tracks of Clés models for the primary star (upper panel) and secondary star (lower panel) corresponding to the best-fit models listed in Table B.1 together with Series II(3) for the primary and Model IIIS for the secondary. All models have αov = 0.30, primary models have ξ = 1 while secondary models have ξ = 0.1. The dots over-plotted on the evolutionary tracks of the primary correspond to the models for which k2, 1 is equal to ${\overline{k}}_2$, while the open circles correspond to the models for which k2, 1 is equal to ${\overline{k}}_2\pm \mathrm{\Delta }{\overline{k}}_2$. The filled 5-branch stars over-plotted on the secondary tracks correspond to the models for which the age is equal to the age of the best-fit model of the primary star of same X and Z. The observational value is represented by the red point, and its error bars are represented by the dark red parallelogram.

In the text