© ESO 2020

1. Introduction

Located at the core of the very young and rich open cluster NGC 6231 (e.g. Sung et al. 1998; Reipurth et al. 2008; Kuhn et al. 2017), the eccentric massive binary HD 152248 is an ideal target for studying tidally induced apsidal motion. Made of two O7.5 III-II(f) stars (Sana et al. 2001; Mayer et al. 2008; Rosu et al. 2020), the twin property of the system offers a unique opportunity to probe the internal structure of the stars composing this binary. The fundamental parameters of HD 152248 have recently been re-determined by Rosu et al. (2020), who established an anomalistic orbital period of $5.{816498}_{-0.000018}^{+0.000016}$ d, an eccentricity of 0.130 ± 0.002, and an inclination of $67.6_{-0.1}^{+0.2}°$. Effective temperatures of 34 000 ± 1000 K and bolometric luminosities of (2.73 ± 0.32)×10⁵ L_⊙ were inferred for both stars, while absolute masses of $29.5_{-0.4}^{+0.5}$ M_⊙, mean stellar radii of $15.{07}_{-0.12}^{+0.08}$ R_⊙, and surface gravities of log g = 3.55 ± 0.01 (in cgs units) were derived for both stars (Rosu et al. 2020). A secular variation of the argument of periastron ω was detected by several authors (Sana et al. 2001; Nesslinger et al. 2006; Mayer et al. 2008), and Rosu et al. (2020) used an extensive set of both spectroscopic and photometric observations to re-determine the rate of apsidal motion as $1.{843}_{-0.083}^{+0.064}°$ yr⁻¹.

In a close eccentric binary, the tidal deformation of the stars leads to stellar gravitational fields that are no longer spherically symmetric. This situation induces a slow precession of the line of apsides, known as apsidal motion (Schmitt et al. 2016, and references therein). The apsidal motion occurs at a rate $\dot{\omega }$, which is a function of the internal structure of the two stars (e.g. Shakura 1985; Claret & Giménez 2010). More specifically, the contribution of a binary component to the binary’s rate of apsidal motion is directly proportional to the star’s internal structure constant k₂ (Shakura 1985). This parameter is a sensitive indicator of the density stratification inside a star, and its value changes significantly as the star evolves away from the main-sequence.

Measuring $\dot{\omega }$ thus provides precious information about the internal structure of binary components. This is even more important for massive stars for which such information is rather scarce (Bulut & Demircan 2007). Therefore, the work performed here consists in constructing stellar evolution models in accordance with the fundamental parameters of the massive binary HD 152248 as re-determined by Rosu et al. (2020). With these models in hand, theoretical rates of apsidal motion are inferred and compared with the observational value. In this context, the fact that the system hosts two twin stars, hence two stars with identical k₂ values, offers a rare chance to directly infer the internal structure constant. In this way, and taking into account the small relativistic contribution to $\dot{\omega }$, Rosu et al. (2020) derived an observational value for the k₂ of both stars of 0.0010 ± 0.0001. In principle, the comparison between observational results and theoretical expectations allows us to get an estimate of the age of the stars as well as a quantitative estimate of the internal mixing processes occurring inside the stars (Mazeh 2008).

The present work builds upon the results of Rosu et al. (2020) and aims at performing a detailed study of the apsidal motion of HD 152248 from a stellar evolution point of view. In accordance with previous studies (e.g. Stickland et al. 1996; Penny et al. 1999; Sana et al. 2001; Mayer et al. 2008), we call the star that is eclipsed during primary minimum the primary star. The paper is organised as follows. In Sect. 2, we present the stellar evolution code Clés (Code Liégeois d’Évolution Stellaire) and give an overview of the models. We briefly review the basic equations that lie behind the internal structure constant concept and the apsidal motion in Sects. 3 and 4, respectively. We solve for the apsidal motion constants and the theoretical rate of apsidal motion for the set of Clés models that best reproduce the observed stellar parameters and investigate the influence of various prescriptions for the internal mixing in Sect. 5. In Sect. 6, we build single- and binary-star evolution models with GENEC in order to assess the impact of tidal interactions and stellar rotation on the internal structure of the stars. We discuss the impact of a rotation axis misalignment, of higher-order terms in the computation of the apsidal motion rate, and of a putative ternary component in Sect. 7. Finally, we present our conclusions in Sect. 8.

2. Clés stellar evolution models

In this section, we present the stellar evolution models computed with the Code Liégeois d’Évolution Stellaire¹ (Clés, Scuflaire et al. 2008). We give an overview of the main features of Clés (Sect. 2.1) and present the models built in the context of the present study (Sect. 2.2).

2.1. Characteristics of the code

The Clés code allows us to build stellar structure and evolution models from the Hayashi track to the beginning of the helium combustion phase for low mass stars, and up to the asymptotic giant branch phase for intermediate mass and massive stars. For the present study, we used version 19.1 of the Clés code.

The standard version of Clés uses OPAL opacities (Iglesias & Rogers 1996) combined with those of Ferguson et al. (2005) at low temperatures. The solar chemical composition is adopted from Asplund et al. (2009). The equation of state is implemented according to Cassisi et al. (2003), and the rates of nuclear reaction are implemented according to Adelberger et al. (2011). Mixing in convective regions is parameterised according to mixing length theory (Cox & Giuli 1968). The code computes the interior models; the atmosphere is computed separately and added to the model as a boundary condition. For massive stars, the internal mixing is restricted to the convective core unless overshooting or/and turbulent mixing is/are introduced in the models. In Clés, instantaneous overshooting displaces the boundary of a mixed region from r_c to r_ov, as given by the following step-function:

$$\begin{array}{c}\hfill r_{}\text{ov}=r_{\mathrm{c}}+{\alpha }_{}\text{ov}min(H_{\mathrm{p}}(r_{\mathrm{c}}),h),\end{array}$$(1)

where h is the thickness of the convective zone,

$$\begin{array}{c}\hfill H_{\mathrm{p}}(r_{\mathrm{c}})=\frac{P(r_{\mathrm{c}})}{\rho (r_{\mathrm{c}})g(r_{\mathrm{c}})}\end{array}$$(2)

is the pressure scale height taken at the edge of the convective core as given by the Schwarzschild or Ledoux criterion, P(r_c), ρ(r_c), and g(r_c) are the pressure, density, and local gravity at the edge of the convective core, respectively, and α_ov is the overshooting parameter chosen by the user. This classical step-function is adopted to prevent having a region of extra mixing larger than the initial size of the convective core. Adding core overshooting results in the enlargement of the centrally mixed core leading, therefore, to an extended main-sequence lifetime of the star as more fuel is available for the hydrogen fusion. Stellar wind mass-loss is included through the prescription of Vink et al. (2001), to which the user can apply a multiplicative scaling factor ξ. Rotational mixing is not included in the code, but its effect on the internal structure can be simulated by means of the turbulent diffusion (see Sect. 6). Turbulent transport can be taken into account in the models. It is modelled as a diffusion process by adding a pure diffusion term of the form

$$\begin{array}{c}\hfill D_{\mathrm{T}}\frac{\partial ln{X}_i}{\partial r},\end{array}$$(3)

with a negative turbulent diffusion coefficient D_T, to each element’s diffusion velocity, which tends to reduce its abundance gradient (Richer et al. 2000). In this expression, r is the radius and X_i is the mass fraction of element i at that point. The turbulent diffusion coefficient D_T (in cm² s⁻¹) takes the following expression:

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

where D_turb, n, and D_ct are chosen by the user. No microscopic diffusion is introduced in the models.

In the context of our present study, the main limitation of the code resides in the fact that the star is assumed to be single. Therefore, the putative influence of binarity on the internal structure of the stars is not taken into account. We will return to this point in Sect. 6.

2.2. Models

For the purpose of this work, we built only one set of stellar evolution models for the two stars given the twin characteristic of the system. As noted in Sect. 1, the absolute present-day mass of each star is $29.5_{-0.4}^{+0.5}{M}_{\odot }$. As the stars show clear evidence of mass-loss through stellar winds, the initial mass and the mass-loss prescription are clearly key parameters. While observational determinations of mass-loss rates in massive stars, accounting for the effects of clumping, often yield values that are lower than those predicted by the Vink et al. (2001) formalism, we consider the scaling factor of the mass-loss recipe ξ as a model parameter to be adjusted, with the main constraint being the upper limit on the present-day mass-loss rate Ṁ_unclumped (without clumping) determined by Rosu et al. (2020). At this stage, we note that the mass-loss rate is poorly constrained by the observations. Indeed, Rosu et al. (2020) could only determine an upper-bound limit of Ṁ_clumped = 8 × 10⁻⁷ M_⊙ yr⁻¹ (or Ṁ_unclumped = 2.5 × 10⁻⁶ M_⊙ yr⁻¹).

For our reference models, we adopted X = 0.715 and Z = 0.015, in accordance with Asplund et al. (2009), but we also investigated the influence of metallicity, testing Z values ranging from 0.013 to 0.017 (see Sect. 5.2). We adopted the Ledoux convection criterion. For the overshooting, we adopted α_ov = 0.20 as a reference. Since there are no obvious constraints on α_ov, we considered values from 0.10 to 0.40 with a step of 0.05 and identified the impact of this parameter on the evolution of the stars and its interrelation with the turbulent diffusion. We assumed the turbulent diffusion to be independent of the local density, meaning that n = 0 in Eq. (4). Several values of D_T were tested and explicated when necessary.

Best-fit models were computed through a Levenberg-Marquardt minimisation implemented as in Press et al. (1992) in a Fortran routine called min-Clés. This routine determines the combination of the Clés models’ free parameters, including, notably, the initial mass and the stellar age, which allow for the best reproduction of a set of the star’s currently observed properties. Other parameters, such as ξ or D_T, were also used as free parameters depending on the case (see Sect. 5). The constraints were fixed by the observational values of the mass, radius, and, depending on the case, the effective temperature and luminosity.

3. Internal structure constant k₂

The internal structure constant, also known as the apsidal motion constant, is obtained from the relation

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

where η₂(R_*) is the logarithmic derivative of the surface harmonic of the distorted star expressed in terms of the ellipticity ϵ₂ and evaluated at the stellar surface R_*:

$$\begin{array}{c}\hfill {\eta }_2(R_{\ast })={\frac{\mathrm{d}ln{ϵ}_2}{\mathrm{d}lnr}|}_{r=R_{\ast }},\end{array}$$(6)

which is the solution of the Clairaut-Radau differential equation

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

with the boundary condition η₂(0) = 0 (Hejlesen 1987). In this expression, r is the current radius at which the equation is evaluated, ρ(r) is the density at distance r from the centre, and $\overline{\rho }(r)$ is the mean density within the sphere of radius r. The internal structure constant represents the density contrast between the core and the external layers of the stars. It is equal to 0.75 for a homogeneous sphere of constant density and decreases towards values of the order of 10⁻⁴ for massive stars. As the star evolves, the density of its external layers decreases compared to the core density; accordingly, k₂ also decreases, making this quantity a good indicator of the evolutionary state of the stars.

We solved the differential equation Eq. (7) for the set of stellar structure models considered in this paper. For this purpose, we used a FORTRAN code that implements a fourth-order Runge-Kutta method with step doubling. The code has been validated against polytropic models and had already been used in the study of the apsidal motion of the massive binary HD 152218 (Rauw et al. 2016).

As the Clés models do not account for the rotation of the stars, we followed Claret (1999) in correcting the derived k₂ by an amount

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

with

$$\begin{array}{c}\hfill {\lambda }_{\mathrm{s}}=\frac{2{\mathrm{\Omega }}_{\mathrm{s}}^2{R}_{\ast }^3}{3Gm},\end{array}$$(9)

where Ω_s is the observational angular rotational velocity of the star, computed based on the projected rotational velocity of the star (taken from Rosu et al. 2020), the inclination given in Sect. 1, and the observational stellar radius. By default, this correction was applied to the outcome of all models presented in this paper, except for those models that explicitly account for the effect of rotation (see Sect. 6).

3.1. A look inside the star

The internal structure constant is a simple algebraic equation of the η₂ parameter evaluated at the stellar surface (see Eq. (5)). Since η₂ is computed from the integration of the Clairaut-Radau equation throughout the stellar interior (see Eq. (7)), we can look at the evolution of η₂ inside the star, from the centre to the surface, in order to identify the layers inside the star that contribute the most to η₂. In this context, it is instructive to consider an approximation suggested by Fitzpatrick (2016), which will help us highlight the parameters contributing the most to η₂. It has been shown (Fitzpatrick 2016) that the Clairaut-Radau equation, Eq. (7), can be written in the form known as Radau’s equation:

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}r}\left[\overline{\rho }(r)r^5{[1+{\eta }_2(r)]}^{1/2}\right]=5\overline{\rho }(r)r^4\psi ({\eta }_2),\end{array}$$(10)

where

$$\begin{array}{c}\hfill \psi ({\eta }_2)=\frac{1+{\eta }_2/2-{\eta }_2^2/10}{{(1+{\eta }_2)}^{1/2}}·\end{array}$$(11)

If η₂ ≲ 1.5, then ψ(η₂) ≈ 1 is roughly constant. Hence, Eq. (10) reduces to

$$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}r}\left[\overline{\rho }(r)r^5{[1+{\eta }_2(r)]}^{1/2}\right]\approx 5\overline{\rho }(r)r^4,\end{array}$$(12)

where the right-hand side is now independent of η₂. In this case,

$$\begin{array}{c}\hfill {[1+{\eta }_2(R)]}^{1/2}=\frac{5}{\overline{\rho }(R)R^5}{\displaystyle {\int }_0^R\overline{\rho }(r)r^4\mathrm{d}r}\end{array}$$(13)

(Fitzpatrick 2016). Replacing the mean density with its expression in terms of the mass and the radius, and expressing the term $\overline{\rho }(r)$ that appears in the integral through its integral expression, we get

$$\begin{array}{c}\hfill {[1+{\eta }_2(R)]}^{1/2}=\frac{5}{M{R}^2}{\displaystyle {\int }_0^R\left({\int }_0^{r}4\pi r^{\prime 2}\rho (r\prime )\mathrm{d}r\prime \right)r\mathrm{d}r.}\end{array}$$(14)

Performing the integration by parts, we get

$$\begin{array}{c}\hfill {[1+{\eta }_2(R)]}^{1/2}=\frac{10\pi }{M{R}^2}{\displaystyle {\int }_0^R\rho (r)r^2(R^2-r^2)\mathrm{d}r.}\end{array}$$(15)

In terms of the mean density, this last equation is equivalent to

$$\begin{array}{cc}\hfill {[1+{\eta }_2(R)]}^{1/2}& =\frac{15}{2\overline{\rho }(R)R^5}{\displaystyle {\int }_0^R\rho (r)r^2(R^2-r^2)\mathrm{d}r}\hfill \\ \hfill & =\frac{15}{2\overline{\rho }(R)}{\displaystyle {\int }_0^1\rho (r){\left(\frac{r}{R}\right)}^2(1-{\left(\frac{r}{R}\right)}^2)\mathrm{d}\frac{r}{R}.}\hfill \end{array}$$(16)

This equation expresses that a change in η₂ – and hence a change in k₂ – can be related to a change in the function ρ(r)r²(R² − r²). The evolution of the normalised² density inside the star as well as the normalised function $\rho (r){\left(\frac{r}{R}\right)}^2(1-{\left(\frac{r}{R}\right)}^2)$ are presented in the top panels of Fig. 1, for models with an initial mass of 31.0 M_⊙, overshooting parameters from 0.10 to 0.40, no turbulent diffusion, ξ = 1, Z = 0.015, and ages of 3, 4, and 5 Myr. As a reminder, the approximation under which Eq. (16) is valid is η₂ ≲ 1.5, which does not hold for all the radii inside the star; this is shown in the bottom left-hand panel of Fig. 1, which represents the evolution of η₂ inside the star as computed with Eq. (7). Nonetheless, this relation can be used as a first approximation, at least for the centre of the star, that is to say, up to a radius of ≈0.2 − 0.4 depending on the age of the star. In the bottom right-hand panel of Fig. 1, we show the evolution of the derivative of η₂ with respect to the radius, $\frac{\mathrm{d}{\eta }_2}{\mathrm{d}r}$. We see that the peak in $\frac{\mathrm{d}{\eta }_2}{\mathrm{d}r}$ happens at a slightly higher radius than the peak in the function ρ(r)r²(R² − r²), indicating the direct – though not perfect – correlation between the two quantities. This peak happens at a radius close to the junction between the convective core and the radiative envelope in the overshooting region.

Fig. 1. Evolution of the normalised density (top left panel), the normalised function $\rho (r){\left(\frac{r}{R}\right)}^2(1-{\left(\frac{r}{R}\right)}^2)$ (top right panel), η2 (bottom left panel), and $\frac{\mathrm{d}{\eta }_2}{\mathrm{d}r}$ (bottom right panel) inside the star for Clés models with an initial mass of 31.0 M⊙, assuming no turbulent diffusion, ξ = 1, and Z = 0.015. Each colour range corresponds to a different age: The models of 3, 4, and 5 Myr are depicted in blue, green, and red declinations, respectively. For each colour range, the darkest shade corresponds to αov = 0.10 while the lightest one corresponds to αov = 0.40; all the intermediate shades span the range from αov = 0.10 to 0.40 in steps of 0.05. Bottom panels: results obtained from the numerical integration of Eq. (7).

4. Apsidal motion

In the simple two-body case, the rate of apsidal motion is given by the sum of a classical Newtonian contribution (N) and a general relativistic correction (GR):

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

Provided that the stellar rotation axes are aligned with the normal to the orbital plane³, the Newtonian term can be expressed according to Sterne (1939), where we only consider the contributions arising from the second-order harmonic distorsions of the gravitational potential:

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

where a is the semi-major axis, q = m₁/m₂ is the mass ratio, k_2, 1 and k_2, 2 are the apsidal motion constants of the primary and secondary stars, respectively, P_rot,1 and P_rot,2 are the rotational periods of the primary star and secondary stars, respectively, and f and g are functions of the orbital eccentricity expressed by the following relations:

$$\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}$$(19)

The Newtonian term is itself the sum of the effects induced by tidal deformation on the one hand and the rotation of the stars on the other hand.

The general relativistic contribution to the rate of apsidal motion, when only the quadratic corrections are taken into account, is given by the expression

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\dot{\omega }}_{\text{GR}}& =\frac{2\pi }{P_{\text{orb}}}\frac{3G(m_1+m_2)}{c^{2}a(1-e^2)}\hfill \\ \hfill & ={\left(\frac{2\pi }{P_{\text{orb}}}\right)}^{5/3}\frac{3{(G(m_1+m_2))}^{2/3}}{c^2(1-e^2)},\hfill \end{array}\end{array}$$(20)

where G is the gravitational constant and c is the speed of light (Shakura 1985).

Using their re-derived parameters of HD 152248, Rosu et al. (2020) inferred values of $(1.{680}_{-0.083}^{+0.064})°{\text{yr}}^{-1}$ and (0.163 ± 0.001)° yr⁻¹, respectively, for the Newtonian and general relativistic contributions to the total observational rate of apsidal motion. Taking advantage of the twin nature of the system, a value for the observational internal structure constant k_2, obs of 0.0010 ± 0.0001 was then inferred for both stars (Rosu et al. 2020).

Before we turn to the computation of the theoretical rate of apsidal motion for the stellar structure models, some remarks are necessary. Our aim is to compare the theoretical rate of apsidal motion as the stars evolve with the observed value. Therefore, the stellar radii in Eq. (18) and the masses in Eq. (20) have been taken from the Clés models at a given evolutionary stage rather than from the photometric analysis. Since the same model is adopted for both stars, both the radius and the internal structure constant are identical for both stars. Adopting the following notations,

$$\begin{array}{c}\hfill \{\begin{array}{cc}& R_1=R_2\equiv R_{\ast },\hfill \\ \hfill & m_1=m_2\equiv m,\hfill \\ \hfill & k_{2,1}=k_{2,2}\equiv k_2,\hfill \end{array}\end{array}$$(21)

and replacing q with 1, the expression of the theoretical rate of apsidal motion thus simplifies to

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \dot{\omega }=& \frac{4\pi }{P_{\text{orb}}}{k}_2{\left(\frac{R_{\ast }}{a}\right)}^5[15f(e)+g(e)P_{\text{orb}}^2(\frac{1}{P_{\text{rot,1}}^2}+\frac{1}{P_{\text{rot,2}}^2})]\hfill \\ \hfill & +{\left(\frac{2\pi }{P_{\text{orb}}}\right)}^{5/3}\frac{3{(2Gm)}^{2/3}}{c^2(1-e^2)}.\hfill \end{array}\end{array}$$(22)

In this expression, a, P_orb, P_rot,1, and P_rot,2 are taken from the spectroscopic and photometric analysis.

5. Comparison with observations

Our observational study of HD 152248 led to the determination of a number of stellar parameters that can be compared to the predictions of the evolutionary models. These quantities are the present-day mass, radius, effective temperature, luminosity, rate of apsidal motion, and internal structure constant. As an illustration, we present in Fig. 2 the evolution, as a function of the stellar age, of the mass, the radius, the effective temperature, the luminosity, the internal structure constant, and the ensuing rate of apsidal motion computed with Eq. (22) for three Clés models with initial masses M_init = 30.0, 31.0, and 32.0 M_⊙, mass-loss according to the Vink et al. (2001) recipe (i.e. with ξ = 1), Z = 0.015, α_ov = 0.20, and no turbulent diffusion. The observational value of the corresponding parameter and its error bars are represented. The evolution of the parameters for a model with an initial mass of 30.0 M_⊙ according to Claret (2019) is also presented for comparison. We note that Claret (2019) applied a scaling factor ξ = 0.1 to the mass-loss rate computed with the Vink et al. (2001) prescription. Using the same ξ and adopting the same initial mass in our Clés models, we recovered the same behaviour as predicted by the Claret (2019) model.

Fig. 2. Evolution as a function of stellar age of the mass (top left panel), radius (top right panel), effective temperature (middle left panel), luminosity (middle right panel), internal structure constant of the star (bottom left panel), and apsidal motion rate of the binary (as computed with Eq. (22) assuming both stars are described with the same Clés model, bottom right panel) for Clés models with initial masses of 30.0 (green), 31.0 (light blue), and 32.0 M⊙ (pink), overshooting parameter of 0.20, Z = 0.015, and no turbulent diffusion. Stellar mass-loss was computed according to the formalism of Vink et al. (2001), with ξ = 1. The observational value of the corresponding parameter and its error bars are represented by the solid red line and the dashed dark red horizontal lines, respectively. The dark blue line represents the Claret (2019) model with an initial mass of 30.0 M⊙ and ξ = 0.1.

We observe that, for a given Clés model, the crossing between the model value and the observational value (taking into account the uncertainties) does not happen at the same age for all the parameters. We note that it is difficult to reconcile the luminosity and the k₂ on the one hand and all other parameters on the other hand. Indeed, both the luminosity and the k₂ suggest a minimum stellar age of ∼5.3 Myr, while the other parameters suggest ages younger than ∼ 4.9 Myr. This is even clearer when looking at the evolutionary tracks of the three models in the Hertzsprung-Russell diagram in Fig. 3. Indeed, none of these tracks cross the observational box defined by the observational radius and effective temperature and their respective error bars. The three models that fit the k₂, represented by the three dots over-plotted on the corresponding tracks, lie well above the two lines delimiting the range of acceptable radii. This observation motivated us to use the min-Clés routine to search for a set of model parameters that best reproduce the observed values simultaneously (i.e. for a single value of the age).

Fig. 3. Hertzsprung-Russell diagram: evolutionary tracks of Clés models with initial masses of 30.0 (green), 31.0 (light blue), and 32.0 M⊙ (pink), overshooting parameter of 0.20, Z = 0.015, and no turbulent diffusion. Stellar mass-loss was computed according to the formalism of Vink et al. (2001), with ξ = 1. The three dots over-plotted on the corresponding tracks correspond to the models that fit the observational k2. The observational value is represented by the red point and its error bars, computed through the effective temperature and the radius, are represented by the dark red parallelogram.

As a first step, we requested the Clés model to simultaneously reproduce the observed radius and mass with the mass-loss computed according to the Vink et al. (2001) recipe (i.e. with ξ = 1), and with Z = 0.015, α_ov = 0.20, and no turbulent diffusion. This model gave us M_init = 31.5 M_⊙ and an age of 4.65 Myr (Model I). However, this model clearly fails to reproduce the observed effective temperature, luminosity, and rate of apsidal motion (see Table 1). This confirms our observation regarding the difficulties involved in simultaneously fitting all the observed quantities. Adopting ξ = 0.5 leads to a greater discrepancy (Model III), whereas adopting ξ = 2.0 slightly improves the results (Model II).

As a second step, we therefore requested the Clés model to simultaneously reproduce the observed radius, mass and location in the Hertzsprung-Russell diagram for Z = 0.015, α_ov = 0.20, and no turbulent diffusion, but now allowing ξ to vary. These trials gave us M_init = 44.9 M_⊙, an age of 3.65 Myr, and a mass-loss scaling factor of ξ = 4.14 (Model IV). While this model predicts a rate of apsidal motion closer to the observational value than previous models (see Table 1), the corresponding mass-loss rate of 4.48 × 10⁻⁶ M_⊙ yr⁻¹ exceeds the observational upper limit of the unclumped mass-loss (2.5 × 10⁻⁶ M_⊙ yr⁻¹) by nearly a factor of two.

5.1. A ξ – α_ov – D_T degeneracy

To overcome this problem, we investigated the influence of the turbulent diffusion and overshooting on the best-fit models. We thus requested the Clés models to simultaneously reproduce the observed mass, radius, and location in the Hertzsprung-Russell diagram for Z = 0.015. In a first series of trials (Series V), we tested the impact of turbulent diffusion – leaving D_T as a free parameter – fixing ξ to 1, 2, and 3. Doing this, we found that models of equal quality, as far as the adjustment of the stellar parameters and of $\dot{\omega }$ is concerned, can be obtained for different pairs of D_T and ξ. To confirm this degeneracy between the two parameters, we performed six other series of trials for different values of the overshooting parameter α_ov: 0.10 (Series VI), 0.15 (Series VII), 0.25 (Series VIII), 0.30 (Series IX), 0.35 (Series X), and 0.40 (Series XI). The left-hand panel of Fig. 4 illustrates the degeneracy between these two parameters: Larger ξ values yield lower absolute values of D_T for a given α_ov. Assuming a higher mass-loss rate means that the initial mass of the star is also higher. This in turn implies that when the star reaches the observed mass value, its convective core is larger than when a smaller mass-loss rate is applied. With a larger mixing, there is no need for a large turbulent mixing. Hence, if ξ increases, D_T decreases in absolute value. The dispersions of the D_T values (for a given value of ξ) are large, but they also stem from the ranges of initial masses and ages that are found for each value of ξ.

Fig. 4. Degeneracy between the various mixing processes in the stellar interior for the best-fit min-Clés models: turbulent diffusion coefficient DT (left panel) and current age of the star (right panel) as a function of the mass-loss rate scaling parameter ξ for different values of αov.

These seven series of models reveal the well-known degeneracy between α_ov and D_T: the larger α_ov, the lower the absolute value of D_T for a given value of ξ. If the overshooting parameter is increased, the mixing is larger and there is no need for a large turbulent mixing. These two parameters affect the stellar structure in a similar way, that is to say, they increase the mixing, which lowers the central concentration of metals. Hence, an increase in one of these leads to a decrease in the other. This degeneracy is also illustrated in the left-hand panel of Fig. 4. The results are given in Table 1.

The only difference between the models lies in the initial mass and current age of the star: the higher the ξ, the higher the initial mass of the star and the lower its current age. This is expected since the mass-loss rate increases with ξ. For a given ξ, the couples of values (α_ov, D_T) giving the best-fit to M, R, T_eff, and L all yield stars of the same initial mass. But the higher the α_ov is, the lower the current age of the star is (see the right-hand panel of Fig. 4). We observe that for a given ξ, the higher the overshooting parameter α_ov, the lower the k₂-value and the $\dot{\omega }$ of the best-fit model (see Fig. 5). Altogether, these results highlight two main effects of a higher overshooting parameter: It accelerates the evolution of the star and increases the density contrast between the core and the external layers of the star, thereby decreasing the internal structure constant k₂.

Fig. 5. Internal structure constant k2 (left panel) and apsidal motion rate $\dot{\omega }$ (right panel) as a function of the mass-loss rate scaling parameter ξ for different values of αov. The colour code is the same as in Fig. 4.

This degeneracy of the best-fit models is presented in the Hertzsprung-Russell diagram in Fig. 6; here, the Series V(1), Series V(3), and Series XI(1) models are presented together with the observational box, which is defined by the observational radius and effective temperature and their respective error bars. A sequence of M_init = 33.0 M_⊙, α_ov = 0.40, ξ = 1, Z = 0.017, and D_T = −1.54 × 10⁷ cm² s⁻¹ is also presented for comparison (see discussion in Sect. 5.2). While the sequences corresponding to the Series V(1) and XI(1) models are very similar, the sequence associated with Series V(3) model follows a very different track. Despite these differences, they all cross the observational box at some stage of their evolution.

Fig. 6. Hertzsprung-Russell diagram: evolutionary tracks of Clés models of Minit = 40.0 M⊙, αov = 0.20, ξ = 3, Z = 0.015, and DT = −0.74 × 107 cm2 s−1 (green); Minit = 32.6 M⊙, αov = 0.20, ξ = 1, Z = 0.015, and DT = −1.91 × 107 cm2 s−1 (light blue); Minit = 32.6 M⊙, αov = 0.40, ξ = 1, Z = 0.015, and DT = −1.23 × 107 cm2 s−1 (pink); and Minit = 33.0 M⊙, αov = 0.40, ξ = 1, Z = 0.017, and DT = −1.54 × 107 cm2 s−1 (orange). The dots over-plotted on the corresponding tracks correspond to the models that fit the observational k2. The observational value is represented by the red point, and its error bars, computed through the effective temperature and the radius, are represented by the dark red parallelogram.

However, as can be seen in Table 1, none of these models are able to reproduce the observational k₂ and $\dot{\omega }$, since they predict values that are too high for these parameters. Since the rate of apsidal motion is the product of k₂, the fifth power of the radius, and a term which is identical for all previously computed models (since they all have the same mass and radius), the remaining discrepancy in $\dot{\omega }$ comes only from the discrepancy in k₂. This means that the models discussed so far are too homogeneous and predict too large a value of k₂. This situation is reminiscent of the conclusion reached by several authors (Claret & Giménez 1992, 1993; Claret 1995, 1999, 2004, 2019), who also found that the theoretical k₂ in their models was often too large.

Hence, as a third step, we relaxed our constraints on the mass, radius, effective temperature, and luminosity to find a model that reproduces these four constraints within the error bars, as well as k₂ and $\dot{\omega }$. So far, we have considered the influence of three parameters: ξ, α_ov, and D_T. Regarding ξ, the Vink et al. (2001) recipe has been found to often overpredict the actual mass-loss rates of O-type stars (Sundqvist et al. 2019, and references therein). Therefore, having a ξ factor larger than unity is rather unlikely, and from now on, we will fix ξ = 1. Regarding α_ov, while some studies suggest that α_ov could take a value of 0.30 for stars more massive than 4 M_⊙ (Claret 2007; Claret & Torres 2016, 2019), there is currently no consensus on the value of this parameter for massive stars, except that it is most likely larger than or equal to 0.20. Likewise, the values of D_T in massive stars are unconstrained.

All other parameters fixed, we need to increase α_ov and/or D_T to obtain a less homogeneous star and, hence, a lower k₂-value and a lower $\dot{\omega }$. We performed five series (Series XII, XIII, XIV, XV, and XVI) of tests with ξ = 1 and the value of α_ov fixed (to 0.20, 0.25, 0.30, 0.35, and 0.40, respectively). For each series, we ran six models with different values of D_T, starting from the best-fit value determined previously and adopting increasingly negative values of D_T. We then used min-Clés to find the initial mass and current age of the model that best reproduces the observed mass, radius, effective temperature, and luminosity for this pair of α_ov and D_T. The resulting values of M, R, T_eff, L, k₂, and $\dot{\omega }$ are presented in Fig. 7. For each individual model, we then computed a χ² based on the six parameters (M, R, T_eff, L, k₂, and $\dot{\omega }$). For each series of models with the same α_ov, we determined the minimum of the χ² as a function of D_T by fitting a cubic spline to the χ² values of the six models. We re-normalised this function by its minimum so that the reduced ${\chi }_{\nu }^2={\chi }^2/{\chi }_{\text{min}}^2$ is now equal to 1 at minimum. The 1σ uncertainty on D_T was then estimated from the D_T values corresponding to ${\chi }_{\nu }^2=2$ (but see the caveats in Andrae 2010 and Andrae et al. 2010). The results are presented in Fig. 8. In this way, we inferred couples of values for α_ov and D_T, which allowed us to obtain the best-fit models to the six parameters (Models XII, XIII, XIV, XV, and XVI). These final models are summarised in Table 1.

Fig. 7. Position of the best-fit min-Clés models: in the (M, R)-plane (left panel), in the (Teff, L/L⊙)-plane (middle panel), and in the (k2, $\dot{\omega }$)-plane (right panel) for the five Series (see text): XII (αov = 0.20, green), XIII (αov = 0.25, blue), XIV (αov = 0.30, violet), XV (αov = 0.35, pink), and XVI (αov = 0.40, orange). Each symbol corresponds to a different value of DT: The lowest absolute value (which corresponds to the best-fit of M, R, Teff, and L) is depicted by an open circle, while increasing absolute values of DT are depicted, in ascending order, by a star, an open triangle, a filled dot, a cross, and an open square. The red square represents the observational box, that is to say, the range of observational values taking into account the error bars.

Fig. 8. Reduced ${\chi }_{\nu }^2$ as a function of DT for the five Series (see text): XII (αov = 0.20, green), XIII (αov = 0.25, blue), XIV (αov = 0.30, violet), XV (αov = 0.35, pink), and XVI (αov = 0.40, orange).

We note that the profiles of the normalised density, of the normalised function $\rho (r){\left(\frac{r}{R}\right)}^2(1-{\left(\frac{r}{R}\right)}^2)$, of η₂, and of $\frac{\mathrm{d}{\eta }_2}{\mathrm{d}r}$ inside the star are not significantly affected by a modification of D_T with all other parameters fixed. As was already mentioned in Sect. 3.1 and outlined in Fig. 1, a change in α_ov with all other parameters fixed does not significantly affect these profiles either. These profiles are indeed mainly affected by the evolution of the star (i.e. its age).

These five final models all give an initial mass for the star of 32.8 ± 0.6 M_⊙ and a current age ranging from 5.12 to 5.18 Myr, which gives an age estimate of 5.15 ± 0.13 Myr for the stars. In order for the mass-loss rate and the turbulent diffusion to take less extreme values, the overshooting parameter has to take a high value of 0.40. If the mass-loss rates were higher than predicted by the Vink et al. (2001) recipe, the initial masses of the stars would be higher and their current age would be younger.

5.2. Influence of the metallicity

There is a priori no observational evidence that HD 152248 would have a metallicity that significantly deviates from solar although Baume et al. (1999) drew attention to the work of Kilian et al. (1994) who found that the CNO mass fraction of a sample of ten B-type stars in NGC 6231 was sub-solar (between 0.7 and 1.3% compared to the then-accepted value of 1.8% for the Sun). Hence, in order to highlight the influence of the metallicity on the binary stars’ parameters, we tested values for Z of 0.013, 0.014, 0.016, and 0.017 for models that have ξ = 1 and different values of the overshooting parameter, namely 0.20, 0.25, 0.30, 0.35, and 0.40. We left the age, initial mass, and D_T as free parameters and constrained M, R, T_eff, and L. The resulting values of these four parameters were all in agreement with the observational values.

The initial mass and the stellar age as well as the turbulent diffusion coefficient D_T are presented as a function of the overshooting parameter α_ov for different values of the metallicity Z in Fig. 9. For a given overshooting parameter α_ov, an increase in Z leads to an increase in the absolute value of D_T, a decrease in the current age, and an increase in the initial mass of the star. The increase in M_init with Z is expected since a higher Z-value implies a higher mass-loss rate, and, combined with the fact that all models have similar ages, a higher initial mass is required. The impacts on the initial mass and on the age are nonetheless quite modest.

Fig. 9. Initial mass Minit of the star (left panel), current age of the star (middle panel), and turbulent diffusion coefficient DT (right panel) as a function of αov for different values of the metallicity Z.

However, the impact of a modification of Z on the internal structure constant k₂ and on the apsidal motion rate $\dot{\omega }$ are much more significant, as Fig. 10 indicates. For a given mass, increasing the metallicity increases the opacity. Hence, the radiative transport is less effective and both the effective temperature and luminosity decrease. Therefore, to satisfy the constraints on the mass, radius, effective temperature, and luminosity, the adjustments show an increase, in absolute value, in D_T. Indeed, an increase in D_T (in absolute value) induces an increase in the luminosity, which compensates for the decrease in luminosity induced by an enhanced metallicity. The decrease in k₂ with metallicity directly follows from the increased mass of the convective core with turbulent diffusion, which in turns induces a higher density contrast inside the star (see Sect. 6, Fig. 14, where the extent of the convective core is seen as the central region of constant X). The decrease in $\dot{\omega }$ with Z directly follows the decrease in k₂. We note that the modification of the various parameters with metallicity is not a direct consequence of the modification of the metallicity, but rather an indirect consequence following from the modification of the turbulent diffusion with metallicity. The degeneracy between Z and D_T is illustrated in the Hertzsprung-Russell diagram in Fig. 6 for an evolutionary sequence of M_init = 33.0 M_⊙, α_ov = 0.40, Z = 0.017, and D_T = −1.54 × 10⁷ cm² s⁻¹. The evolutionary track of this sequence overlaps almost perfectly with the evolutionary track of a sequence of M_init = 32.6 M_⊙, α_ov = 0.40, Z = 0.015, and D_T = −1.23 × 10⁷ cm² s⁻¹.

Fig. 10. Internal structure constant k2 (left panel) and apsidal motion rate $\dot{\omega }$ (right panel) as a function of αov for different values of the metallicity Z. The colour code is the same as in Fig. 9.

6. Binary star evolution models

In this section, we go one step further and build stellar evolution models with the GENEC code (Eggenberger et al. 2008). GENEC is a one-dimensional code that accounts for the two-dimensional stellar surface deformation due to rotation. In its binary version, this code allows the inclusion of the effects of tidal interactions on the stellar structure, and thus on the k₂ value.

The GENEC code exists in two versions: the “single” version for an isolated star (Eggenberger et al. 2008) and the “binary” version for a star belonging to a binary system (Song et al. 2013). The binary version of the code currently follows the evolution of the stars until the onset of mass transfer. It accounts for the mixing induced by tidal interactions (Song et al. 2013) but does not account for the temporal modulation of these effects due to the eccentricity. The mass-loss rate is implemented through the Vink et al. (2001) prescription, as in the Clés code. In GENEC, like in Clés, we used the Ledoux criterion for convection.

We built three sequences of stellar evolution models: two sequences with the single version (GENEC single) and one sequence with the binary version (GENEC binary). All sequences assume an initial mass of 32.8 M_⊙, α_ov = 0.20, ξ = 1, Z = 0.015, and no turbulent diffusion. For the binary version and one single version (hereafter “single-v₁₀₀”), an initial equatorial rotational velocity v_eq,ini of 100 km s⁻¹ was chosen considering that the orbital period corresponds to an equatorial rotational velocity of approximately 60 km s⁻¹. For the second single version (hereafter “single-v₀”), we adopted v_eq,ini = 0 km s⁻¹. The derived k₂ value of the single-v₀ model was corrected by the amount given by the empirical correction proposed by Claret (1999, see Eq. (8), Sect. 3). The binary star sequence reaches the mass transfer state at time t = 5.34 Myr. At this stage, we emphasise that the Clés and GENEC codes adopt different definitions for the zero age of the stars. While the age in Clés is computed including the pre-main sequence phase, in GENEC the ages are taken from the ZAMS, where a model is considered to be on the ZAMS when 3‰ of the initial hydrogen content has been burned in the core.

In Fig. 11 we present the evolution as a function of the stellar age of the mass, radius, effective temperature, luminosity, and internal structure constant of the star for these three GENEC models as well as the evolution of $\dot{\omega }$ as a function of the age as computed with Eq. (22), assuming both stars are described by the same GENEC model. Two Clés models with an initial mass of 32.8 M_⊙, α_ov = 0.20, ξ = 1, and Z = 0.015 – one with no turbulent diffusion and one with a turbulent diffusion coefficient D_T = −2.41 × 10⁷ – are also presented for comparison. The observational value of the corresponding parameter and its error bars are also represented.

Fig. 11. Evolution as a function of stellar age of the mass (top left panel), radius (top right panel), effective temperature (middle left panel), luminosity (middle right panel), internal structure constant of the star (bottom left panel), and apsidal motion rate of the binary (as computed with Eq. (22) assuming both stars are described with the same model, bottom right panel) for: three GENEC models with an initial mass of 32.8 M⊙, overshooting parameter 0.20, and Z = 0.015; the single version (blue) and the binary version (green) with initial rotational velocities veq,ini = 100 km s−1; and the single version (violet) with zero initial rotational velocity (veq,ini = 0 km s−1). Two Clés models with an initial mass of 32.8 M⊙, αov = 0.20, and Z = 0.015, one with no turbulent diffusion (pink) and one with DT = −2.41 × 107 cm2 s−1 (orange), are also presented for comparison. Stellar mass-loss was computed according to the Vink et al. (2001) formalism with ξ = 1 for all models. The observational values of the corresponding parameters and their error bars are represented by the solid red line and the dashed dark red horizontal lines, respectively.

Large differences between the GENEC single-v₁₀₀ and the GENEC single-v₀ models are observed. Indeed, for a given age, the GENEC single-v₁₀₀ model is less massive, smaller, more homogeneous, and has a higher effective temperature and luminosity. The higher degree of homogeneity of the model with rotation is explained by the fact that the rotation induces additional mixing inside the star, which flattens the density and η profiles, from the core to the external layers. Indeed, for a given age, the peak of the function $\rho (r){\left(\frac{r}{R}\right)}^2(1-{\left(\frac{r}{R}\right)}^2)$ happens further away from the centre of the star and is less steep when rotation is included in the model. Consequently, the same is observed for the $\frac{\mathrm{d}{\eta }_2}{\mathrm{d}r}$ profile inside the star. These effects are even more evident as the age of the model increases. The higher k₂ value does not compensate for the smaller radius of the star, and, consequently, the apsidal motion rate is significantly reduced compared to the GENEC single-v₀ model.

One important question to address here is whether the empirical correction proposed by Claret (1999) for k₂ (see Eq. (8), Sect. 3) to account for the effect of stellar rotation is appropriate. As explained by Claret (1999), the comparison of k₂ values between models with and without rotation applies to models of the same evolutionary state as measured by the hydrogen fraction in the core X_c. For models involving different physics, this implies that we have to compare models at different ages, as the age when a given value of X_c is reached strongly depends on the internal mixing processes. Hence, Fig. 12 shows the evolution of k₂ as a function of X_c for the GENEC single-v₁₀₀ and GENEC single-v₀ models as well as the GENEC single-v₀ model to which the empirical correction of Claret (1999) has not been applied (single-v₀, uncorrected). From these evolutions, we can assert that the empirical correction proposed by Claret (1999) for k₂ to account for the effect of stellar rotation is not sufficient for the stars considered here. The empirical correction effectively leads to a small reduction of k₂ for a given X_c, but, as Fig. 12 indicates, this reduction is significantly smaller than the actual effect of rotation on k₂, especially during the early stages of the evolution.

Fig. 12. Evolution of k2 as a function of the hydrogen mass fraction Xc inside the star for three GENEC models with an initial mass of 32.8 M⊙, overshooting parameter 0.20, Z = 0.015, and ξ = 1; the single version with initial rotational velocity veq,ini = 100 km s−1 (blue); the single version with zero initial rotational velocity (veq,ini = 0 km s−1, violet); and the single version with zero initial rotational velocity to which the empirical correction for k2 of Claret (1999) has not been applied (red).

Compared to the Clés model without turbulent diffusion, the GENEC single-v₀ model gives a slightly more massive, more voluminous, cooler, more luminous, and less homogeneous star. The higher radius dominates the smaller k₂, hence the higher apsidal motion rate. The differences between the GENEC single-v₀ and Clés models are most probably due to inherent differences in the implementation of the codes. Comparing the Clés models with the GENEC single-v₁₀₀ model gives an interesting result: The GENEC single-v₁₀₀ model has an evolution in terms of absolute parameters of the stars intermediate between Clés models with and without turbulent diffusion.

This is even clearer when looking at the evolutionary tracks of the five models in the Hertzsprung-Russell diagram in Fig. 13. Indeed, all the tracks, except the GENEC single-v₀ one, cross the observational box defined by the observational radius and effective temperature and their respective error bars. The five models that fit the k₂ are represented by the five dots over-plotted on the corresponding tracks. All of these, except the one corresponding to the Clés model with enhanced turbulent diffusion, lie well above the range of acceptable radii. This highlights, again, the need for an enhanced turbulent mixing inside the stars to reproduce the small observational k₂ value. Nonetheless, the two GENEC models with rotation that fit the k₂ observational value are closer to the observational box than the GENEC model without rotation, suggesting that, at first glance, the stellar rotation acts on the stellar interior in the same way as turbulent diffusion does.

Fig. 13. Hertzsprung-Russell diagram: evolutionary tracks of three GENEC models with an initial mass of 32.8 M⊙, overshooting parameter 0.20, Z = 0.015, and ξ = 1 for the single version (blue) and binary version (green) – both with initial rotational velocity veq,ini = 100 km s−1 – and the single version with zero initial rotational velocity (veq,ini = 0 km s−1, violet). Two Clés models with an initial mass of 32.8 M⊙, αov = 0.20, Z = 0.015, and ξ = 1 – one without turbulent diffusion (pink) and one with DT = −2.41 × 107 cm2 s−1 (orange) – are also presented for comparison. The dots on the corresponding tracks correspond to the models that fit the observational k2. The red point indicates the observationally determined value, while the dark red parallelogram indicates the associated error box.

Figure 14 illustrates the hydrogen mass fraction profile inside the star for these five models at an age when the central fraction of hydrogen reaches a value of X_c = 0.2. The hydrogen profile of the Clés model with turbulent diffusion is closer to that of the GENEC models with rotation than the model without rotation. However, while the Clés model with enhanced turbulent diffusion has a hydrogen abundance at the stellar surface identical to the ones of the Clés model without turbulent diffusion and the GENEC model without rotation, the surface hydrogen abundance of the two GENEC models with rotation are lower: Stellar rotation induces a decrease in X at the stellar surface. This suggests that there is a qualitative similarity in the effects caused by turbulent diffusion inside the star on the one hand and stellar rotation on the other hand. We therefore suspect that the high turbulent diffusion coefficient found for the stars in Sect. 5.1 could be partly attributed to the effect of stellar rotation, but further investigation is required to confirm this hypothesis. Comparing the GENEC single-v₁₀₀ and GENEC binary curves, we observe that the inclusion of binarity in the code induces a reduction in the size of the convective core and an increase in the hydrogen mass fraction at the stellar surface.

Fig. 14. Evolution of the hydrogen mass fraction X as a function of the mass fraction inside the star for: three GENEC models with an initial mass of 32.8 M⊙, overshooting parameter 0.20, Z = 0.015, and ξ = 1; the single version (blue) and the binary version (green) with initial rotational velocities of veq,ini = 100 km s−1; and the single version (violet) with zero initial rotational velocity (veq,ini = 0 km s−1). Two Clés models with an initial mass of 32.8 M⊙, αov = 0.20, Z = 0.015, and ξ = 1 – one without turbulent diffusion (pink) and one with DT = −2.41 × 107 cm2 s−1 (orange) – are also presented for comparison. All models have a central fraction of hydrogen Xc of 0.2.

We now compare the single-v₁₀₀ and binary GENEC models and observe that the mass and the luminosity are not affected much: Until the onset of mass-transfer, both models have essentially identical masses and luminosities. For a given age, the binary model gives a slightly bigger and less homogeneous star with a lower effective temperature than the single model. Again, the smaller value of k₂ does not compensate for the higher radius, and, consequently, the apsidal motion rate at a given age is slightly higher for the binary version. Due to the presence of a companion and the induced tidal interactions, a binary star does not evolve on the same timescale as a single star. Therefore, it is interesting to look at the evolution of k₂ as a function of the radius of the star rather than as a function of the age. In Fig. 15, we present the evolution of k₂ as a function of the radius for the three GENEC models and the two Clés models. We observe that the single-v₁₀₀ and binary GENEC models overlap perfectly, the Clés model with turbulent diffusion is closer to the two GENEC models with rotation, and the GENEC single-v₀ and the Clés model without turbulent diffusion overlap almost perfectly. This means that whatever the ages of the single-v₁₀₀ and binary GENEC models, if they have the same radius then they have the same k₂, and, hence, the same density profile inside the star. As a consequence, the binarity property implemented in the code acts as if it would only change the age at which the star reaches a given state characterised by a given couple of values for R and k₂.

Fig. 15. Evolution of k2 as a function of stellar radius for: three GENEC models with an initial mass of 32.8 M⊙, overshooting parameter 0.20, Z = 0.015 and ξ = 1; the single version (blue) and the binary version (green) with initial rotational velocities of veq,ini = 100 km s−1; and the single version with zero initial rotational velocity (veq,ini = 0 km s−1, violet). We note that the green and blue curves overlap perfectly until the green curve stops around ∼20 R⊙. Two Clés models with an initial mass of 32.8 M⊙, αov = 0.20, Z = 0.015, and ξ = 1 – one without turbulent diffusion (pink) and one with DT = −2.41 × 107 cm2 s−1 (orange) – are also presented for comparison.

7. Additional effects

In this section, we analyse the impact of some effects that could bias our interpretation of the apsidal motion rate in terms of the internal structure constant. We consider the inclination of the stars’ rotation axes with respect to the normal to the orbital plane (Sect. 7.1), higher-order terms in the expression of the apsidal motion rate (Sect. 7.2), and the possible action of a ternary component (Sect. 7.3).

7.1. Rotation axes misalignment

The general expression of the Newtonian contribution to the total rate of apsidal motion when the rotation axes are randomly oriented is 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)\{\frac{k_{2,1}}{q}{\left(\frac{R_1}{a}\right)}^5+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}\frac{1+q}{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}(1+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}$$(23)

where α₁ (resp. α₂) is the angle between the primary (resp. secondary) star rotation axis and the normal to the orbital plane, and β₁ (resp. β₂) is the angle between the primary (resp. secondary) star rotation axis and the line joining the binary centre and the observer. Again, we adopt the notations given by Eq. (21). If we further assume the rotation axis of the two stars to be parallel, we have

$$\begin{array}{c}\hfill \{\begin{array}{cc}& {\alpha }_1={\alpha }_2\equiv \alpha ,\hfill \\ \hfill & {\beta }_1={\beta }_2\equiv \beta .\hfill \end{array}\end{array}$$(24)

Then Eq. (23) simplifies to

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\dot{\omega }}_{\text{N}}=& \frac{4\pi }{P_{\text{orb}}}{k}_2{\left(\frac{R_{\ast }}{a}\right)}^5[15f(e)-\frac{g(e)}{{sin}^{2}i}{P}_{\text{orb}}^2(\frac{1}{P_{\text{rot,1}}^2}+\frac{1}{P_{\text{rot,2}}^2})\hfill \\ \hfill & [cos\alpha (cos\alpha -cos\beta cosi)+\frac{1}{2}{sin}^{2}i(1-5{cos}^2\alpha )]].\hfill \end{array}\end{array}$$(25)

We define the function F_α as

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

We define θ as the azimutal angle of the rotation axes of the stars, such that θ = 0° and 180° correspond to the cases where the rotation axes r, the line of sight l, and the normal to the orbital plane n lie in the same plane (see Fig. 16). Using the cosine relationship in spherical trigonometry for β, namely

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

Fig. 16. Definition of the azimutal angle θ. The plane depicted is the plane of the binary orbit; the binary orbit is depicted in orange. The two stars are symbolised by the two orange dots. The r, l, and n are the rotation axis of the star, the line of sight, and the normal to the orbital plane, respectively, while the angles i and α are, respectively, the orbital inclination and the angle between the stellar rotation axis and the normal to the orbital plane.

we get

$$\begin{array}{c}\hfill \begin{array}{cc}& F_{\alpha }=cos\alpha (cos\alpha -{cos}^{2}icos\alpha -cosisinisin\alpha cos\theta )\hfill \\ \hfill & +\frac{1}{2}{sin}^{2}i(1-5{cos}^2\alpha ).\hfill \end{array}\end{array}$$(28)

F_α takes its extremum values for θ = 0° and 180°, in which cases it reduces to

$$\begin{array}{c}\hfill F_{\alpha ,\pm }={sin}^{2}i(\frac{1}{2}-\frac{3}{2}{cos}^2\alpha )\pm \frac{1}{4}sin(2i)sin(2\alpha ),\end{array}$$(29)

where the plus and minus signs correspond to θ = 180° and 0°, respectively. We represent F_α, + and F_α, − as a function of the angle α in Fig. 17 for the specific case of i = 67.6°. The minimum and maximum values of F_α, − are reached for α = 7.7° and α = 97.7°, respectively, while the minimum and maximum values of F_α, + are reached for α = 172.3° and α = 82.3°, respectively.

Fig. 17. Behaviour of Fα, + (in red) and Fα, − (in blue), given by Eq. (29), as a function of α, the angle between the rotation axes of the stars and the normal to the orbital plane.

We compute the total (Newtonian plus relativistic) apsidal motion rate as a function of the age for Clés models with an initial mass of 31.0 M_⊙ and α_ov = 0.20 for different angles α of the stars’ rotation axes and assuming θ = 180°. The results are presented in Fig. 18. A small misalignment angle (up to approximately 30°) has nearly no impact on the rate of apsidal motion value. For a higher value of the misalignment angle, an effect on the rate of apsidal motion is clearly seen. In the extreme case (α = 82°), the inferred age is increased by 0.5 Myr compared to the aligned case (α = 0°). The small effect on the rate of apsidal motion induced by the misalignment of the stars’ rotation axes is related to the fact that it only contributes to the term generated by the stellar rotation, a term that itself contributes, in the case of HD 152248, approximately 15% to the total Newtonian term.

Fig. 18. Evolution of the apsidal motion rate as a function of stellar age for Clés models with an initial mass of 31.0 M⊙ and αov = 0.20 for different misalignment angles α of the stellar rotation axes with respect to the normal to the orbital plane in the case where θ = 180°. The cases of α = 0° ,15° ,30° ,45° ,60° and 82° are depicted in green, light green, cyan, blue, purple, and pink, respectively. The observational value of the apsidal motion rate and its error bars are represented by the solid red line and the dashed dark red horizontal lines, respectively.

One may wonder whether the condition of sub-critical rotation velocity leads to a restriction on the values of the putative misalignment angle α between the rotational and orbital axes. In fact, the stellar rotation rate Ω_rot cannot exceed the critical value Ω_crit, which corresponds to a situation where the centrifugal force at the equator equals the effective gravity (i.e. the gravity corrected for the effect of radiation pressure):

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

Here, the Eddington factor is given by

$$\begin{array}{c}\hfill {\mathrm{\Gamma }}_{\mathrm{e}}=\frac{{\kappa }_{\mathrm{e}}{L}_{\ast }}{4\pi cG{M}_{\ast }}.\end{array}$$(31)

We take the electron scattering opacity to be κ_e ≃ 0.34 cm² g⁻¹ (e.g. Sanyal et al. 2015), implying ${\mathrm{\Gamma }}_{\mathrm{e}}=2.6{10}^{-5}\frac{L_{\ast }{M}_{\odot }}{L_{\odot }{M}_{\ast }}$.

Condition (30) translates into an upper limit on the equatorial rotational velocity

$$\begin{array}{c}\hfill v_{\mathrm{eq}}\le \sqrt{\frac{G{M}_{\ast }}{R_{\ast }}(1-\frac{{\kappa }_{\mathrm{e}}{L}_{\ast }}{4\pi cG{M}_{\ast }})}\end{array}$$(32)

or

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

The stellar parameters found by Rosu et al. (2020) lead to Γ_e = 0.2 and v_eq ≤ 546 km s⁻¹. Combining this result with the value of v_eq sin β = 138 km s⁻¹ found by Rosu et al. (2020) for the primary, we obtain β ≥ 14.6°.

The orbital plane is seen under an inclination of i = 67.6° (Rosu et al. 2020). For an external observer, the impact of a putative misalignment angle α depends on the azimuthal angle θ between the plane defined by the orbital axes and our line of sight as well as the plane defined by the rotational and orbital axes through Eq. (27).

As shown above, the largest impact of α on the rate of apsidal motion corresponds to situations where θ is either 0° or 180°. For an azimuthal angle of θ = 0°, the condition on β translates into α <  53.0° or α >  82.2°. Likewise, for θ = 180°, we obtain α <  97.8° or α >  127.0°. These conditions do not rule out the values of α that produce the largest impact on $\dot{\omega }$. Hence, the condition of sub-critical rotation cannot be used to infer a general constraint on α.

7.2. Higher-order terms in $\dot{\omega }$

In this section, we look at the higher-order contributions to the rate of apsidal motion due to tidal distortion. The first two contributions are ${\dot{\omega }}_3$ and ${\dot{\omega }}_4$, which can be expressed as

$$\begin{array}{c}\hfill {\dot{\omega }}_3=\frac{56\pi }{P_{\text{orb}}}{f}_3(e)(k_{3,1}\frac{1}{q}{\left(\frac{R_1}{a}\right)}^7+k_{3,2}q{\left(\frac{R_2}{a}\right)}^7)\end{array}$$(34)

and

$$\begin{array}{c}\hfill {\dot{\omega }}_4=\frac{90\pi }{P_{\text{orb}}}{f}_4(e)(k_{4,1}\frac{1}{q}{\left(\frac{R_1}{a}\right)}^9+k_{4,2}q{\left(\frac{R_2}{a}\right)}^9),\end{array}$$(35)

where

$$\begin{array}{c}\hfill f_3(e)=\frac{1+\frac{15}{4}{e}^2+\frac{15}{8}{e}^4+\frac{5}{64}{e}^6}{{(1-e^2)}^7}\end{array}$$(36)

and

$$\begin{array}{c}\hfill f_4(e)=\frac{1+7e^2+\frac{35}{4}{e}^4+\frac{35}{16}{e}^6+\frac{7}{128}{e}^8}{{(1-e^2)}^9},\end{array}$$(37)

as shown by Sterne (1939). In these two expressions, the internal structure constants of the stars are defined by

$$\begin{array}{c}\hfill k_j=\frac{j+1-{\eta }_j(R_{\ast })}{2(j+{\eta }_j(R_{\ast }))},\end{array}$$(38)

where the η_j are solutions of the second-order differential equation

$$\begin{array}{c}\hfill r\frac{\mathrm{d}{\eta }_j(r)}{\mathrm{d}r}+{\eta }_j^2(r)-{\eta }_j(r)+6\frac{\rho (r)}{\langle \rho \rangle (r)}({\eta }_j(r)+1)-j(j+1)=0\end{array}$$(39)

with the boundary condition η_j(0) = j − 2 (Hejlesen 1987). We have computed ${\dot{\omega }}_3$ and ${\dot{\omega }}_4$ for all the models presented previously. For each of them, ${\dot{\omega }}_3$ and ${\dot{\omega }}_4$ contribute to a maximum amount of 6 × 10⁻⁵° yr⁻¹ and 6 × 10⁻⁶° yr⁻¹, respectively, to the rate of apsidal motion. These values are, respectively, three and four orders of magnitude below the actual precision on the value of the observed rate of apsidal motion. Hence, we can ignore these two contributions, ${\dot{\omega }}_3$ and ${\dot{\omega }}_4$, when computing the theoretical rate of apsidal motion.

7.3. Action of a putative ternary star?

As pointed out by Borkovits et al. (2019), the picture of apsidal motion can become significantly more complicated in the case of a hierarchical triple system. What happens in these systems strongly depends on the relative orientation of the inner and outer orbital planes. If the two planes are essentially co-planar, the action of the ternary star on the inner binary can be summarised as an increase in the rate of apsidal motion compared to the value expected from the tidal interactions and relativistic effects (Borkovits et al. 2019; Bozkurt & Deǧirmenci 2007). In this case, the contribution of the ternary star to $\dot{\omega }$ would be positive (i.e. in the same direction as the Newtonian and general relativistic effects).

However, things can change significantly if the outer orbit is not co-planar with the inner one. Such systems are subject to the Kozai-Lidov mechanism. In these cases, the eccentricity of the inner binary undergoes cyclic variations (Söderhjelm 1982; Naoz et al. 2013; Borkovits et al. 2019). In addition, the line of apsides of the inner binary may either rotate or librate (Borkovits et al. 2019). In the former case, the contribution of the third star to the $\dot{\omega }$ of the inner binary can be positive or negative, giving either too fast or too slow an apsidal motion compared to the tidal and relativistic effects. Therefore, the action of a putative ternary star could actually bias the observational determination of the rate of apsidal motion (Borkovits et al. 2005, 2019).

Naoz et al. (2013) investigated the secular perturbations in hierarchical triple systems up to the octupole order. These authors showed that the relative importance of the octupole term of the Hamiltonian compared to the quadrupole term is proportional to

$$\begin{array}{c}\hfill {ϵ}_M=\frac{m_1-m_2}{m_1+m_2}\frac{a_{\mathrm{in}}}{a_{\mathrm{out}}}\frac{e_{\mathrm{out}}}{1-e_{\mathrm{out}}^2},\end{array}$$(40)

where the quantities a_in, a_out, and e_out stand for the semi-major axis of the inner binary, the semi-major axis of the ternary star motion around the inner binary’s centre of mass, and the eccentricity of the ternary star’s orbit. Since the eclipsing binary of HD 152248 consists of two nearly identical stars (i.e. m₁ = m₂), we conclude that, in the present case, only the quadrupole term should be relevant. As a result, the rate of apsidal motion due to a putative ternary component would be given by Eq. (A.29) of Naoz et al. (2013):

$$\begin{array}{cc}\hfill \dot{{\omega }_1}=6C_2& \{\frac{1}{G_1}[4{cos}^2{i}_{\mathrm{tot}}+(5cos2{\omega }_1-1)(1-e_1^2-{cos}^2{i}_{\mathrm{tot}})]\hfill \\ \hfill & +\frac{cos{i}_{\mathrm{tot}}}{G_2}[2+e_1^2(3-5cos2{\omega }_1)]\}.\hfill \end{array}$$(41)

In this equation, the argument of periastron of the inner (short-period) binary is denoted with ω₁, while e₁ is the eccentricity of this orbit. The G₁ and G₂ stand for the angular momentum of the inner and outer orbits. The angle i_tot is the mutual inclination between the angular momentum vectors of the inner and outer orbits. The constant C₂ gathers all the dependencies on the masses of the components and semi-major axes of the inner and outer orbits. This equation indicates that the rate of apsidal motion induced by the third star depends on the masses of all three components (via C₂) and the relative inclination of the inner and outer orbits, as well as on the value of ω itself and on the amplitudes of the angular momentum vectors of the two orbits. Solving this equation therefore requires information on the third star and the outer orbit.

In Rosu et al. (2020), we searched for spectroscopic evidence of a ternary star in HD 152248. It was shown that (1) the disentangled spectra and the photometric lightcurves could be interpreted without the need for a third light contribution, and that (2) the existing radial velocity data fail to reveal a reflex motion around the centre of mass of a triple system. Yet, both the residuals of the radial velocities about the orbital solution of the system accounting for the presence of apsidal motion and the dispersion of ω values inferred from the photometric data were quite large (Rosu et al. 2020). While the former effect was attributed to the dynamical tidal deformations and the impact of wind interactions on the stellar atmospheres, we cannot completely rule out the possibility of the presence of a third star that would be less massive and less bright than the components of the eclipsing binary. However, we currently lack the information regarding the orbit of such a putative star that would be needed to compute its impact on the apsidal motion of the inner binary.

8. Conclusion

We have analysed the eccentric massive binary HD 152248 and its apsidal motion from a theoretical point of view. We first built Clés stellar evolution models by constraining the mass, radius, effective temperature, and luminosity to the observational values. Assuming no turbulent diffusion occurs inside the star and fixing the overshooting parameter α_ov to a value of 0.20 gave rise to a model requiring an unrealistically high mass-loss rate scaling factor ξ. Considering different prescriptions for the internal mixing occurring inside the stars, we highlighted the well-known degeneracy between the effects of overshooting and turbulent diffusion. To reproduce k₂ and $\dot{\omega }$ along with the mass, radius, effective temperature, and luminosity required a significant internal mixing enhancement, either through overshooting or turbulent diffusion as both mechanisms lead to an enlargement of the convective core. However, further investigation is required to determine whether this enhanced turbulent diffusion, necessary to reproduce the low k₂ observational value, is physically representative of the stars. Indeed, the standard models, that is to say, models without enhanced turbulent mixing, that simultaneously fit the observational mass, radius, effective temperature, and luminosity all predict too large a k₂ value compared to the observations. The fact that standard models seem too homogeneous has already been suggested by previous studies (Claret 2019, and reference therein), although mostly for lower-mass systems.

Assuming the mass-loss rate is described by the Vink et al. (2001) recipe, these results allowed us to determine an initial mass of both stars of 32.8 ± 0.6 M_⊙ and a current age of 5.15 ± 0.13 Myr. A higher mass-loss rate would imply higher initial masses and younger current age of the stars. Rauw et al. (2016) performed a similar, but more limited, study for the massive binary system HD 152218, also located in the NGC 6231 cluster, and inferred an age estimate of 5.8 ± 0.6 Myr, which is compatible with our age estimate of HD 152248 within the error bars. If the stars appeared to have a sub-solar metallicity, as advocated by some authors (Kilian et al. 1994; Baume et al. 1999), then their initial masses would be slightly smaller, their age slightly older, and the required turbulent diffusion slightly lower.

We also investigated the impact of binarity and stellar rotation by means of single and binary GENEC models that account for stellar rotation. The binary version of this code takes into account the mixing induced by tidal interactions occurring inside a binary system. These tests revealed that the current evolutionary stage of HD 152248 can be described by single-star models and that the system has not yet gone through a mass-exchange episode; this is in line with the conclusions of Rosu et al. (2020). Furthermore, we suggest that the high turbulent diffusion required to reproduce the observations can be partly attributed to stellar rotation, but further investigation is required to confirm or infirm this hypothesis.

Finally, we analysed the impact of some effects that could bias our interpretation of the apsidal motion in terms of the internal structure constant. We considered the impact of a misalignment of the stellar rotation axes and deduced that unless the axes are inclined by an angle of at least ∼60° with respect to the normal to the orbital plane, the apsidal motion rate is not affected much. This comes from the fact that the rotation axis orientation with respect to the orbital plane only appears in the rotational term of the Newtonian contribution to the total apsidal motion rate, a term which contributes an amount of only ∼15% in the case of HD 152248. Such a high inclination angle seems highly unlikely in a close eccentric binary where the two stars are on the point of reaching pseudo-synchronisation. Indeed, Rosu et al. (2020) infers values of 0.87 ± 0.03 and 0.86 ± 0.04 for the ratios of rotational angular velocity to instantaneous orbital angular velocity at periastron for the primary star and secondary star, respectively. We further show that higher-order terms in the expression of the apsidal motion rate are completely negligible. A far more severe impact could arise from the action of a third star if HD 152248 turned out to be a triple system, although there is currently no observational evidence for the existence of a third component in this system.

Acknowledgments

The Liège team acknowledges support from the Fonds de la Recherche Scientifique (F.R.S. – FNRS, Belgium). The authors thank the referee for his/her suggestions and comments towards the improvement of the manuscript.

References

All Tables

All Figures

Fig. 1. Evolution of the normalised density (top left panel), the normalised function $\rho (r){\left(\frac{r}{R}\right)}^2(1-{\left(\frac{r}{R}\right)}^2)$ (top right panel), η2 (bottom left panel), and $\frac{\mathrm{d}{\eta }_2}{\mathrm{d}r}$ (bottom right panel) inside the star for Clés models with an initial mass of 31.0 M⊙, assuming no turbulent diffusion, ξ = 1, and Z = 0.015. Each colour range corresponds to a different age: The models of 3, 4, and 5 Myr are depicted in blue, green, and red declinations, respectively. For each colour range, the darkest shade corresponds to αov = 0.10 while the lightest one corresponds to αov = 0.40; all the intermediate shades span the range from αov = 0.10 to 0.40 in steps of 0.05. Bottom panels: results obtained from the numerical integration of Eq. (7).

In the text

Fig. 2. Evolution as a function of stellar age of the mass (top left panel), radius (top right panel), effective temperature (middle left panel), luminosity (middle right panel), internal structure constant of the star (bottom left panel), and apsidal motion rate of the binary (as computed with Eq. (22) assuming both stars are described with the same Clés model, bottom right panel) for Clés models with initial masses of 30.0 (green), 31.0 (light blue), and 32.0 M⊙ (pink), overshooting parameter of 0.20, Z = 0.015, and no turbulent diffusion. Stellar mass-loss was computed according to the formalism of Vink et al. (2001), with ξ = 1. The observational value of the corresponding parameter and its error bars are represented by the solid red line and the dashed dark red horizontal lines, respectively. The dark blue line represents the Claret (2019) model with an initial mass of 30.0 M⊙ and ξ = 0.1.

In the text

Fig. 3. Hertzsprung-Russell diagram: evolutionary tracks of Clés models with initial masses of 30.0 (green), 31.0 (light blue), and 32.0 M⊙ (pink), overshooting parameter of 0.20, Z = 0.015, and no turbulent diffusion. Stellar mass-loss was computed according to the formalism of Vink et al. (2001), with ξ = 1. The three dots over-plotted on the corresponding tracks correspond to the models that fit the observational k2. The observational value is represented by the red point and its error bars, computed through the effective temperature and the radius, are represented by the dark red parallelogram.

In the text

	Fig. 4. Degeneracy between the various mixing processes in the stellar interior for the best-fit min-Clés models: turbulent diffusion coefficient DT (left panel) and current age of the star (right panel) as a function of the mass-loss rate scaling parameter ξ for different values of αov.
In the text

	Fig. 5. Internal structure constant k2 (left panel) and apsidal motion rate $\dot{\omega }$ (right panel) as a function of the mass-loss rate scaling parameter ξ for different values of αov. The colour code is the same as in Fig. 4.
In the text

Fig. 6. Hertzsprung-Russell diagram: evolutionary tracks of Clés models of Minit = 40.0 M⊙, αov = 0.20, ξ = 3, Z = 0.015, and DT = −0.74 × 107 cm2 s−1 (green); Minit = 32.6 M⊙, αov = 0.20, ξ = 1, Z = 0.015, and DT = −1.91 × 107 cm2 s−1 (light blue); Minit = 32.6 M⊙, αov = 0.40, ξ = 1, Z = 0.015, and DT = −1.23 × 107 cm2 s−1 (pink); and Minit = 33.0 M⊙, αov = 0.40, ξ = 1, Z = 0.017, and DT = −1.54 × 107 cm2 s−1 (orange). The dots over-plotted on the corresponding tracks correspond to the models that fit the observational k2. The observational value is represented by the red point, and its error bars, computed through the effective temperature and the radius, are represented by the dark red parallelogram.

In the text

Fig. 7. Position of the best-fit min-Clés models: in the (M, R)-plane (left panel), in the (Teff, L/L⊙)-plane (middle panel), and in the (k2, $\dot{\omega }$)-plane (right panel) for the five Series (see text): XII (αov = 0.20, green), XIII (αov = 0.25, blue), XIV (αov = 0.30, violet), XV (αov = 0.35, pink), and XVI (αov = 0.40, orange). Each symbol corresponds to a different value of DT: The lowest absolute value (which corresponds to the best-fit of M, R, Teff, and L) is depicted by an open circle, while increasing absolute values of DT are depicted, in ascending order, by a star, an open triangle, a filled dot, a cross, and an open square. The red square represents the observational box, that is to say, the range of observational values taking into account the error bars.

In the text

	Fig. 8. Reduced ${\chi }_{\nu }^2$ as a function of DT for the five Series (see text): XII (αov = 0.20, green), XIII (αov = 0.25, blue), XIV (αov = 0.30, violet), XV (αov = 0.35, pink), and XVI (αov = 0.40, orange).
In the text

	Fig. 9. Initial mass Minit of the star (left panel), current age of the star (middle panel), and turbulent diffusion coefficient DT (right panel) as a function of αov for different values of the metallicity Z.
In the text

	Fig. 10. Internal structure constant k2 (left panel) and apsidal motion rate $\dot{\omega }$ (right panel) as a function of αov for different values of the metallicity Z. The colour code is the same as in Fig. 9.
In the text

Fig. 11. Evolution as a function of stellar age of the mass (top left panel), radius (top right panel), effective temperature (middle left panel), luminosity (middle right panel), internal structure constant of the star (bottom left panel), and apsidal motion rate of the binary (as computed with Eq. (22) assuming both stars are described with the same model, bottom right panel) for: three GENEC models with an initial mass of 32.8 M⊙, overshooting parameter 0.20, and Z = 0.015; the single version (blue) and the binary version (green) with initial rotational velocities veq,ini = 100 km s−1; and the single version (violet) with zero initial rotational velocity (veq,ini = 0 km s−1). Two Clés models with an initial mass of 32.8 M⊙, αov = 0.20, and Z = 0.015, one with no turbulent diffusion (pink) and one with DT = −2.41 × 107 cm2 s−1 (orange), are also presented for comparison. Stellar mass-loss was computed according to the Vink et al. (2001) formalism with ξ = 1 for all models. The observational values of the corresponding parameters and their error bars are represented by the solid red line and the dashed dark red horizontal lines, respectively.

In the text

Fig. 12. Evolution of k2 as a function of the hydrogen mass fraction Xc inside the star for three GENEC models with an initial mass of 32.8 M⊙, overshooting parameter 0.20, Z = 0.015, and ξ = 1; the single version with initial rotational velocity veq,ini = 100 km s−1 (blue); the single version with zero initial rotational velocity (veq,ini = 0 km s−1, violet); and the single version with zero initial rotational velocity to which the empirical correction for k2 of Claret (1999) has not been applied (red).

In the text

Fig. 13. Hertzsprung-Russell diagram: evolutionary tracks of three GENEC models with an initial mass of 32.8 M⊙, overshooting parameter 0.20, Z = 0.015, and ξ = 1 for the single version (blue) and binary version (green) – both with initial rotational velocity veq,ini = 100 km s−1 – and the single version with zero initial rotational velocity (veq,ini = 0 km s−1, violet). Two Clés models with an initial mass of 32.8 M⊙, αov = 0.20, Z = 0.015, and ξ = 1 – one without turbulent diffusion (pink) and one with DT = −2.41 × 107 cm2 s−1 (orange) – are also presented for comparison. The dots on the corresponding tracks correspond to the models that fit the observational k2. The red point indicates the observationally determined value, while the dark red parallelogram indicates the associated error box.

In the text

Fig. 14. Evolution of the hydrogen mass fraction X as a function of the mass fraction inside the star for: three GENEC models with an initial mass of 32.8 M⊙, overshooting parameter 0.20, Z = 0.015, and ξ = 1; the single version (blue) and the binary version (green) with initial rotational velocities of veq,ini = 100 km s−1; and the single version (violet) with zero initial rotational velocity (veq,ini = 0 km s−1). Two Clés models with an initial mass of 32.8 M⊙, αov = 0.20, Z = 0.015, and ξ = 1 – one without turbulent diffusion (pink) and one with DT = −2.41 × 107 cm2 s−1 (orange) – are also presented for comparison. All models have a central fraction of hydrogen Xc of 0.2.

In the text

Fig. 15. Evolution of k2 as a function of stellar radius for: three GENEC models with an initial mass of 32.8 M⊙, overshooting parameter 0.20, Z = 0.015 and ξ = 1; the single version (blue) and the binary version (green) with initial rotational velocities of veq,ini = 100 km s−1; and the single version with zero initial rotational velocity (veq,ini = 0 km s−1, violet). We note that the green and blue curves overlap perfectly until the green curve stops around ∼20 R⊙. Two Clés models with an initial mass of 32.8 M⊙, αov = 0.20, Z = 0.015, and ξ = 1 – one without turbulent diffusion (pink) and one with DT = −2.41 × 107 cm2 s−1 (orange) – are also presented for comparison.

In the text

Fig. 16. Definition of the azimutal angle θ. The plane depicted is the plane of the binary orbit; the binary orbit is depicted in orange. The two stars are symbolised by the two orange dots. The r, l, and n are the rotation axis of the star, the line of sight, and the normal to the orbital plane, respectively, while the angles i and α are, respectively, the orbital inclination and the angle between the stellar rotation axis and the normal to the orbital plane.

In the text

	Fig. 17. Behaviour of Fα, + (in red) and Fα, − (in blue), given by Eq. (29), as a function of α, the angle between the rotation axes of the stars and the normal to the orbital plane.
In the text

Fig. 18. Evolution of the apsidal motion rate as a function of stellar age for Clés models with an initial mass of 31.0 M⊙ and αov = 0.20 for different misalignment angles α of the stellar rotation axes with respect to the normal to the orbital plane in the case where θ = 180°. The cases of α = 0° ,15° ,30° ,45° ,60° and 82° are depicted in green, light green, cyan, blue, purple, and pink, respectively. The observational value of the apsidal motion rate and its error bars are represented by the solid red line and the dashed dark red horizontal lines, respectively.

In the text