Apsidal motion in (massive) binaries

• The majority of massive stars belong to binary systems:
  • Considerably affects the evolution of the stars;
  • Offers possibilities to constrain the properties of the stars.

• Interesting systems: double-line spectroscopic eclipsing binaries
  • Combine the photometric eclipses and the radial velocities obtained with spectroscopy;
  • Determine the masses and radii of the stars in a model independent way.

• Most interesting systems: binaries showing a significant apsidal motion
  • Slow precession of the line of apsides in an eccentric binary;
  • Arises from tidal interactions occurring between the stars of a close binary, interactions which are responsible for the non-spherical gravitational field of the stars;
  • Non-Keplerian movement (with a general relativistic contribution but small in the cases we study here).

• The rate of apsidal motion is directly related to the internal structure of the stars. Measuring the rate of apsidal motion hence
  
  • Provides a diagnostic of the internal mass-distribution of the stars, which is otherwise difficult to constrain;
  • Offers a test of our understanding of stellar structure and evolution.

The apsidal motion has been derived in hundreds of systems. Here, we report two large-scale studies (left and middle panels) as well as individual studies each focused on one single system (right panel).

Right panel: data from Baroch et al. 2021, A&A, 649, A64; Baroch et al. 2022, A&A, 665, A13; Claret et al. 2021, A&A, 654, A17; Marcussen & Albrecht 2022, ApJ, 933, 227; Rauw et al. 2016, A&A, 594, A33; Rosu et al. 2020, A&A, 635, A145; Rosu et al. 2022, A&A, 660, A120; Rosu et al. 2022, A&A, 664, A98; Rosu et al. 2023, MNRAS, 521, 2988; Torres et al. 2010, A&A Rev., 18, 67; Wolf et al. 2006, A&A, 456, 1077; Wolf et al. 2008, MNRAS, 388, 1836; Wolf et al. 2010, A&A, 509, A18). 

Check out our most recent works!

• If you want to read a review on the topic, this link is for you: 

• If you are interested in specific results for specific massive binaries:

• If you want a thorough analysis, check out my PhD thesis here: 

There are two different scenarii, depending if the binary is a double-line spectroscopic and/or  an eclipsing binary. We here review both cases. 

Radial velocities (RVs) of the primary (\(P\)) and secondary (\(S\)) stars: 
$$RV_P(t) = K_P [\cos(\phi(t) + \omega(t)) + e\,\cos(\omega(t))] + \gamma_P,$$

$$RV_S(t) = -K_S [\cos(\phi(t) + \omega(t)) + e\,\cos(\omega(t))] + \gamma_S,$$

where 

\(K_P\), \(K_S\): Amplitudes of the RV curves; 

\(e\): Eccentricity of the orbit;

\(\gamma_P\), \(\gamma_S\): Apparent systemic velocities;

\(\phi(t)\): True anomaly computed from Kepler's equation using the orbital period (\(P_{orb}\)), eccentriciy, and time of periastron passage (\(T_0\)). 

With apsidal motion, the longitude of periastron is given by:

$$\omega(t) = \omega_0 + \frac{d\omega}{dt} (t-T_0),$$

where 

\(\omega_0\): Value of \(\omega\) at \(t=T_0\).

We computed a set of RV curves for different values of \(\omega\) (left panel). We clearly see how the morphology of the RV curve changes with \(\omega\). On the right panel, we show an example of RV curve adjustment to data for which the change of morphology of the RV curve with time is evident. 

We computed a set of 9 theoretical lightcurves for the twin binary system HD 152248. 

Twin means that the two stars share the same fundamental parameters, that is to say, the same mass, radius, and effective temperature.

Each lightcurve is computed assuming a different value of the longitude of periastron \(\omega\), all other parameters identical. In particular, the eccentricity is equal to 0.134 while the inclination is equal to 67.6°.

These lightcurves are defined in such a way that the primary eclipse happens at phase 0. In addition, for these values of \(\omega\), the primary eclipse happens close to periastron passage. As a natural consequence, the secondary eclipse, around phase 0.5, happens close to apastron passage for this range of \(\omega\) values.  

The difference in the primary and secondary eclipses depths as well as the differences between the primary and secondary eclipses depths cannot be attributed to a difference in effective temperature of the stars. 

It is entirely attributed to the change in \(\omega\)! 

But why is that so?

If the inclination is different from 90°, the depth of the eclipses are determined by the orbital separation at conjuction phase.
Because the primary eclipse happens close to periastron passage for these values of \(\omega\) (between 0° and 180°), when the orbital separation is smaller, the primary eclipse is deeper than the secondary eclipse.
The orbital separation at conjunction phase itself depends on the orientation of the ellipse with respect to our line of sight - that is to say, \(\omega\). This is why the depths of both eclipses change with the value of \(\omega\).

Finally, the phase at which the eclipse occurs also depends upon the orientation of the ellipse, that is to say, \(\omega\).  

The phase difference between the primary and secondary eclipses can be computed through Equation (3) of Giménez & Bastero 1995, Ap&SS, 226, 99.

In the special case where the inclination \(i=90^\circ\), the equations take a very simple form: 
$$\Delta\phi = \frac{t_2-t_1}{P_{orb}} = \frac{\Psi-\sin\Psi}{2\pi},$$

with

$$ \Psi = \pi+2 \arctan \left(\frac{e}{\sqrt{1-e^2}}\cos\omega\right).$$

In the two-body problem, and neglecting any contribution from a third body, the total rate of apsidal motion is the sum of the Newtonian (N) contribution and the General Relativistic (GR) correction:

$$\dot\omega = \dot\omega_N + \dot\omega_{GR}.$$ 

The General Relativistic contribution is given by: 

$$ \dot\omega_{GR} = \left(\frac{2\pi}{P_{orb}}\right)^{5/3} \frac{3(G(m_1+m_2))^{2/3}}{c^2(1-e^2)},$$

where \(G\) is the gravitational constant, \(c\) is the speed of light, \(m_1\) and \(m_2\) are the masses of the primary and secondary star, respectively, \(e\) is the eccentricity of the orbit, and \(P_{orb}\) is the orbital period of the system (Shakura 1985, Sov. Astron. Lett., 11, 224). 

Assuming that the rotation axes of the stars are aligned with the normal to the orbital plane, the Newtonian term takes the form adopted by Sterne (1939, MNRAS, 99, 451), where only the contribution arising from the second-order harmonic distorsions of the potential is considered: 

$$ \dot\omega_N = \frac{2\pi}{P_{orb}} 15f(e) \left(k_{2,1} q \left(\frac{R_1}{a}\right)^5 + \frac{k_{2,2}}{q} \left(\frac{R_2}{a}\right)^5 \right) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~$$

$$+ \frac{2\pi}{P_{orb}} g(e)\left(k_{2,1}(1+q)\left(\frac{R_1}{a}\right)^5 \left(\frac{P_{orb}}{P_{rot,1}}\right)^2+k_{2,2}\frac{1+q}{q}\left(\frac{R_2}{a}\right)^5 \left(\frac{P_{orb}}{P_{rot,2}}\right)^2\right), $$

where \(R_1\) and \(R_2\) are the radii of the primary and secondary star, respectively, \(P_{rot,1}\) and \(P_{rot,2}\) are the rotational periods of the primary and secondary star, respectively, \(q=m_2/m_1\) is the mass ratio of the system, \(k_{2,1}\) and \(k_{2,2}\) are the apsidal motion constants, also called internal stellar structure constants, of the primary and secondary star, respecively, and \(f\) and \(g\) depend on the eccentricity through the following relations: 

$$ f(e) = \frac{1+\frac{3e^2}{2}+\frac{e^4}{8}}{(1-e^2)^5},$$

$$ g(e) = \frac{1}{(1-e^2)^2}.$$

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

In practise: 

We can derive \(\dot\omega\) as well as all terms in \(\dot\omega_{GR}\) and \(\dot\omega_{N}\) from observations, except for \(k_{2,1}\) and \(k_{2,2}\). 

We have an undetermined system... 

BUT: 

In the case of a twin binary system (two stars of same mass, radius, effective temperature), we can assume \(k_{2,1} = k_{2,2} \equiv k_2\).

In this case, we have \(\dot\omega = constant * k_2\). 

We can therefore put stringent constraints on the internal mass distribution.

The internal stellar structure constant of a star, also called the apsidal motion constant, is obtained from the relation 

$$k_2 = \frac{3-\eta_2(R_*)}{4+2\eta_2(R_*)},$$

where \(\eta_2(R_*)\) is the logarithmic derivative of the surface harmonic of the distorted star expressed in terms of the ellipticity \(\epsilon_2\) and evaluated at the stellar surface \(R_*\): 

$$\eta_2(R_*) = \left(\frac{d\ln \eta_2}{d \ln r}\right)_{r=R_*},$$

which is the solution of the Clairaut-Radau differential equation 

$$ r\frac{d\eta_2(r)}{dr}+\eta_2^2(r)-\eta_2(r)+6\frac{\rho(r)}{\bar\rho(r)}(\eta_2(r)+1)-6 = 0$$

with the boundary condition \(\eta_2(0)=0\) (Hejlesen 1987, A&AS, 69, 251).

In this expression, \(r\) is the current radius at which the equation is evaluated, \(\rho(r)\) is the density at a distance \(r\) from the centre of the star, and \(\bar\rho(r)\) is the mean density within the sphere of radius \(r\).

This can be better appreciated if we look at the evolution of \(k_2\) as a function of the hydrogen mass fraction \(X_c\) inside the core (left panel) and as a function of the radius of the star (right panel).

We here do not enter into much details regarding the different models shown, and we rather refer to Sections 3.3 and 4.1 for a description of these models. 

We computed stellar structure and evolution models with two different codes: Clés and GENEC. 

Clés:

Clés (Code Liégeois d'Evolution Stellaire) is developed and maintained at the University of Liège, Belgium (Scuflaire et al. 2008, Ap&SS, 316, 83).

Wind mass-loss is included in the models through the Vink et al. (2001, A&A, 369, 574) formalism. 

Overshooting is included in the models through the \(\alpha_{ov}\) parameter.

Turbulent transport is taken into account in the models and modelled as a diffusion process by adding a pure diffusion term of the form 

$$ D_T \frac{\partial \ln X_i}{\partial r},$$ 

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, ApJ, 529, 338). In this expression, \(r\) is the radius and \(X_i\) is the mass fraction of the element \(i\) at that point. 

The turbulent diffusion coefficient \(D_T\) (in \(cm^2\,s^{-1}\)) takes the following expression: 

$$D_T = D_{turb} \left(\frac{\rho}{\rho_0}\right)^n + D_{ct}, $$

where \(D_{turb}, n, D_{ct}\) are chosen by the user. 

The main advantage of Clés is its routine called min-Clés in which a Levenberg-Marquardt minimisation is implemented to search for best-fit models. min-Clés determines the combination of the Clés models' free parameters, including, notably, the initial mass and stellar age, which allow for the best reproduction of a set of the star's currently observed properties (mass, radius, effective temperature, apsidal motion constant). 

GENEC:

GENEC code is developed and maintained at the University of Geneva, Switzerland (Eggenberger et al. 2008, Ap&SS, 316, 43).

Wind mass-loss and overshooting are included in GENEC as in Clés.

The main advantage of GENEC is that it accounts for the two-dimensional stellar surface deformation due to rotation. In its binary version (Song et al. 2013, A&A, 556, A100), this code allows the inclusion of the effects of tidal interactions on the stellar structure, and thus on the apsidal motion constant. It accounts for the mixing induced by both the stellar rotation and the tidal interactions. 

We combined RV data and photometric data to derive a consistent apsidal motion rate for the system using the PHOEBE 2 code (Prsa et al. 2016, ApJS, 227, 29; Jones et al. 2020, ApJS, 247, 63).

We derived \(\dot\omega = 1.843 \pm 0.083 \, ^\circ\,yr^{-1}\). 

We computed stellar structure and evolution models with the min-Clés routine asking for the models to reproduce the mass, radius, and effective temperature of the stars. The best-fit models obtained in this way predicted stars having too high a \(k_2\) value, that is to say, stars having too homogeneous a density stratification between their cores and external layers, and hence too high an apsidal motion rate. 

If we ask for the models to also reproduce the \(k_2\) value of the stars and hence the apsidal motion rate of the system, the best-fit models have an enhanced mixing (measured through the overshooting parameter and the turbulent diffusion in the Clés models) compared to the "standard" models obtained previously. 

This can be seen in the following three panels showing the radius of the star as a function of its mass, its luminosity as a function of its effective temperature, and the apsidal motion rate of the system as a function of the internal structure constant of the stars. Models of increasing mixing are depicted from open circle to open square symbols.

It can also be visualised in the Hertzsprung-Russell diagram showing different evolutionary tracks with different values of the turbulent diffusion (left panel). The dots overplotted on the tracks correspond to the models having the correct value of \(k_2\). These models are found outside of the observational boxes of the stars, showing that the best-fil models have too high a \(k_2\) value (given that \(k_2\) decreases with time as the star evolves). 

As a preliminary investigation, we built stellar evolution models with the GENEC code accounting for both the stellar rotation and binarity. We show that rotation can be, at least partially, responsible for the enhanced mixing required in the Clés models (right panel). On the contrary, binarity does not seem to change the picture much. 

We derived the apsidal motion rate of the system based on the RVs: \(\dot\omega = 1.198 \pm 0.300\, ^\circ\,yr^{-1}\). 

We derived values for \(\omega\) for the lightcurves of the system and checked their consistency with the values of \(\omega\) expected from the RVs adjustement (left panel). 

We derived values for the phase differences between primary and secondary eclipses for the different photometric data and checked their consistency with the phase differences expected from the RVs adjustement (right panel). 

We searched for best-fit models for the primary and secondary stars separately. The need for enhanced mixing established for HD 152248 is confirmed here, especially for the primary star, as can be seen in the following HR diagram.

We derived the apsidal motion rate of the system based on the times of minima of the eclipses (left panel). We obtained \(\dot\omega = 15.38 \pm 0.51 ^\circ\,yr^{-1}\).

We adjusted the light curves of the system and checked that the values obtained for \(\omega\) are in agreement with those inferred from the times of minima (right panel). 

We derived the apsidal motion rate of the system based on the RV curve of the primary star (left panel). We obtained \(\dot\omega = 2.04 \pm 0.24 ^\circ\,yr^{-1}\). 

The fact that there is a single eclipse leads to a degeneracy between the radii of the primary and secondary stars and the orbital inclination that can shift the formally best-fit solution to values of the radii that are inconsistent with the spectroscopic brightness ratio. We used the spectroscopic brightness ratio to constrain the radii of the stars and derive an inclination for the system. Two TESS light curves of the system are shown (right panel). 

The massive eccentric binary HD 165052 is a well-known early-type double-lined spectroscopic binary system of the very young open cluster NGC 6530, located in the Lagoon Nebula. 

The low inclination of the system leads to the absence of photometric eclipses. 

We combined our RVs derived for a set of spectroscopic observations with RVs from the literature, with a total timespan of 46 years, to derive the rate of apsidal motion for the system: \(\dot\omega=11.30 \pm 0.64 \,^\circ yr^{-1}\).

We computed stellar structure and evolution tracks with Clés to derive evolutionary masses for the stars and put a constraint on the value of the inclination of the orbit through the confrontation with the minimum stellar masses derived from the RVs. 

We assumed that HD 165052 has an age of \(2.0 \pm 0.5 \) Myr as derived in the literature. 

Assuming that the primary model with \(M_{init} = 25 M_\odot\) and \(D_T =  2\times 10^7 cm^2 s^{-1}\) and the secondary model with \(M_{init} = 21 M_\odot\) and \(D_T =  2\times 10^7 cm^2 s^{-1}\), both at ages 2 Myr, are representative of the stars, we infer evolutionary masses \(M_{ev,P} = 24.8 \pm 1.0 M_\odot\) and \(M_{ev,S} = 20.9 \pm 1.0 M_\odot\) and evolutionary radii \(R_{ev,P} = 7.0 \pm 0.5 R_\odot\) and \(R_{ev,S} = 6.2 \pm 0.4 R_\odot\).

The evolutionary masses put a constraint on the inclination of the orbit: \(22.1^\circ \le i \le 23.3^\circ\). This small value of the inclination confirms that eclipses are very unlickely to be seen in photometric observations of the binary. 

The combined analysis of observational data (both spectroscopic and photometric observations) of (massive) eclipsing and spectroscopic binary systems allows us to derive

1. The fundamental (global) properties of the stars: effective temperatures, radii, masses, and luminosities;
2. The properties of the binary system: eccentricity, inclination, orbital period, semi-major axis, apsidal motion rate. 

Through the apsidal motion equations and their dependence on the internal stellar structure of the stars, the apsidal motion rate allows us to probe the interior of the stars. In particular, this internal stellar structure constant \(k_2\) gives us an indication of how the mass is distributed inside the star, in other words, of the density stratification inside the star. 

With the fundamental stellar parameters, together with the \(k_2\)-parameter, we can build dedicated stellar structure and evolution models to constrain the internal mixing occurring inside the stars. 

In this project, we have demonstrated that standard stellar evolution models predict stars that have too low a density constrast between their core and external layers compared to what is expected from the observations. To reproduce the properties of the stars, and in particular the apsidal motion rate, enhanced mixing needs to be included inside the models. 

This conclusion is a very important discovery that paves the way for the future stellar evolution modelling and demonstrates why and how the study of apsidal motion is so precious in this context. Indeed, without the knowledge of the apsidal motion rate, we would never have identified that the standard models underestimate the mixing inside the models. 

The enhanced mixing being entirely due to stellar rotation is still under investigation. However, ongoing research with the MoBiDICT non-perturbative code shows that this need for enhanced mixing is not inherent to the 1D stellar modelling and is rather even exacerbated when the non-perturbative formalism is adopted in the 3D stellar modelling, especially in close eccentric binaries (see figure).   

• From the observational side, we are currently increasing the number of systems for which an accurate apsidal motion rate is derived together with the fundamental stellar properties and orbital parameters like the analyses performed in this project.

• From a stellar evolution side, we need to confront more systems with observationally determined properties to stellar evolution models to constrain the mixing processes in the models. 

• We plan to investigate the origin of this enhanced mixing: in which proportion is that due to stellar rotation or tidal interactions? 

• Finally, at this stage, we cannot support more the need for 3D hydrodynamic simulations to capture the proper deformations of the stars, as well as the need for including the complete complex binary interations in stellar evolution models but also in observational codes used to fit data such as light curves. 

The need for enhanced mixing inside the models to reproduce the stellar properties of massive stars is supported by 2D hydrodynamical simulations. 

Simulations for the binaries HD 152218, HD 152219, and CPD-41° 7742 show that high values of the overshooting parameter (in red) are always necessary to reproduce the position of the stars in the HR-diagram (Baraffe et al. 2023, MNRAS, 519, 5333). 

The same study shows that models without overshooting or with a small overshooting are unable to reproduce the main sequence width, especially for the more massive stars. Additional mixing is indeed necessary.

These results are in agreement with and support the works of Castro et al. 2014, A&A, 570, L13 and Martinet et al. 2021, A&A, 648, A126. 

To answer this question, we used the MoBiDICT code developed and maintained at the University of Liège, Belgium (Fellay and Dupret 2023, A&A, 676, A22).

MobiDICT is a non-perturbative method. 

MoBiDICT stands for Modelling binaries Deformations Induced by Centrifugal and Tidal forces. 

Main differences between perturbative and non-perturbative models

Perturbative models 

• Centrifugal and tidal forces = small perturbations of spherical symmetry only accounting for leading terms;
• Significantly underestimate the deformation of binaries, like Roche models do too; 
• Dipolar term of the gravitational potential is neglected, while it is of the same order of magnitude as the leading tidal term in most distorted cases! 

MoBiDICT: non-perturbative model

• Computes the entire precise 3D structure of each component;
• Highest deformation observed for stars having significant envelope mass; 
• Calculates instantaneous tidal acceleration perturbation and its consequence on the apsidal motion. 

We computed models for low-mass (green), solar (orange), and massive (blue) stars on the main-sequence phase as well as for a red giant branch star (purple). We here explicitly assume the stars belong to twin binaries.

The impact on \(\eta_2\) is important for small orbital separation but is almost independent of the type of star (left panel). 

It has a direct impact on the apsidal motion computation. 

These results show that the perturbative model assumption is not justified in high distorsion cases, when the semi-major axis is smaller than 5 times the radius of the primary star (assumed to be the most evolved star, right panel). 

We computed the theoretical apsidal motion rate for twin binaries assuming an eccentricity of 0.1. We show the difference in apsidal motion rate between MoBiDICT and the perturbative case scaled by the apsidal motion rate obtained in the perturbative case (left panel), as a function of the semi-major axis of the system.

We observe that the impact on the apsidal motion rate is, as expected, more important for close binary systems and for stars of highest envelope mass. 

We also computed the theoretical apsidal motion rate for a twin binary system made of massive stars and assuming different values of the eccentricity (right panel).

We observe that the impact on the apsidal motion rate is more important for higher eccentricities. This is expected given that a higher eccentricity means the stars get closer during periastron passage. Hence, the tidal interactions are more important.

We conclude that the perturbative models underestimate the apsidal motion rate when \(a(1-e^2)/R_1 < 7\). 

We compare the results obtained for perturbative models and MoBiDICT for four twin binary systems (see Table).

We show that instead of solving this "enhanced mixing issue", the MoBiDICT models require even more mixing than the perturbative models to reproduce the observations. 

In addition, we show that the non-perturbative correction is on the order of the General Relativistic correction. Hence, we definitely cannot ignore the non-perturbative term.