© The Authors 2024

1. Introduction

Tidal forces in binary systems result from the gravitational interaction between non-point-like bodies and can be categorised into two components: the equilibrium tides, corresponding to large-scale circulations resulting from a hydrostatic readjustment (i.e. stellar deformations) and the dynamical tides corresponding to the excitation of eigenmodes of oscillations from the periodic perturbations from a companion (Zahn 1975). These tidal forces and the subsequent stellar deformations lead to various phenomena such as orbital synchronisation, mass transfer, and tidal dissipation (Jeans 1929), the latter corresponding to the dissipation of orbital energy within the stellar bodies. A direct consequence of the tidal interactions is the apsidal motion in eccentric binaries (Sterne 1939). It represents the time variation of the argument of periastron ω or, in other words, the motion of the apsidal line with time. Apsidal motion originates from equilibrium tides, dynamical tides (Gimenez & Margrave 1985; Willems & Claret 2002), and general relativistic corrections¹ (Gimenez & Margrave 1985). It serves as a valuable observational constraint to understand the structure of deformed stars, providing insights into the microphysics of stellar models (Rosu et al. 2020a, 2022a,b, 2023). Accurate modelling of tidal and centrifugal deformations is essential to determine the precise stellar structure of each component in binary systems and, consequently, to predict the apsidal motion.

When modelling the structure and deformations of binary stars, the perturbative approach is one of the most sophisticated and widely used methods in the literature. The principle of this methodology is to consider that stellar deformations are small perturbations to the spherical symmetry of stars, only accounting for the leading orders terms in the developments (Sterne 1939; Kopal 1959, 1978). This approach allows one to obtain simplified structural deformations from the unperturbed spherically symmetric structure through the well-known stellar structure constant, the k₂ coefficient that is obtained from solving the Clairaut-Radau equation (Kopal 1959). With the perturbative formalism, the dipolar components of the tidal and centrifugal forces cancel out (Kopal 1959; Fitzpatrick 2012), implying that no asymmetry exists in the stellar models, in contradiction with the fundamental Roche geometry. Moreover, our previous work (Fellay & Dupret 2023) already demonstrated that the perturbative approach underestimates the deformations when compared to our non-perturbative models. Therefore, it becomes crucial to verify the consequences of the assumptions inherent to the perturbative approach. Finally, it is worth mentioning that the tidal distortion and dissipation formalism (Zahn 1977; Hut 1981) employed by widely recognised and publicly available codes such as MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019) are based on the perturbative approach.

In this study, we expand on the work done in Fellay & Dupret (2023) by further investigating the limitations imposed by the assumptions of the perturbative approach, using our non-perturbative method. Our tool, called Modelling Binary Deformations Induced by the Centrifugal and Tidal forces (MoBiDICT), allows us to determine the deformed structure of close binaries in a non-perturbative manner, accounting for the redistribution of mass in the bodies. By utilising the deformed structure obtained from MoBiDICT, we could calculate the instantaneous tidal acceleration perturbation and its consequence on the apsidal motion of binary systems. Finally, MoBiDICT’s formalism enabled us to locate and identify the specific spherical order responsible for the discrepancies seen.

In this article, we start, in Sect. 2, by introducing the physics behind MoBiDICT and the methodology developed to compute the apsidal motion of binaries in a non-perturbative way. In Sect. 3, we explain how we developed the perturbative approach formalism usually used to model the apsidal motion of binaries. Sect. 4 is dedicated to exploring the discrepancies in models for different types of theoretical stars. In Sect. 5 we study binary systems for which the apsidal motion has been measured and provided in the literature. In Sect. 6, we discuss the implication of our findings on the stellar modelling methods used for binary stars. Finally, in the Appendices, we develop all the equations of the perturbative approach starting from a single star’s gravitational potential to arrive at the apsidal motion in binaries.

2. Non-perturbative modelling

2.1. Set-up of the problem and first models

In this article, we investigate the deformations in binary systems resulting from three primary sources: the gravitational force, the tidal force arising from a companion, and the centrifugal force originating from the orbital rotation of the system and the individual stellar rotations. Accurate modelling of these phenomena necessitates a non-perturbative treatment of deformations in three dimensions (3D).

We considered a binary system consisting of two distinct stars separated by a distance r and with an orbital rotation rate n. Assuming the stars (i = 1, 2) within the system are in hydrostatic equilibrium, their structures are governed by the equation of hydrostatic equilibrium:

(1)

and the Poisson equation,

(2)

linking the gravitational potential Ψ_i(r_i, μ, ϕ) of a body to its density distribution ρ_i(r_i, μ, ϕ) in 3D. In Eq. (1), P is the pressure, g_eff is the effective gravity, Ψ_c(r_i, μ, ϕ) is the centrifugal potential on a 3D grid, where r_i denotes the radial coordinate of the star i, μ = cos(θ) its colatitude, and ϕ its azimuthal angle. The centrifugal force can be solely derived from a potential if the rotations exhibit cylindrical symmetry, with the solid body rotation being a specific case. Throughout this entire article, we made the assumption of solid body rotations to simplify the problem (see Sect. 6 for a discussion on the validity of this assumption).

In our recent work (Fellay & Dupret 2023), we introduced a new tool called MoBiDICT, specifically designed to precisely model close binaries using a non-perturbative approach to treat the deformations. Our method iteratively solves the Poisson equation on a 3D structure, allowing to obtain the gravitational and tidal potential accounting for the redistribution of mass within the bodies. MoBiDICT takes advantage of the conservative nature of all the forces considered, resulting in a barotropic model structure (i.e. the densities are constant on equipotentials). By assuming that in a given direction (μ_crit, ϕ_crit) each star is described by a 1D spherically symmetric input model, we can recompose the entire 3D density profile from the total potential (Ψ_tot = Ψ₁ + Ψ₂ + Ψ_c) of each star. This method is thoroughly explained in Fellay & Dupret (2023).

Figure 1 illustrates the equipotential lines within a deformed primary star, that were obtained using our non-perturbative modelling method. In this case, we applied our tool to a binary system composed of a primary 1 M_⊙ main-sequence (MS) star and a 0.2 M_⊙ MS companion separated by r = 2R₁, R₁ being the radius of the primary given by the 1D spherically symmetric model.

In the scenario depicted in Fig. 1, we assumed that the rotational velocities of both stars are identical and equal to the orbital rotation rate, in other terms the system is tidally locked. However, when studying eccentric binaries, pseudo-synchronisation becomes the end result of tidal interactions, making impossible to model the system as entirely synchronised. Furthermore, mass transfers between binary components often result in significant transfer of angular momentum (Packet 1981), leaving the accretor star close to its critical rotation while the donor slowly rotates. Generalising our model to non-synchronised binary systems is necessary to model eccentric binaries and binary evolution.

2.2. Non-perturbative modelling of non-synchronised binary systems

We considered a scenario where each star i composing a system has its individual rotation rate Ω_⋆, i, independent from the orbital rotation rate n. In this situation, assuming alignment of all rotational axes, two planes of symmetry emerge and are exploited by our method. The first plane of symmetry is the orbital plane, the second is the plane including all the rotation axes. A 3D representation of this configuration is shown in Fig. 2 for the same system as in Fig. 1.

With our non-perturbative method, unsynchronised binary systems can be directly modelled by modifying, and generalising the centrifugal potential. The previous expression of the centrifugal potential, which relied on the assumption of solid body synchronised rotation and the circular orbit was

(3)

and becomes in the non-synchronised cylindrical symetric rotation and circular orbit case:

(4)

where s is the distance from the stellar rotation axis. In Eqs. (3) and (4), we used Cartesian coordinates centred on the star considered. Here, x represents the coordinate along the apsidal line (i.e. the major-axis of the orbit) pointing towards the companion, x_CM denotes the distance between the system’s centre of mass and the star’s centre, and y is the coordinate along the perpendicular axis to x included in the orbital plane.

As expected from the comparison of Eqs. (3) and (4), for a synchronised system (Ω_⋆, i = n), these expressions are equivalent up to a constant factor. The effective gravity resulting from the centrifugal force can be obtained by differentiating Eq. (4).

Our model can be generalised to eccentric systems by decomposing the motion of each element as a Keplerian motion around the system’s centre of mass plus a cylindrical rotation of the star. To obtain the centrifugal potential at a given point of the orbit where the instantaneous distance between the centre of the two binary components is r, n² has to be replaced, in Eq. (4), by ñ² defined as

(5)

so that ñ^2r is the norm of the Keplerian instantaneous relative acceleration of the bodies in the system. In Eq. (5), M_tot is the total mass of the system. With this decomposition and our assumption of the cylindrical rotation, the acceleration of an element in a star is entirely given by the gradient of a potential. The acceleration induced from our decomposition of the mouvement is the gradient of the centrifugal potential that is expressed for the primary star as

(6)

where the bold symbols denote vectorial quantities, is the unit vector linked to the x-axis of the primary and is the unit vector linked to the gradient of s. In all this article, we symplify further the problem by only consiering solid body rotations.

2.3. Force perturbation

In classical celestial mechanics theory, the relative acceleration of two bodies in a binary system is entirely given by its radial component expressed as

(7)

where the bold symbols denote vectorial quantities, F₂₁ is the force exerted by the secondary star on the primary (noting F_2 → 1 ≡ F₂₁), F₁₂ its opposite, is the radial unit vector centred on the secondary and pointing towards the primary, M_i are the individual masses of the bodies and m_μ is the reduced mass of the system defined as

(8)

The last step of Eq. (7) was deduced from the third Newton Law (F₁₂ = −F₂₁). The total force per unit of mass exerted by the secondary on its primary can be decomposed as

(9)

In Eq. (9), F_R is the radial perturbed acceleration resulting from considering stars as non-spherical extended interacting bodies. This perturbation is responsible for all the equilibrium tidal-induced dynamical effects of binary evolution, including orbital migration (i.e. variation of systems semi-major axis), eccentricity dissipation (i.e. reduction of orbital eccentricty), and apsidal motion (i.e. time variation of the argument of the periastron).

In the framework of our modelling, the force exerted by the secondary on the primary can be expressed as

(10)

where F₂ is the gravitational force from the secondary inside the primary. By decomposing the densities and the forces in unperturbed spherical symmetric terms, respectively ρ₁₀ and F₂₀, and non-spherical perturbations, respectively and , the force exerted by the secondary on the primary can be rewritten as

(11)

The spherically symmetric component of the force can be simply identified as being

(12)

Therefore, the perturbed acceleration is expressed as

(13)

In MoBiDICT, all quantities are projected on spherical harmonics to facilitate the solving of the Poisson equation (Eq. (2)), for example the densities and potentials are decomposed as

(14)

and

(15)

where L is the free parameter of our modelling defining the maximum degree of the spherical harmonics to be accounted for, p = 0 if ℓ is even, p = 1 otherwise, and m = 2k + p. In this formalism, the decomposition in spherically symmetric and non-symmetric terms is straightforward. The ℓ = 0,  m = 0 term is the spherically symmetric component of each quantity while all other spherical orders are non-symmetric terms. In our formalism the perturbed acceleration can thus be obtained with

(16)

Finally, each force term projected on the spherical harmonic basis can be decomposed and re-expressed in a spherical coordinate frame (to be integrated) with

(17)

The spherical harmonics Y_ℓ, m in Eq. (17) are functions of μ₂ and ϕ₂. The expressions of the , and terms can be found in Fellay & Dupret (2023). More than having a non-perturbative description of the resulting force, the main advantage of our formalism is the possibility of naturally compare the contributions of the different spherical harmonic orders by selecting the desired orders in the sums.

2.4. Non-perturbative apsidal motion computation

In this section, we developed the formalism used to model the apsidal motion with our non-perturbative method. Figure 3 illustrates the setup of an eccentric binary system, the apsidal motion of a system corresponding to dω/dt. In eccentric binary systems, aspidal motion is a direct observational constraint on the deformations of distorted bodies. Recent investigations (Rosu et al. 2020a, 2022a,b) highlight its utility in constraining physical processes occurring during stellar evolution, such as internal mixing.

The apsidal motion is induced by the equilibrium tidal perturbed acceleration F_R, the general relativity and dynamical tidal perturbed acceleration. The general relativistic component of the apsidal motion can be expressed as

(18)

where a is the system’s semi-major axis, e its eccentricity, P_orb its orbital period, c is the celerity, and G the gravitational constant. In this article, we studied the equilibrium tides component of the apsidal motion and neglected the dynamical tides component, the latter will be the focus of a forthcoming article.

From the Gauss planetary equations for a perturbed acceleration along the apsidal line (see e.g. Kopal 1959), the apsidal motion of a binary system is classically given by

(19)

where φ is the true anomaly and ⟨F_R cos φ⟩ denotes the time average of the perturbed acceleration expressed as

(20)

To obtain the apsidal motion of a system through a non-perturbative approach, the sole quantity to be computed by our non-perturbative modelling is the time average perturbed acceleration. To achieve this, we took advantage of the fast computing capabilities of our method, we captured snapshots of stellar deformations at various points along the binary orbits, as depicted in Fig. 3. Using these snapshots, we can derive the time-averaged perturbed acceleration by reformulating its expression with the second Kepler’s law:

(21)

where r is the separation between the stars expressed, following the first Kepler’s Law, as

(22)

By combining Eqs. (19)–(22), the apsidal motion obtained with our non-perturbative modelling is given by

(23)

where the integrals are only taken over half of an orbit by symmetry. In practice, with our method’s implementation, merely fifteen equally spaced values of φ over half of an orbit are sufficient to reach a relative precision of 10⁻³ − 10⁻⁴ on the apsidal motion in the vast majority of close binaries.

3. Perturbative model

3.1. The Clairaut-Radau equation

Perturbative models are obtained by assuming that the deformations of the stellar models are small (more precisely, it is assumed that in Eq. (A.1)), retaining only first order terms in derivations. Similar to our modelling approach, the gravitational potential of each star i is obtained through the solution of the Poisson equation projected on a spherical harmonics basis,

(24)

The main distinction in this modelling is to treat stars as lightly perturbed spherically symmetric bodies. This assumption greatly simplifies the problem and allows us to relate the deformations of a body directly to its unperturbed structure and system properties. In this framework, the deformations of each body are linked to its structural coefficients η_ℓ (see Eqs. (A.1) and (A.20)). These coefficients can be obtained by solving the Clairaut-Radau first-order differential equation Kopal (1959), and derived in Appendix A) expressed as

(25)

using the boundary conditions

(26)

In the Clairaut-Radau equation, r̅ is the average radius of a chosen isobar while is the average density under the chosen isobar. The relationship between the deformed structure and the η_ℓ coefficients is explained and derived in Appendix A.

3.2. Perturbation of the gravitational potential and force by a companion

We considered the perturbation induced by a point-mass secondary on the primary. As developed in Appendix B, the potential perturbation caused by the secondary at the surface of the primary is given by

(27)

where is the (ℓ, m) component of the gravitational potential generated by the secondary evaluated at the surface of the primary, and the centrifugal components of the potential are noted . The first terms of each of these quantities can be found in Appendix C. The coefficients k_ℓ, 1 also known as Love numbers, characterise the surface non-spheroidal deformation response of a body subjected to a perturbative potential. Mathematically, these coefficients are defined as:

(28)

Equation (27) provides a direct comparison between the models, where the quantity is equivalent to of our non-perturbative iterative method. One particular application to Eq. (27) concerns the dominant contribution to deformations and perturbations, the quadrupolar terms (ℓ = 2 terms). For those terms, the sum of the centrifugal and tidal potentials generated at the surface of the primary by the secondary treated as a point-mass body is expressed, in the non-synchronised case, as

(29)

and

(30)

In both synchronised and non-synchronised systems, the dipolar component of the total potential is null when the secondary is treated as a point-mass, as gravitational and centrifugal forces cancel out according to the third Kepler law. As highlighted earlier, when computing the apsidal motion of a system, it is imperative to account for the tidally induced acceleration perturbation F_R. By only considering the spherical order lower than ℓ = 3, the perturbed acceleration can be decomposed as follows:

(31)

where F_2′10 is the force exerted by the quadrupolar deformation of the secondary on the unperturbed primary, F_201′ is the force exerted by the unperturbed secondary on the already perturbed primary, and F_2′1′ is the force exerted by the perturbed secondary on the already perturbed primary. The last term, F_2′1′, is a second order term, thus neglected in the pertubative approach. In accordance with Newton’s third law, the force exerted by the unperturbed primary on the already perturbed secondary is the opposite of the force exerted by the perturbed secondary on the unperturbed primary, mathematically F_201′ = −F_1′20. As shown in Appendix D, the force exerted by the perturbed secondary on the unperturbed primary is given by

(32)

similarly, F_1′20 can be obtained by permuting the primary and secondary in Eq. (32). Finally, combining Eqs. (31) and (32) the perturbed acceleration is expressed as

(33)

with q = M₂/M₁.

3.3. Perturbative apsidal motion computation

In the case of an eccentric orbit, the mean orbital rotation rate n is distinct from the instantaneous orbital rotation rate ñ experienced by each individual binary component. Using this formalism, the perturbed acceleration can be rewritten by grouping the tidal and centrifugal terms as

(34)

where P_⋆, 1 and P_⋆, 2 are respectively the rotation period of the primary and secondary, and R_tides and R_c can be identified in Eq. (34). The apsidal motion of a system is dependent on the time average of the perturbed acceleration over an orbit ⟨F_R cos φ⟩. Using Eq. (20), the perturbed acceleration averaged over an orbit is given by

(35)

The integrals in Eq. (35) can be re-expressed with the zero-order Hansen coefficients defined for every positive m as

(36)

Without entering into the detailed expressions of these coefficients, the two orders that are of interest in our case are the n = −4, m = 1 and n = −7, m = 1 that are given by

(37)

and,

(38)

The average perturbed acceleration is finally given by

(39)

and can be inserted in Eq. (19) to obtain the apsidal motion of a binary system with the perturbative approach:

(40)

with,

(41)

and,

(42)

4. General model comparisons

This section is dedicated to a theoretical comparison of different types of stars encompassing a range of potential binary systems across different mass regimes. Here, we present the input stellar models used to compare the perturbative and non-perturbative modelling approaches. We chose four models to represent a broad spectrum of stars and structural characteristics. The physical ingredients of the chosen models are detailed in Table 1, those models are the same models compared in Fellay & Dupret (2023).

Table 1 in Fellay & Dupret (2023) has only been extended to include the k₂ value of each model. As shown by Eq. (27), to a higher k₂ corresponds a greater surface response from a body when subjected to a potential perturbation to its surface. Moreover, k₂ can be used as an indicator of a body’s deformability. In addition, our results from Fellay & Dupret (2023), showed that the low-mass and red giant branch (RGB) stars are the most distorted by a companion while the massive and solar type stars are expected to be less deformed and impacted by our non-perturbative method for a given a/R₁ ratio. In the following sections, we compared the surface deformations and apsidal motion obtained with our non-perturbative modelling method and the perturbative modelling for twin binary systems (q = 1) composed of the stars presented in Table 1. Twin binary systems being systems composed of two identical stars (i.e. stars having the same mass, radius, effective temperature, and rotational period).

4.1. Potential and forces discrepancy

The perturbation of the surface potential is a reliable indicator of the deformations undergone by a star. Within either our formalism or the perturbative formalism, the perturbation of the surface potential is given by all spectral terms of the surface potential with ℓ > 0. In the perturbative method, this quantity is given by Eq. (27), which is applied to the ℓ = 2 terms in Eqs. (29) and (30). Within our non-perturbative method, are obtained in an iterative process, solving Poisson’s equation at each step and for each spherical order. Consequently, there is no analytical expressions available for this quantity in our method. To compare the surface deformations with different modelling procedures we introduce the discrepancy of surface potential, denoted and defined as:

(43)

where is the spectral surface potential obtained with our non-perturbative method and is the same quantity obtained with the perturbative approach. In Fig. 4, we show this surface potential discrepancy as a function of the orbital separation for twin binary systems composed of the stars presented in Table 1. Our focus in this analysis is on the dominant terms, the ℓ = 2 terms.

The surface potential discrepancy, varies significantly depending on the stellar type and the orbital separation. For both ℓ = 2 components, stars with the higher k₂ exhibit more pronounced effects from our non-perturbative treatment, resulting in differences of up to ∼40% in the ℓ = 2, m = 2 term for low mass stars. This implies that the higher the deformations, the more substantial the underestimation of the surface deformations by the perturbative approach. Furthermore, orbital separation plays a critical role, with closer star configurations yielding stronger deformations and more notable disparities in the surface potential.

By comparing the two panels in Fig. 4, we observe a significantly smaller difference in surface potential for the ℓ = 2, m = 0 term, regardless of the model. The ℓ = 2, m = 0 component arises from both centrifugal and tidal deformations, while the ℓ = 2, m = 2 term is solely originating from the tidal deformation. Both our modelling and perturbative modelling are sharing the same assumptions on the centrifugal potential, therefore, the difference originating from tidal contribution in the ℓ = 2, m = 0 term is diluted by the centrifugal contribution. The major difference between the perturbative and non-perturbative modelling approaches is thus appearing in the treatment of the tidal deformation.

4.2. Apsidal motion

In Sect. 2.4 we showed that the apsidal motion of an eccentric binary system is proportional to the time average of the perturbed acceleration. As detailed in Sect. 4.1, the perturbed acceleration denotes the tidal acceleration originating from a nearby companion. A modification of the surface ℓ = 2 potential as illustrated in Fig. 4 introduces variations in the perturbed acceleration. In this section, we investigated the perturbed acceleration discrepancy resulting from the potential differences seen and the consequent impact on the apsidal motion of eccentric binaries.

In the classical perturbative model, the perturbed acceleration F_R is solely originating from the quadrupolar component of the deformation. The detailed expression of F_R is given in Eq. (33) and derived in Appendix D. In opposition, our methodology incorporates a non-perturbative treatment of all spherical orders spanning ℓ = 1, 2, ...L. The exact expression of F_R is given by Eq. (16). To quantify the disparities between perturbative and non-perturbative perturbed accelerations, we define their difference ΔF_R as follows:

(44)

The difference of perturbed accelerations as a function of the separation between the twin binary systems components presented in Table 1 is illustrated in Fig. 5.

The perturbed acceleration discrepancy observed in Fig. 5 is in agreement with the results showed in Fig. 4. Models exhibiting the greatest surface potential discrepancies are also exhibiting notable perturbed acceleration differences. Moreover, models with higher k₂, namely the low mass and RGB stars, are the most impacted by our modelling. Solar-like stars and massive stars are comparatively less affected, even if in close orbit F_R discrepancies can still reach up to, respectively, 20% and 12%.

To assess the discrepancy of apsidal motion resulting from the perturbed acceleration differences, we need to consider binaries with eccentric orbits. To limit the free parameters of the system we compare binary systems with synchronised rotations and a fixed eccentricity, e = 0.1. The impact of the eccentricity on the apsidal motion is explored for systems in Sect. 4.3. We introduce the apsidal motion difference:

(45)

to compare the apsidal motion originating from the non-perturbative and perturbative approaches. In the perturbative approach, the apsidal motion is given by Eq. (40). With our non-perturbative method we required to capture snapshots of the binary system along half of an orbit to estimate the perturbed acceleration time average and hence the apsidal motion (see Sect. 2.4). For each binary system composed of the stars presented in Table 1, Fig. 6 illustrates as the function of the binaries orbital separation normalised by the stellar radii.

The apsidal motion discrepancy is a direct result of the perturbed acceleration illustrated in Fig. 5. Binary systems with greater F_R discrepancis exhibit the highest apsidal motion differences. For the low mass and RGB stars, the apsidal motion discrepancy can reach up to respectively 70% and 45% when the stars are close to contact, in their periastron. Similarly, for the solar type and massive stars we found an apsidal motion discrepancy up to respectively 30% and 15%. These discrepancies vary depending on the exact structure of the stars or the architecture and the orbital parameters of the system.

4.3. Dependency on the orbital eccentricticty

In the previous section, we fixed the eccentricity of the systems to infer the modelled apsidal motion discrepancy. In this section, we study the impact of varying the orbital eccentricity on the apsidal motion discrepancy. We expect that binaries with higher eccentricities are closer in their periastron, consequently more deformed and more impacted by our methodology. We focussed on the twin binary system composed of 20 M_⊙ stars as this system is less impacted by our modelling and will consequently give a lower limit to the modelling discrepancies for a given semi-major axis. We computed a grid of models with different eccentricities and orbital separations, looking at the apsidal motion relative difference. The results from this computation are illustrated in Fig. 7². The eccentricity is indeed a crucial parameter directly impacting the modelling discrepancies. As expected, binary systems with higher eccentricities are more significantly impacted by our modelling method for a given orbital separation. For our least deformed model, we found that the modelling discrepancy reaches at least 5% when a(1 − e²)/R₁ ≲ 4.5 and 2% when a(1 − e²)/R₁ ≲ 6.5. For all our other theoretical models these thresholds are higher, by applying the same methodology to the twin system composed of 0.2 M_⊙ stars, we found that the modelling discrepancy thresholds of 5% and 2% are respectively located at a(1 − e²)/R₁ ≲ 6.5 and a(1 − e²)/R₁ ≲ 4.5.

In addition to this work, we verified the precision of our method at higher orbital separation and found that for a(1 − e²)/R₁ = 20 the discrepancies between the models reach 0.1% indicating a slow decrease of the apsidal motion discrepancy at high orbital separation and a sufficient precision to impose the thresholds given in this section.

4.4. Dependency on k₂ and the orbital separation

In the previous sections, we saw that a hierarchy exists between the perturbed acceleration of different stellar models and the resulting apsidal motion discrepancy. We observed that stellar models with higher k₂ values were more affected by our modelling approach. Another crucial quantity impacting the discrepancies is the separation of the binary components. In this section, we aim to comprehensively assess modelling discrepancies across a diverse range of stellar models and orbital separation.

To explore this, we constructed a grid of models of MS stars with masses ranging from 0.2 M_⊙ to 30 M_⊙ using the Code Liégeois d’Evolution Stellaire (CLES, Scuflaire et al. 2008). For each mass in our grid we took several models along the MS to have diverse models and stellar structures with similar k₂ to explore whether discrepancies are solely dependent on k₂ and the orbital separation.

Figure 8, illustrates the evolution of the acceleration perturbation relative difference ΔF_R/F_R, pert between our modelling and the perturbative approach. The top panel of the figure represents a group of stars that are either fully convective or composed of a radiative core and convective envelope, specifically stars with masses below 1.25 M_⊙. The lower panel corresponds to stars with a convective core, which includes stars with masses greater than 1.5 M_⊙.

Figure 8 reveals a distinct behaviour between the two groups of stars. While the lower k₂ limit of the low mass stars and the upper limit of the intermediate-high mass stars have similar k₂, their acceleration perturbation discrepancy is not similar. This result indicates that discrepancies between the models are not only a function of k₂ but also involve a component related to the intrinsic stellar structure. From our model grid, we observed that stars with different structures, particularly those with radiative or convective cores, displayed significant variations in discrepancy despite similar k₂ values. This phenomenon is probably related to the differences in density profiles between stars having or not a convective core. If looking for a possible parametrisation solely dependent on k₂ and orbital separation to correct the modelling discrepancies in the perturbative approach, the non-linear behaviours of the models greatly reduces its feasibility. We attempted an MCMC analysis using linear combinations of power laws of the orbital separation, for both the intermediate-high mass stars and low-mass stars. However, the results were unsatisfactory to propose an empirical parametrisation. As detailed in Sect. 5.2, the modelling discrepancies indeed arise from a combination of discrepancies in the ℓ = 1, 2, 3 terms. Non-perturbative modelling is necessary to recover the discrepancies seen in this section. However, Fig. 8 can be used to get an order of magnitude of the underestimation of tidal forces introduced by the perturbative modelling. The updated grid of k₂ values provided by Claret (2023) could for example be used to easily locate observed stellar systems in this diagram.

4.5. Dependency on the mass ratio q

Our analysis so far focussed on the dependency on R₁/a as this latter parameter has the most important impact on the results. However, the perturbative apsidal motion has a direct dependency on 8 parameters: (n, (n/Ω₁)², (n/Ω₂)², (R₁/a)⁵, (R₂/a)⁵, q, e, k₂₁, k₂₂), see Eq. (40). With our non-pertubative modelling the dependencies on k₂₁, k₂₂ become dependencies on the entire density profile of each star. Characterising the dependency of our result on mass ratio q is useful as this quantity can strongly impact the results and is a direct sub-product of the observations of binaries. To study and modify the q parameter in our modelling, different stellar models with the same age and initial chemical properties have to be used for the secondary star, impacting 4 of the parameters of our modelling for a given orbital separation: (n, (R₂/a)⁵, q, k₂₂ or ρ₂(r)). Consequently the effect of a change of q will be seen indirectly, through, in particular (R₂/a)⁵. For a given q parameter different than one, the behaviour of the apsidal motion discrepancy as a function of the orbital separation is dependent on the regime of (R₂/R₁)⁵, as discussed hereafter.

4.5.1.

This regime describes the case of a system composed of either a star plus a compact object or an evolved star plus a MS star. To simulate this situation, we chose to model the primary star as being the 1.5 M_⊙ RGB star described in Table 1. The secondary is a MS star with a lower mass, ranging between 0.2 and 1.25 M_⊙, the same age (2 Byr) and initial chemical composition (X₀ = 0.72, Z₀ = 0.010). In Fig. 9, we show the evolution of the apsidal motion discrepancy for the different systems considered. Based on Eq. (40), we chose the x-axis of Fig. 9 to be q^−1/5R₁/a to maintain the scale between (R₁/a)⁵ and q.

We can see that all the systems are superposed in Fig. 9 which indicates a clear proportionality of the apsidal motion discrepancy on q. When , as seen in Eq. (40), the terms with (R₁/a)⁵ are dominating for a given orbital separation, and the contribution from the secondary is negligible. In this situation, only the deformations of the primary impact the apsidal motion. As both the deformations and the apsidal motion are proportional to q for a given orbital separation, the apsidal motion discrepancy is also proportional to q.

4.5.2.

This regime describes systems composed of stars in similar evolutionary stages and structural properties. To simulate this situation we combined a primary MS star with a mass of 1.5 M_⊙, an age of 1.5 Byr, an initial chemical composition of X₀ = 0.72, Z₀ = 0.010, with a secondary star with the same properties but a lower or equal mass. In Fig. 10, we explore the apsidal motion discrepancy as a function of q^−1/5R₁/a for these binary systems.

In the case where the primary and secondary have similar radii despite different q, the properties of the secondary have an impact on the apsidal motion discrepancy, as illustrated by comparing Figs. 9 and 10. The apsidal motion discrepancy is still dominated by the contribution from the primary, consequently a partial proportionality to q can be seen. However, additional contributions from the other parameters of the secondary, in particular its density profile, have to be accounted for.

In our determination of the limit after which the modelling discrepancy reaches 2% and 5% (see Sect. 4.3), we assumed that q = 1. In practice, observed systems often have q < 1, and, to include the dependency of these limits on q we modified them. As for both cases presented above the dependency on q dominates the apsidal motion discrepancy, we need that q^−1/5a(1 − e²)/R₁ ≲ 4.5 to have at least 5% modelling discrepancy and q^−1/5a(1 − e²)/R₁ ≲ 6.5 to have at least 2% modelling discrepancy. By adopting these modifications we introduced a little bias in the regime where however, as showed before, the dependency on q still dominates.

5. Application to observed systems

5.1. Modelling of observed binaries and quantification of the extra mixing

In this section, we apply our non-perturbative modelling technique to existing twin binary systems, comparing the stellar parameters obtained with our approach to those from the perturbative theory. As we study twin binary systems, modelling one star is sufficient to describe the entire system. We integrated MoBiDICT to a minimisation method following a Levenberg-Marquardt algorithm to compute on the fly the non-perturbative apsidal motion resulting from twin binaries which unperturbed structure are given by our 1D single star stellar evolution code, CLES. By using this methodology we followed the procedure of Rosu et al. (2020a, 2022a,b, 2023), neglecting the interactions of the binaries during their evolution. To avoid such assumptions, our modelling method would have to be directly coupled to binary stellar evolution codes, this work will be conducted in the future with the code BINSTAR (Siess et al. 2013) and will be presented in a forthcoming article.

We selected four well-known close observed twin binary systems from the literature: PV Cas (Torres et al. 2010; Claret et al. 2021; Marcussen & Albrecht 2022), IM Per (Lacy et al. 2015; Claret et al. 2021), Y Cyg (Gimenez et al. 1987; Harmanec et al. 2014; Claret et al. 2021; Marcussen & Albrecht 2022), and HD152248 (Rosu et al. 2020a,b). The first binary system is composed of MS intermediate mass stars, the second one of sub-giant stars, and the last two ones of massive stars. The properties of these systems are given in the second and fifth line blocks of Table 2. For each of these binary systems, an apsidal motion was determined in the literature and is used in this section as a constraint on the stellar structure. The systems were selected to be twin systems with accurate determination of the orbital and stellar parameters, and an apsidal motion primarily driven by tidal forces.

For this modelling we were inspired by the findings of Rosu et al. (2020a, 2022a,b, 2023), who demonstrated that modelling twin binary systems using 1D stellar evolution codes with apsidal motion as a constraint necessitates significant additional mixing during stellar evolution. In our modelling, we used the same constraints for all observed systems: the masses, radii, effective temperatures, and apsidal motions. The constraints for each system are listed in the second line block of Table 2. The objective of our modelling technique is to generate an evolutionary track such as the end model is fitting all the observed constraints of the system. As we have four constraints, we used four free parameters for our stars, their initial mass M₀, initial hydrogen mass fraction X₀, and age, as well as the overshooting parameter α_over. In our modelling, we used the overshooting to control the additional mixing necessary to reproduce the observed constraints as a strong degeneracy is existing between the overshooting and the mixing (see e.g. Rosu et al. 2020a, 2022a). For all the models, we used the AGSS09 abundances (Asplund et al. 2009), the FreeEOS equation of state (Irwin 2012), the OPAL opacities (Iglesias & Rogers 1996), the T(τ) relation from Model-C of Vernazza et al. (1981) for the atmosphere, the mixing length theory of convection implemented as in Cox & Giuli (1968), and the nuclear reaction rates of Adelberger et al. (2011). For the massive stars, the mass loss rate was computed, accounting for the metallicity of the models, with the prescription of Vink et al. (2001) assuming that the scaling factor for the mass loss rate equals ξ = 1. The equivalent averaged mass loss along the entire evolution is given in the fourth line block of Table 2.

In our modelling procedure, we started by finding models reproducing the observations with the perturbative approach. We then fixed the initial mass and hydrogen mass fraction, only using the constraints on the radii and apsidal motion to model the systems with MoBiDICT. This approach aims to directly evaluate the impact of our non-perturbative modelling on the required additional mixing and age. The objective of our modelling is not to provide an accurate estimation of the stellar properties of each system but rather to compare the stellar properties obtained when modelling the apsidal motion with the perturbative and non-perturbative approaches. The results from our modelling procedure are given both for the non-perturbative and perturbative modelling approaches in the third line block of Table 2.

For the modelling of each system, the parameters presented in the fourth line block of Table 2 were empirically selected to facilitate the convergence of our Levenberg-Marquardt method. For the first system, Z₀ was chosen from a model grid to minimise the χ² and be able to reproduce the observations as we could not achieve a coherent fit with a solar metallicity. For the other systems, a solar metallicity was adopted as no constraints were given from the observations.

Table 2 shows that both the perturbative and non-perturbative approaches can accurately reproduce all observed system constraints. For the four systems, we find that extremely large extra-mixing is required to reproduce the observational properties of stars. Comparing the stellar parameters obtained with MoBiDICT and the perturbative approach, our non-perturbative technique necessitates even more extra mixing for the first three systems. While for the three first systems, the leading contribution to the χ² is the reproduction of the apsidal motion and the radius of the stars, for HD152248, the major difficulty was to reproduce the observed effective temperature. This can partly explain the smaller impact of the non-perturbative treatment on the overshooting parameter found in this system. Nonetheless, our approach does not reduce the extremely high extra-mixing required to reproduce the apsidal motion motion inferred by Rosu et al. (2020a) but rather exacerbate it. This additional mixing could stem from various sources, including rotation or dynamical tides. We do not think that it should be concluded from our study that such large overshooting is present in these stars. Other physical processes neglected in our study such as dynamical tides are known to significantly affect the apsidal motion (Willems & Claret 2002).

5.2. Apsidal motion contribution from the different sources

In this section, we focus on exploring the origin of the apsidal motion discrepancy seen for the observed targets presented in Sect. 5. We reused the structures from the modelling done in Table 2 and looked at the different contributions to the apsidal motion for each binary system. Table 3 provides a decomposition of the apsidal motion for each observed binary system using the models obtained with our non-perturbative method. In the right part of Table 3, we also decomposed the difference between the perturbative and non-perturbative approaches based on the considered spherical orders, revealing the key contributors to the apsidal motion discrepancies.

The analysis of Table 3 reveals that, for all chosen targets, the contributions from our non-perturbative modelling approach surpass the uncertainties in the observational apsidal motion determinations. Generally, the corrections provided by MoBiDICT are comparable in magnitude to the general relativistic corrections. Then, going into more details, the perturbed acceleration discrepancy between the non-perturbative and perturbative methods can be attributed to several main components:

• ℓ = 1: Constituting ∼23% of the modelling discrepancy, this component cannot be recovered by the perturbative model due to its inherent formalism and assumptions;

• ℓ = 2: Representing ∼40% of the modelling discrepancy, this component cannot either be recovered by the perturbative approach. The differences arise from the first-order assumption of the perturbative approach. We tried taking into account the effect from the quadrupolar deformation of the primary star on the quadrupolar component of its companion and the reverse impact of this modification on the quadrupolar deformation of the primary star. However this effect is negligible as each of such additional effect considered has an increase dependency to (R/a)³ of a power three;

• ℓ = 3: Accounting for ∼35% of the modelling discrepancy, this component can be at least partially recovered by the perturbative model by going at higher spherical orders than usually considered. In Sect. 5.3, we expand the computation of the pertubative apsidal motion to the ℓ = 3 assessing the extent to which these discrepancies can be rectified;

• ℓ > 3: The contribution of these components only represents 2% of the discrepancies. Our analysis suggests that computing spherical orders higher than three is unnecessary, in particular when comparing these discrepancies with the contributions from other spherical orders;

• F′ρ′: This term corresponds to the assumption made to neglect the third contribution to the apsidal motion in Eq. (31). With our modelling method we see that this assumption is totally valid even in the most distorted cases.

For the case of IM Per, we examine the individual contributions to this discrepancy. Figure 11 illustrates the evolution of the components of the perturbed acceleration with orbital separation normalised by the stellar radii for the perturbative and non-perturbative models.

For the ℓ = 2 terms, Fig. 11 shows that most of the discrepancies are originating from the ℓ = 2, m = 2 term. This aspect is explained by the fact that in the ℓ = 2, m = 0 term both tidal force and centrifugal force are included while in ℓ = 2, m = 2 only tidal force in present. This indicates that the perturbative approach lacks precision when modelling tidal forces at low orbital separations, corroborating findings from Sect. 4.2 and Fig. 4.

5.3. Higher order perturbative method

An advantage of the perturbative approach lies in the possibility to extend the modelling to higher spherical orders. In the previous section we showed that the ℓ = 3 component of the acceleration perturbation represents ∼35% of the discrepancies seen between our non-perturbative and the perturbative models. In Appendix E, we developed the equations necessary to extend the perturbative approach to the ℓ = 3 term. In this section, we verify the relevance of adding the ℓ = 3 term in the computation of the apsidal motion in the perturbative approach. As shown in Appendix E, the apsidal motion arising from the ℓ = 3 component of the perturbed acceleration is given by

(46)

with

(47)

We took back the model of observed stars extensively presented in Sect. 5.1 and computed the apsidal motion using our method and the perturbative approach including the ℓ = 3 term. The results from this modelling are presented in Table 4.

By comparing Tables 3 and 4 we conclude that by including the ℓ = 3 component of the apsidal motion we corrected ∼20% of the discrepancies seen in Sect. 5.2. However, the ℓ = 3 term is only partially corrected when included in the perturbative approach. Indeed, we found that for all the systems modelled this term stills represent ∼10% of the observed model discrepancies. Due to the reduction of the contribution from the ℓ = 3 term, the quadrupolar term now represents ∼50% of the observed discrepancies while the dipolar term contributes to ∼30% of the differences. The exact repartition of these contributions can vary depending on the orbital separation and the stellar type, in particular, we saw a highly non-linear comportment of the dipolar term. Nonetheless, in the future, the ℓ = 3 contribution to the apsidal motion have to be included when modelling the apsidal motion of close binaries with the perturbative approach.

6. Discussion

In the previous sections, we showed that non-pertubative modelling is necessary for a precise interpretation of observations in close binary systems. Aside from the orbital evolution considerations, the detection of tight systems with apsidal motion is favoured due to the (R/a)⁵ dependency of the tidal component. The combination of high-precision observations and strong tidal effects on apsidal motion motivates their investigation as stellar laboratories and, consequently, the need to model them comprehensively. Our modelling has conclusively shown that, for close binary systems, a non-perturbative treatment of deformations is necessary for accurately characterising stellar global properties and stellar structure.

To assess the impact of our modelling method on observed binary systems with apsidal motion, we constructed a comprehensive catalogue of such systems. In total, our catalogue includes 61 systems, encompassing all systems with constraints on their stellar parameters and observed apsidal motion found in the literature. These systems were obtained by combining two recent catalogues of observed binary systems with apsidal motion (Claret et al. 2021; Marcussen & Albrecht 2022) to systems from various literature sources where stellar parameter determinations were available (Gimenez et al. 1987; Benvenuto et al. 2002; Wolf et al. 2006, 2008, 2010; Torres et al. 2010; Pablo et al. 2015; Baroch et al. 2021, 2022; Rosu et al. 2020b, 2022a,b, 2023). With this catalogue, Fig. 12 illustrates the fraction of observed apsidal motion arising from relativistic corrections against the parameter q^−1/5a(1 − e²)/R₁ of each observed system.

The observed binary systems in our catalogue can be broadly categorised into two equal-sized groups: close binaries with most of the apsidal motion driven by tidal effects and wider binaries where apsidal motion is mostly due to relativistic corrections. The boundary between these two groups is located around q^−1/5a(1 − e²)/R₁ = 8 − 9. Because of their low tidal contribution to the apsidal motion, the group of systems with high relativistic contribution is not favoured when wanting to constrain the stellar structure with apsidal motion. For such systems, as shown in Sect. 4, no significant correction is expected from our non-perturbative method. On the other hand, for binary systems where apsidal motion is dominated by tidal effects (close systems), which is the group preferred for drawing constraints on stellar structure with apsidal motion, our non-perturbative models show discrepancies compared to perturbative models. In Sect. 4.4 we showed that, for a given orbital separation, intermediate-high mass stars are the least impacted by our non-perturbative method, independently from the k₂ value. In Fig. 7 we also showed that for a theoretical massive star the apsidal motion relative difference reaches 5% when q^−1/5a(1 − e²)/R₁ ≲ 4.5 and 2% when q^−1/5a(1 − e²)/R₁ ≲ 6.5. Therefore, we can estimate that in general when q^−1/5a(1 − e²)/R₁ ≲ 4.5 the relative difference of tidally induced apsidal motion reaches at least 5% which corresponds to about 10% of the stars in the catalogue presented in Fig. 12. Similarly, for 42% of the catalogue we expect at least a 2% apsidal motion correction from our methodology, corresponding to almost all the systems with the apsidal motion highly dominated by the tidal component. Among the binaries impacted by our methodology, a majority are high-mass or intermediate-mass binaries. In these catalogues and in general, low-mass binaries are under-represented, and not commonly observed with apsidal motion probably due to observational biases.

The discrepancies observed in the apsidal motion are a consequence of the underlying tidal force discrepancy. However, various aspects and phenomena of binary systems’ evolutions are impacted by an underestimation of the tidal force. One key quantity affected by this underestimation is tidal dissipation, which plays a crucial role in the evolution of binary systems through orbital circularisation, transfer of angular momentum between bodies, and orbital migrations. Current models of the tidal dissipation rely on the perturbative approach (Hut 1981) combined to a free parameter that can be adjusted to control the dissipation. Even if this parameter can be tuned to reproduce the dissipation obtained by our non-pertrubative models with enhanced tidal interactions at a given point of the evolution, it cannot capture the non-linearities that arise during the evolution of binary systems. Consequently, our methodology predicts that tidal dissipation models used by the majority of binary stellar evolution codes underestimate the dissipation at low orbital separation. This leads to imprecise predictions regarding the evolution of global orbital parameters, stellar rotation, and the overall outcomes of these models for systems with close orbits. One specific area impacted by these findings is population synthesis, where the modelled population strongly depends on stellar interactions. Therefore, population models focused on the outcomes of close massive binary systems and the remnants of these systems face significant uncertainties, even before considering mass exchange and processes that rely on imprecise models with free parameters. These findings highlight the need for improved models that better account for the complex interplay of tidal forces in close binary systems.

The discrepancies revealed by our modelling have wide-ranging implications that extend beyond the effects on global stellar and orbital parameters, given that the entire structures of these distorted bodies are affected. The enhanced stellar interaction obtained with our modelling drastically increases the stellar surface deformations, as shown in Fellay & Dupret (2023). Furthermore, the internal structure of stars, which is typically modelled using the perturbative approach, is deeply influenced by our non-perturbative methodology. To precisely capture the interior deformations, mass redistributions, and structural properties of deformed stars, a transition to 3D modelling becomes imperative. Having highly accurate structural models is crucial for the precise study and modelling of the pulsation properties of deformed celestial bodies. Among the oscillations propagating in deformed stars, gravito-inertial waves, called dynamical tides, significantly contribute to the orbital evolution of binary systems. However, the current dynamical tides models rely on structural modelling techniques based on the perturbative approach combined with 1D non-adiabatic spherically symmetric oscillations codes to treat this highly sensitive and non-linear problem. As such, the dynamical tide models, as obtained in Ma & Fuller (2023), for example, are only estimations of the impact of dynamical tides on close binary systems. Without having appropriate 3D non-perturbative deformed structures and the associated 3D non-perturbative non-adiabatic oscillation codes this problem cannot be precisely solved.

Our methodology, similar to the perturbative approach, does have a notable limitation: it assumes solid body rotation within stars. For single stars, this assumption is not fully justified because they often exhibit strong differential rotation, which varies over the course of stellar evolution. However, as we showed in Sect. 5.1, interactions between binary components induce significant extra mixing within deformed stars. This additional mixing reduces the differential rotation in interacting bodies, satisfying our assumption. The subject of angular momentum transport processes in binary systems is a complex and extensive area of research that warrants further investigation. MoBiDICT represents an important initial step in this direction. It allows for the first precise 3D computation of stellar interactions and the resulting forces at each point within the deformed stellar structure.

Finally, the modelling of mass transfer and phases of common envelopes in binary systems is still associated with significant uncertainties and often relies on the use of multiple free parameters. In theory, MoBiDICT has the potential to be adapted for the static modelling of common envelope phases or, at the very least, to provide initial 3D deformed structures that can be used in conjunction with more sophisticated codes to simulate common envelope evolution. Furthermore, our 3D non-perturbative models can be valuable in modelling mass transfer processes in binaries or providing initial 3D deformed structures for further studies. In the near future, it would be useful to investigate the oscilation properties of binaries before they fill their Roche lobes as determining their configurations and properties before these phases play a crucial role in determining their ultimate fate.

7. Conclusion

In this work, we investigated the impact of non-perturbative modelling on the theoretical computation of the apsidal motion and tidal force in eccentric binaries. We developed formalisms for apsidal motion in our non-perturbative approach and the perturbative approach in respectively Sects. 2 and 3. Then, we compared the perturbed acceleration and apsidal motion for theoretical models, in Sect. 4, and observed binary systems, in Sect. 5. Finally, we discussed the implications of our results, in Sect. 6, on the observed population of binaries and the general problem of binary modelling.

All theoretical stellar models were significantly impacted at low orbital separation by our methodology. For instance, we observed a maximum increase in tidal force of about 40% for our low-mass stars, 28% for the RGB, 20% for the solar types, and 12% for massive stars. Consequently, the apsidal motion of these systems is impacted. For an orbital eccentricity of e = 0.1, these discrepancies led to a respectively increase in apsidal motion of 70%, 45%, 30%, and 15% for the low-mass, RGB, solar type, and massive twin binary systems. The discrepancies seen may vary depending on the exact structure of the stars or the architecture and the orbital parameters of the system. These results are consistent with Fellay & Dupret (2023): the models with higher envelope mass experience greater impacts from our methodology. From this study we noted that models with higher k₂ are more impacted by our methodology. We attempted to establish an empirical relationship between tidal force discrepancy, the orbital separation, and k₂. However, due to high non-linearities and the numerous free parameters of the problem we were not able to identify such a relationship. We also explored the dependency of the apsidal motion modelling discrepancy on q and found that this difference is proportional to q except when stars have similar radii. In the latter case, the properties of the secondary are modifying this relationship and the exact properties of the system have to be accounted for.

In our analysis of observed systems, we explored in more details the origins of these discrepancies. Specifically, we modelled four observed twin binary systems, PV Cas, IM Per, Y Cyg, and HD152248, composed of respectively, intermediate mass, sub-giants, and massive stars for the two latests. We started by modelling these systems using the pertubative model to reproduce the observed values of the constraints parameters. When best fitting models were found, we started from this solution to model stars using our non-perturbative method to obtain the theoretical apsidal motion. In agreement with the results of Rosu et al. (2020a, 2022a,b, 2023), our results highlighted the necessity for significant extra mixing in the 1D stellar models to accurately reproduce the observational properties of these systems. Moreover, our methodology induced an additional increase of extra mixing due to the modification of apsidal motion, reinforcing the results of Rosu et al. (2020a, 2022a,b, 2023).

For the systems we examined, the discrepancies in apsidal motion exceeded the observed uncertainties, indicating their non-negligible impact. Our investigation into the origins of these discrepancies revealed that the ℓ = 1 term is responsible for ∼23% of the discrepancies, ∼40% for the ℓ = 2, ∼35% for the ℓ = 3, and less than 2% for the ℓ > 3 components. The exact distribution of these discrepancies may vary based on the model. For the ℓ = 1 and ℓ = 2 terms the discrepancies cannot be recovered by the perturbative approach due to its formalism and approximations. We were able to correct ∼20% of the ℓ = 3 discrepancies by employing a higher spherical order perturbative approach. In opposition to the results of Rosu et al. (2020b), the ℓ = 3 contribution to the apsidal motion cannot be neglected in close binaires.

Finally, we constructed a catalogue of observed binary systems with apsidal motion from the literature to assess the fraction of the observed binaries with apsidal motion effectively impacted by our modelling method. Our analysis indicated that differences in tidal force, and consequently apsidal motion, are non-negligible for systems where q^−1/5a(1 − e²)/R₁ ≲ 6.5. As a result, a significant fraction of the observed binaries (42%) with apsidal motion are indeed affected by our modelling approach.

When modelling the dynamics and deformations of binaries the perturbative approach is the most sophisticated method employed. In this article, we demonstrated that this approach loses precision for binaries with q^−1/5a(1 − e²)/R₁ ≲ 6.5. The most immediate observable quantity affected is the apsidal motion. While, in this work, we only considered the contribution from the equilibrium tides and the general relativistic contributions, additional phenomenon could impact the determination of the apsidal motion. Furthermore, including the dynamical tides can highly modify the theoretical determination of the apsidal motion. However, a proper non-perturbative non-addiabatic 3D astreroseimic modelling is required to study the propagation of gravito-inertial waves in the stellar interior and their impact on the dynamics of binary systems. In the future, our aim is to develop such a tool. In addition, in the future, the inclination of binary systems have to be considered when modelling the tidal contribution to the apsidal motion to account for all the parameters of the orbits. An extension of MoBiDICT for non-aligned systems can be developed; however, all the system’s symmetries exploited in this work will disappear, consequently our methodology will require significantly more resources to run.

Importantly, the underestimations seen extends beyond the apsidal motion to encompass tidal interactions between the binary components, thereby impacting tidal dissipation processes. Notably, a vast majority of binary stellar evolution codes are also subject to these imprecisions at low orbital period, which directly influences the predictive accuracy of binary system outcomes. Given the significant tidal force discrepancies of our models at low orbital separations, it is reasonable to anticipate an underestimation of tidal dissipations in such conditions. Therefore, with our methodology, we expect to amplify the interactions and exchanges of angular momentum in close binary systems, ultimately leading to shorter dynamical evolution timescales. To solidify these observations, the next step is to integrate our method into existing stellar binary evolution codes. This implementation will offer a more comprehensive and accurate understanding of the dynamics and evolution of close binary systems.

---

Acknowledgments

The authors are thanking the anonymous referee for their comments. L.F was supported by the Fonds de la Recherche Scientifique F.R.S.-FNRS as a Research Fellow.

References

Appendix A: The Clairaut-Radau equation

A partial derivation of the Clairaut-Radau equation can already be found in Sterne (1939), Kopal (1959), Fitzpatrick (2012). In this section we provide a detailed proper derivation of this equation and focus on the implication of the assumptions of the perturbative approach. Let us express the radial coordinate of a given equipotential in the direction r₁(θ, ϕ) with its spherical expansion:

(A.1)

where r̅ denotes the averaged radial coordinate of the considered equipotential and the spectral terms of the equipotential radii that can be obtained with

(A.2)

where the spherical harmonics were appropriately normalised. Throughout this entire section, Fig. A.1 illustrates the different variables and orbital quantities used to describe the deformations of a binary system.

In the perturbative approach, r̅ can be obtained by expanding Eq. A.1 to the first order ():

(A.3)

The density of an equipotential passing by (r₁, θ, ϕ), denoted ρ(r̅(r₁,θ,ϕ)), can be approximated to the first order using a Taylor expansion of r̅ around r₁:

(A.4)

where ρ(r₁) is the density of the equipotential with r̅ = r₁. The densities of a sphere of radius r₁ projected on the spherical harmonics basis is obtained through

(A.5)

that can be rewritten with

(A.6)

Using the orthogonality properties of the spherical harmonics, at the first order the projected densities are

(A.7)

The spectral development of the gravitational potential on a sphere of radius r₁ is given by Eq. 24. The gravitational potential on a sphere of radius r₁ will be noted in the future Ψ_r(r₁, θ, ϕ) while the gravitational potential on an equipotential is noted Ψr̅(r̅,θ,ϕ) and can be expressed as

(A.8)

that can be developed to the first order using the Taylor approximation around r₁ = r̅:

(A.9)

and therefore

(A.10)

The spectral components of the gravitational potential (Eq. 24) on an equipotential can be rewritten using Eq. A.7 before integrating by parts:

(A.11)

(A.12)

(A.13)

where at the second step the two crossed terms are cancelling each others out.

As the total potential is constant on an equipotential, one can write that:

(A.14)

implying that for all the ℓ and m different from 0,

(A.15)

Eq. A.13 can be developed and multiplied by r̅^ℓ+1/4πG:

(A.16)

At first order and , By differentiating Eq. A.16, we obtain

(A.17)

This expression can then be divided by r̅^2ℓ and again differentiated by r̅ to obtain

(A.18)

where we used the definition of the mean density under an isobar:

(A.19)

To the first order meaning that therefore the last term of the previous equation is cancelling out during the differentiation. Finally by defining

(A.20)

and developing the second derivative of :

(A.21)

we obtain the well known Clairaut-Radau equation expressed as

(A.22)

Appendix B: Perturbation of the surface total potential

In this section our aim is to look at how the perturbative model is modifying the surface gravitational potential of each star composing the system. Let us go back to Eq. A.16 and divide it by r̅^2ℓ+1 and next differentiate it to make fall:

(B.1)

The previous equation is then multiplied by before being injected in Eq. A.13 to eliminate and obtain the spectral projections of the gravitational potentials expressed as

(B.2)

The integral term is negligible at the surface. For an arbitrary chosen primary with a radius R₁ the surface gravitational potential is now expressed as

(B.3)

Eq. A.17 is then evaluated at the stellar surface and to the first order to obtain

(B.4)

that we inject in Eq. B.3 to express the perturbation of the spectral surface gravitational potential

(B.5)

With the definition of the Love numbers given in Eq. 28, we have finally:

(B.6)

Appendix C: First terms of the centrifugal, gravitational, and tidal potentials

C.1. Centrifugal potential in unsynchronised systems

In a non-synchronised system, the centrifugal potential expressed with the Cartesian coordinates centred on the primary is given by

(C.1)

where (x, y, z) are the Cartesian coordinates centred on the star studied, x pointing towards the centre of mass of the system, ñ is the orbital rotation rate of the binary system with a separation r and Ω_⋆ is the rotation rate of this star. The quantity x² + y² can be written when passing in spherical coordinates as with z = r₁ cos θ that can be expressed with the spherical harmonics

(C.2)

The centrifugal potential can therefore be expressed as

and decomposed as

(C.3)

(C.4)

(C.5)

C.2. Gravitational and tidal potential

Treating the secondary as a point source, the spectral component of the gravitational potential exerted by the secondary star at the surface of the primary can be expressed as

(C.6)

where r is the instantaneous orbital separation of the stars and P_ℓ(sin θ₁ cos ϕ₁) is the classical Legendre polynomial that can be decomposed using spherical harmonics as

(C.7)

where are constant coefficients linked to the normalisation of the spherical harmonics. In the following we define and use

(C.8)

Using this notation:

(C.9)

yield . Similarly,

(C.10)

yields . Therefore,

(C.11)

(C.12)

(C.13)

Summing Eqs. (C.3), (C.11) and using the definition of ñ (Eq. 5), we have:

(C.14)

Substituting this result in the right hand side of Eq. (B.6), we see that the perturbative theory predicts that the dipolar component of the gravitational potential generated by a star is equal to zero when its companion is treated as a point source. This is an important result explaining why the dipolar component was systematically neglected in previous studies.

The sums of the quadrupolar terms of the centrifugal and secondary tidal forces at the surface of the primary are finally given by

(C.15)

and

(C.16)

Appendix D: Perturbation of the force

In this section, we develop the expression of the tidal force exerted by a deformed body on its companion, specifically F_2′10 as introduced in Sect. 3.2. We limit ourselves to the dominant term, which arises from the quadrupolar deformation of the secondary. The system of coordinates used in this Section is described in Fig. A.1. Let us consider the secondary, which is deformed due to the presence of an unperturbed primary. In the section, r_i denotes the radial coordinate of the spherical coordinate frame centred on each star. With this formalism the potential perturbation at the surface of the secondary is given by Eq. B.5 applied to the ℓ = 2 terms:

(D.1)

that can be rewritten, using Eq. C.6 as

(D.2)

where P₂ is the Legendre polynomial of order two and λ_i ≡ sin θ_i cos ϕ_i and μ_i = cos θ_i. To study the impact of this potential perturbation on the primary, we can express the ℓ = 2 component of the potential outside of the secondary as

(D.3)

Following Kopal (1959), different coordinates related quantities from the secondary star can be approximated with series developments in the coordinates of the primary star:

(D.4)

(D.5)

(D.6)

and

(D.7)

All those approximations can be inserted back in Eq. D.3 to obtain the perturbation of the secondary gravitational potential in the coordinates of the primary:

(D.8)

(D.9)

that can be rewritten using the third Kepler law including the periods of the system as follows

(D.10)

We can note that all the constants appearing in the expressions have been included in the Ψ_2, 0 as they are not depending on r₁/r.

Let us now compute the force exerted by the perturbation of the secondary on the unperturbed primary, F_2′10, that is given by

(D.11)

when treating the primary as point like (r₁ = 0). In this case, only the gradient of the ℓ = 1 component is non-null at the centre of the star (r₁ = 0) as

(D.12)

Therefore, by only taking the ℓ = 1 component, the force can be expressed as

(D.13)

Appendix E: Higher order perturbations

In this section, we apply the developments done in Appendix D and Sect. 3.3 to obtain the apsidal motion resulting from the ℓ = 3 perturbation of the potential. By combining Eq. (B.6) and Eq. C.6, the surface potential perturbation ℓ of the secondary originating from the presence of the primary is given by

(E.1)

Beyond the surface of the secondary, at a distance r₂ from its centre, the potential perturbation from the ℓ = 3 perturbed potential of the secondary is expressed as

(E.2)

(E.3)

Based on Kopal (1959), λ₂ can be expended as

(E.4)

meaning that

(E.5)

and

(E.6)

which can be inserted back in Eq. (E.3) to obtain the perturbation originating from the ℓ = 3 perturbation of the secondary in the coordinate frame of the primary:

(E.7)

By using the same methodology as in the Appendix D, the force perturbation originating from the ℓ = 3 deformation of the secondary on the primary is expressed as

(E.8)

Following the Sect. 3.2, the acceleration perturbation undergone by the system due to the tidal perturbation is expressed as

(E.9)

(E.10)

when defining a as the semi-major axis of the system. Following the same approach as in Sect. 3.3, the time orbital average of the perturbation is formally expressed as

(E.11)

which can be developed into polynomial functions of the eccentricity using Hansen’s functions:

(E.12)

with

(E.13)

The apsidal motion caused by the ℓ = 3 acceleration perturbation can finally be obtained by inserting Eq. (E.12) in Eq. 19 and simplifying:

(E.14)

with,

(E.15)

All Tables

All Figures