Issue 
A&A
Volume 683, March 2024



Article Number  A210  
Number of page(s)  20  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202348134  
Published online  26 March 2024 
Underestimation of the tidal force and apsidal motion in close binary systems by the perturbative approach: Comparisons with nonperturbative models
^{1}
STAR Institute, University of Liège, 19C Allée du SixAoût, 4000 Liège, Belgium
email: lfellay@uliege.be
^{2}
Department of Physics, KTH Royal Institute of Technology, The Oskar Klein Centre, AlbaNova, 106 91 Stockholm, Sweden
Received:
2
October
2023
Accepted:
23
December
2023
Context. Stellar deformations play a significant role in the dynamical evolution of stars in binary systems, impacting the tidal dissipation and the outcomes of mass transfer processes. The prevalent method for modelling the deformations and tidal interactions of celestial bodies solely relies on the perturbative approach, which assumes that stellar deformations are minor perturbations to the spherical symmetry. An observable consequence of stellar deformations is the apsidal motion in eccentric systems, which has be observationally determined across numerous binary systems.
Aims. Our objective is to assert the reliability of the perturbative approach when applied to close and strongly deformed binary systems.
Methods. We have developed a nonperturbative 3D modelling method designed to account for high stellar deformations. We focus on comparing the properties of perturbatively deformed stellar models with our 3D models, particularly in terms of apsidal motion.
Results. Our research highlights that the perturbative model becomes imprecise and underestimates the tidal force and rate of apsidal motion at a short orbital separation. This discrepancy primarily results from the firstorder treatment in the perturbative approach, and cannot be rectified using straightforward mathematical corrections due to the strong nonlinearity and numerous parameters of the problem. We have determined that our methodology affects the modelling of approximately 42% of observed binary systems with measured apsidal motion, introducing a discrepancy greater than 2% when the normalised orbital separation verifies q^{−1/5}a(1 − e^{2})/R_{1} ≲ 6.5 (q is the mass ratio of the system, a is its semimajor axis, e is its orbital eccentricity and R_{1} is the radius of the primary star).
Conclusions. The perturbative approach underestimates tidal interactions between bodies up to ∼40% for close lowmass binaries. All the subsequent modelling is impacted by our findings, in particular, the tidal dissipation is significantly underestimated. As a result, all binary stellar models are imprecise when applied to systems with a low orbital separation, and the outcomes of these models are also affected by these inaccuracies.
Key words: celestial mechanics / binaries: close / binaries: general / stars: evolution / stars: interiors
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
Tidal forces in binary systems result from the gravitational interaction between nonpointlike bodies and can be categorised into two components: the equilibrium tides, corresponding to largescale 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^{1} (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 wellknown stellar structure constant, the k_{2} coefficient that is obtained from solving the ClairautRadau 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 nonperturbative 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 nonperturbative 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 nonperturbative 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 nonperturbative 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. Nonperturbative modelling
2.1. Setup 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 nonperturbative 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:
and the Poisson equation,
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 nonperturbative 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} = Ψ_{1} + Ψ_{2} + Ψ_{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 nonperturbative modelling method. In this case, we applied our tool to a binary system composed of a primary 1 M_{⊙} mainsequence (MS) star and a 0.2 M_{⊙} MS companion separated by r = 2R_{1}, R_{1} being the radius of the primary given by the 1D spherically symmetric model.
Fig. 1. View of a 1.0 M_{⊙} MS primary star in a binary system, cut in the plane including all the rotations axes of the system. The companion to this star is a 0.2 M_{⊙} MS secondary star with an orbital separation of r = 2R_{1}. The black curve corresponds to the Roche lobe of the primary star. The black lines respectively correspond to the reference direction for θ and to the reference direction in witch our model is given by the 1D spherical symmetric model (θ = μ_{crit}). The blue lines are equipotentials inside the star with the colourbar representing the total potential. On each equipotential the pressure, density, and temperature, are constant. L_{1} is the Lagrangian point and R_{1} is the radius of the initial 1D model in the reference direction. 
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, pseudosynchronisation 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 nonsynchronised binary systems is necessary to model eccentric binaries and binary evolution.
2.2. Nonperturbative modelling of nonsynchronised 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.
Fig. 2. Illustration of a binary system configuration with nonsynchronised rotations. The system is composed of 1 M_{⊙} MS primary star and a 0.2 M_{⊙} MS companion separated by r = 2R_{1}. The two grey planes are the planes of symmetry, each rotation axis is marked with a black line. The orbital rotation is originating from the system centre of mass, located inside the primary. The black curve is the secondary orbital path, a black axis is joining the centre of each star with and emphasis on the Lagrangian point L1. Finally, the colour code corresponds to the surface effective gravity, the local flux emitted by a star being a function of this quantity. 
With our nonperturbative 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
and becomes in the nonsynchronised cylindrical symetric rotation and circular orbit case:
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 majoraxis 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^{2} has to be replaced, in Eq. (4), by ñ^{2} defined as
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
where the bold symbols denote vectorial quantities, is the unit vector linked to the xaxis 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
where the bold symbols denote vectorial quantities, F_{21} is the force exerted by the secondary star on the primary (noting F_{2 → 1} ≡ F_{21}), F_{12} 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
The last step of Eq. (7) was deduced from the third Newton Law (F_{12} = −F_{21}). The total force per unit of mass exerted by the secondary on its primary can be decomposed as
In Eq. (9), F_{R} is the radial perturbed acceleration resulting from considering stars as nonspherical extended interacting bodies. This perturbation is responsible for all the equilibrium tidalinduced dynamical effects of binary evolution, including orbital migration (i.e. variation of systems semimajor 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
where F_{2} is the gravitational force from the secondary inside the primary. By decomposing the densities and the forces in unperturbed spherical symmetric terms, respectively ρ_{10} and F_{20}, and nonspherical perturbations, respectively and , the force exerted by the secondary on the primary can be rewritten as
The spherically symmetric component of the force can be simply identified as being
Therefore, the perturbed acceleration is expressed as
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
and
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 nonsymmetric terms is straightforward. The ℓ = 0, m = 0 term is the spherically symmetric component of each quantity while all other spherical orders are nonsymmetric terms. In our formalism the perturbed acceleration can thus be obtained with
Finally, each force term projected on the spherical harmonic basis can be decomposed and reexpressed in a spherical coordinate frame (to be integrated) with
The spherical harmonics Y_{ℓ, m} in Eq. (17) are functions of μ_{2} and ϕ_{2}. The expressions of the , and terms can be found in Fellay & Dupret (2023). More than having a nonperturbative 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. Nonperturbative apsidal motion computation
In this section, we developed the formalism used to model the apsidal motion with our nonperturbative 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.
Fig. 3. Orbital trajectory of a star in a binary system with an eccentricity of 0.2. The orbital plane is in red, the grey plane is a the reference plane (i.e. the plane of the sky), i is the inclination of the system compared to the reference plane, Ω is the longitude of the ascending node, Ω_{0} is the ascending node, ω is the argument of the periastron, and φ is the true anomaly. The reds dots on the trajectory represent the different points at which the snapshots are taken to compute the apsidal motion. 
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
where a is the system’s semimajor 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
where φ is the true anomaly and ⟨F_{R} cos φ⟩ denotes the time average of the perturbed acceleration expressed as
To obtain the apsidal motion of a system through a nonperturbative approach, the sole quantity to be computed by our nonperturbative 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 timeaveraged perturbed acceleration by reformulating its expression with the second Kepler’s law:
where r is the separation between the stars expressed, following the first Kepler’s Law, as
By combining Eqs. (19)–(22), the apsidal motion obtained with our nonperturbative modelling is given by
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^{−3} − 10^{−4} on the apsidal motion in the vast majority of close binaries.
3. Perturbative model
3.1. The ClairautRadau 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,
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 ClairautRadau firstorder differential equation Kopal (1959), and derived in Appendix A) expressed as
using the boundary conditions
In the ClairautRadau 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 pointmass 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
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 nonspheroidal deformation response of a body subjected to a perturbative potential. Mathematically, these coefficients are defined as:
Equation (27) provides a direct comparison between the models, where the quantity is equivalent to of our nonperturbative 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 pointmass body is expressed, in the nonsynchronised case, as
and
In both synchronised and nonsynchronised systems, the dipolar component of the total potential is null when the secondary is treated as a pointmass, 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:
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
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
with q = M_{2}/M_{1}.
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
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
The integrals in Eq. (35) can be reexpressed with the zeroorder Hansen coefficients defined for every positive m as
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
and,
The average perturbed acceleration is finally given by
and can be inserted in Eq. (19) to obtain the apsidal motion of a binary system with the perturbative approach:
with,
and,
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 nonperturbative 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).
Summary of the stellar properties of the 1D input models used in this work.
Table 1 in Fellay & Dupret (2023) has only been extended to include the k_{2} value of each model. As shown by Eq. (27), to a higher k_{2} corresponds a greater surface response from a body when subjected to a potential perturbation to its surface. Moreover, k_{2} can be used as an indicator of a body’s deformability. In addition, our results from Fellay & Dupret (2023), showed that the lowmass 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 nonperturbative method for a given a/R_{1} ratio. In the following sections, we compared the surface deformations and apsidal motion obtained with our nonperturbative 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 nonperturbative 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:
where is the spectral surface potential obtained with our nonperturbative 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.
Fig. 4. Surface potential discrepancy between perturbative and nonperturbative modelling approaches as a function of the orbital separation normalised by the stellar radii for the different binary components presented in Table 1. The upper panel corresponds to the ℓ = 2, m = 0 component of the potential while the lower panel is the ℓ = 2, m = 2 term. 
The surface potential discrepancy, varies significantly depending on the stellar type and the orbital separation. For both ℓ = 2 components, stars with the higher k_{2} exhibit more pronounced effects from our nonperturbative 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 nonperturbative 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 nonperturbative 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 nonperturbative perturbed accelerations, we define their difference ΔF_{R} as follows:
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.
Fig. 5. Perturbed acceleration discrepancy as a function of the separation normalised by the stellar radii for the twin binary systems composed of the stars presented in Table 1. The colour code denotes the different twin binary systems. 
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_{2}, namely the low mass and RGB stars, are the most impacted by our modelling. Solarlike 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:
to compare the apsidal motion originating from the nonperturbative and perturbative approaches. In the perturbative approach, the apsidal motion is given by Eq. (40). With our nonperturbative 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.
Fig. 6. Evolution of the apsidal motion relative difference as a function of the orbital separation normalised by the stellar radii. The eccentricity of the orbits were set to e = 0.1 and the colour code denotes the different twin binary systems. 
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 semimajor 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^{2}. 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^{2})/R_{1} ≲ 4.5 and 2% when a(1 − e^{2})/R_{1} ≲ 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^{2})/R_{1} ≲ 6.5 and a(1 − e^{2})/R_{1} ≲ 4.5.
Fig. 7. Apsidal motion relative different for a twin binary system composed of two 20 M_{⊙} stars presented in Table 1. Each colour correspond to the variation of discrepancy as a function of the orbital separation scaled by the eccentricity and stellar radii for different values of the orbital eccentricity. The black lines are lines of constant model relative difference, corresponding to 5% and 2% discrepancy. 
In addition to this work, we verified the precision of our method at higher orbital separation and found that for a(1 − e^{2})/R_{1} = 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_{2} 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_{2} 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_{2} to explore whether discrepancies are solely dependent on k_{2} 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_{⊙}.
Fig. 8. Mapping of the acceleration perturbation discrepancy for our grid of lowmass stars and intermediatehigh mass stars. The colour code denotes the difference of perturbed acceleration between the perturbative and our nonpeturbative approach as a function of the models k_{2} and orbital separation. 
Figure 8 reveals a distinct behaviour between the two groups of stars. While the lower k_{2} limit of the low mass stars and the upper limit of the intermediatehigh mass stars have similar k_{2}, their acceleration perturbation discrepancy is not similar. This result indicates that discrepancies between the models are not only a function of k_{2} 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_{2} 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_{2} and orbital separation to correct the modelling discrepancies in the perturbative approach, the nonlinear 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 intermediatehigh mass stars and lowmass 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. Nonperturbative 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_{2} 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_{1}/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/Ω_{1})^{2}, (n/Ω_{2})^{2}, (R_{1}/a)^{5}, (R_{2}/a)^{5}, q, e, k_{21}, k_{22}), see Eq. (40). With our nonpertubative modelling the dependencies on k_{21}, k_{22} 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 subproduct 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_{2}/a)^{5}, q, k_{22} or ρ_{2}(r)). Consequently the effect of a change of q will be seen indirectly, through, in particular (R_{2}/a)^{5}. 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_{2}/R_{1})^{5}, 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} = 0.72, Z_{0} = 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 xaxis of Fig. 9 to be q^{−1/5}R_{1}/a to maintain the scale between (R_{1}/a)^{5} and q.
Fig. 9. Apsidal motion discrepancy as a function of the orbital separation scaled with the q parameter. Each system is composed of a primary 1.5 M_{⊙} RGB star presented in Table 1, and a secondary star with a lower mass. The stars composing the systems have the same age (2 Byr) and initial chemical composition (X_{0} = 0.72, Z_{0} = 0.010). Each colour corresponds to a secondary with a different mass. 
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_{1}/a)^{5} 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} = 0.72, Z_{0} = 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/5}R_{1}/a for these binary systems.
Fig. 10. Apsidal motion discrepancy as a function of the orbital separation scaled with the q parameter. Each system is composed of a primary 1.5 M_{⊙} MS star, and a secondary star with a lower or equal mass. The stars composing the systems have the same age (1.5 Byr) and initial chemical composition (X_{0} = 0.72, Z_{0} = 0.010). Each colour corresponds to a secondary with a different mass. 
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/5}a(1 − e^{2})/R_{1} ≲ 4.5 to have at least 5% modelling discrepancy and q^{−1/5}a(1 − e^{2})/R_{1} ≲ 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 nonperturbative 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 LevenbergMarquardt algorithm to compute on the fly the nonperturbative 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 wellknown 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 subgiant 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.
Observed and model properties of the systems studied.
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_{0}, initial hydrogen mass fraction X_{0}, 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 ModelC 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 nonperturbative 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 nonperturbative approaches. The results from our modelling procedure are given both for the nonperturbative 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 LevenbergMarquardt method. For the first system, Z_{0} was chosen from a model grid to minimise the χ^{2} 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 nonperturbative approaches can accurately reproduce all observed system constraints. For the four systems, we find that extremely large extramixing is required to reproduce the observational properties of stars. Comparing the stellar parameters obtained with MoBiDICT and the perturbative approach, our nonperturbative technique necessitates even more extra mixing for the first three systems. While for the three first systems, the leading contribution to the χ^{2} 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 nonperturbative treatment on the overshooting parameter found in this system. Nonetheless, our approach does not reduce the extremely high extramixing 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 nonperturbative method. In the right part of Table 3, we also decomposed the difference between the perturbative and nonperturbative approaches based on the considered spherical orders, revealing the key contributors to the apsidal motion discrepancies.
Decomposition of the apsidal motion of the modelled observed systems.
The analysis of Table 3 reveals that, for all chosen targets, the contributions from our nonperturbative 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 nonperturbative 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 firstorder 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)^{3} 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 nonperturbative models.
Fig. 11. Decomposition of the perturbed acceleration for the IMPer system with the perturbative and nonperturbative models as a function of the system orbital separation normalised by the stellar radii. Each colour represents one spherical order: the lighter curves correspond to the perturbative model while the darker curves correspond to the with nonperturbative approach. The red curves are the total of all the spherical orders. 
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.
Decomposition of the apsidal motion of the modelled observed systems including the ℓ = 3 term of the perturbative approach.
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 nonperturbative 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
with
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 nonlinear 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 nonpertubative 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)^{5} dependency of the tidal component. The combination of highprecision 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 nonperturbative 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/5}a(1 − e^{2})/R_{1} of each observed system.
Fig. 12. Repartition of observed binary systems from our catalogue affected by our methodology. Each point correspond to an observed system with a given orbital separation and a relativistic contribution to the apsidal motion. The blue dots correspond to massive binaries, the yellow points are intermediate mass binary stars, and the red dots are lowmass stars. For this diagram, the primary was chosen as the most massive star in the system. 
The observed binary systems in our catalogue can be broadly categorised into two equalsized 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/5}a(1 − e^{2})/R_{1} = 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 nonperturbative 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 nonperturbative models show discrepancies compared to perturbative models. In Sect. 4.4 we showed that, for a given orbital separation, intermediatehigh mass stars are the least impacted by our nonperturbative method, independently from the k_{2} value. In Fig. 7 we also showed that for a theoretical massive star the apsidal motion relative difference reaches 5% when q^{−1/5}a(1 − e^{2})/R_{1} ≲ 4.5 and 2% when q^{−1/5}a(1 − e^{2})/R_{1} ≲ 6.5. Therefore, we can estimate that in general when q^{−1/5}a(1 − e^{2})/R_{1} ≲ 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 highmass or intermediatemass binaries. In these catalogues and in general, lowmass binaries are underrepresented, 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 nonpertrubative models with enhanced tidal interactions at a given point of the evolution, it cannot capture the nonlinearities 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 wideranging 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 nonperturbative 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, gravitoinertial 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 nonadiabatic spherically symmetric oscillations codes to treat this highly sensitive and nonlinear 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 nonperturbative deformed structures and the associated 3D nonperturbative nonadiabatic 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 nonperturbative 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 nonperturbative modelling on the theoretical computation of the apsidal motion and tidal force in eccentric binaries. We developed formalisms for apsidal motion in our nonperturbative 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 lowmass 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 lowmass, 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_{2} are more impacted by our methodology. We attempted to establish an empirical relationship between tidal force discrepancy, the orbital separation, and k_{2}. However, due to high nonlinearities 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, subgiants, 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 nonperturbative 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 nonnegligible 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 nonnegligible for systems where q^{−1/5}a(1 − e^{2})/R_{1} ≲ 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/5}a(1 − e^{2})/R_{1} ≲ 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 nonperturbative nonaddiabatic 3D astreroseimic modelling is required to study the propagation of gravitoinertial 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 nonaligned 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.
In the case of a hierarchical triple system, a contribution to the apsidal motion arises from the third body orbiting the inner binary (Naoz et al. 2013). In the present study, we do not account for such an effect as we focus on binary systems.
In Fig. 7, we have chosen the empirical parameter a(1 − e^{2})/R_{1} for the xaxis such as the curves are the most superposed at high orbital separation.
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
 Adelberger, E. G., García, A., Robertson, R. G. H., et al. 2011, Rev. Mod. Phys., 83, 195 [Google Scholar]
 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Baroch, D., Giménez, A., Ribas, I., et al. 2021, A&A, 649, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baroch, D., Giménez, A., Morales, J. C., et al. 2022, A&A, 665, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Benvenuto, O. G., Serenelli, A. M., Althaus, L. G., Barbá, R. H., & Morrell, N. I. 2002, MNRAS, 330, 435 [NASA ADS] [CrossRef] [Google Scholar]
 Claret, A. 2023, A&A, 674, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Claret, A., Giménez, A., Baroch, D., et al. 2021, A&A, 654, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cox, J. P., & Giuli, R. T. 1968, Principles of Stellar Structure (New York: Gordon and Breach) [Google Scholar]
 Fellay, L., & Dupret, M. A. 2023, A&A, 676, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fitzpatrick, R. 2012, An Introduction to Celestial Mechanics (UK: Cambridge University Press) [CrossRef] [Google Scholar]
 Gimenez, A., & Margrave, T. E. 1985, AJ, 90, 358 [NASA ADS] [CrossRef] [Google Scholar]
 Gimenez, A., Kim, C.H., & Nha, I.S. 1987, MNRAS, 224, 543 [NASA ADS] [CrossRef] [Google Scholar]
 Harmanec, P., Holmgren, D. E., Wolf, M., et al. 2014, A&A, 563, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hut, P. 1981, A&A, 99, 126 [NASA ADS] [Google Scholar]
 Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
 Irwin, A. W. 2012, Astrophysics Source Code Library [record ascl:1211.002] [Google Scholar]
 Jeans, J. H. 1929, The Universe Around Us (New York: The Macmillan company) [Google Scholar]
 Kopal, Z. 1959, Close Binary Systems (London: Chapman& Hall) [Google Scholar]
 Kopal, Z. 1978, Dynamics of Close Binary Systems (Dordrecht: Reidel) [Google Scholar]
 Lacy, C. H. S., Torres, G., Fekel, F. C., Muterspaugh, M. W., & Southworth, J. 2015, AJ, 149, 34 [CrossRef] [Google Scholar]
 Ma, L., & Fuller, J. 2023, ApJ, 952, 53 [CrossRef] [Google Scholar]
 Marcussen, M. L., & Albrecht, S. H. 2022, ApJ, 933, 227 [NASA ADS] [CrossRef] [Google Scholar]
 Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., & Teyssandier, J. 2013, MNRAS, 431, 2155 [NASA ADS] [CrossRef] [Google Scholar]
 Pablo, H., Richardson, N. D., Moffat, A. F. J., et al. 2015, ApJ, 809, 134 [NASA ADS] [CrossRef] [Google Scholar]
 Packet, W. 1981, A&A, 102, 17 [NASA ADS] [Google Scholar]
 Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
 Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [Google Scholar]
 Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
 Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10 [Google Scholar]
 Rosu, S., Noels, A., Dupret, M. A., et al. 2020a, A&A, 642, A221 [EDP Sciences] [Google Scholar]
 Rosu, S., Rauw, G., Conroy, K. E., et al. 2020b, A&A, 635, A145 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rosu, S., Rauw, G., Farnir, M., Dupret, M. A., & Noels, A. 2022a, A&A, 660, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rosu, S., Rauw, G., Nazé, Y., Gosset, E., & Sterken, C. 2022b, A&A, 664, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rosu, S., Quintero, E. A., Rauw, G., & Eenens, P. 2023, MNRAS, 521, 2988 [NASA ADS] [CrossRef] [Google Scholar]
 Scuflaire, R., Théado, S., Montalbán, J., et al. 2008, ApSS, 316, 83 [NASA ADS] [Google Scholar]
 Siess, L., Izzard, R. G., Davis, P. J., & Deschamps, R. 2013, A&A, 550, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sterne, T. E. 1939, MNRAS, 99, 451 [NASA ADS] [Google Scholar]
 Torres, G., Andersen, J., & Giménez, A. 2010, A&A Rev., 18, 67 [Google Scholar]
 Vernazza, J. E., Avrett, E. H., & Loeser, R. 1981, ApJS, 45, 635 [Google Scholar]
 Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Willems, B., & Claret, A. 2002, A&A, 382, 1009 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wolf, M., Kučáková, H., Kolasa, M., et al. 2006, A&A, 456, 1077 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wolf, M., Zejda, M., & de Villiers, S. N. 2008, MNRAS, 388, 1836 [Google Scholar]
 Wolf, M., Claret, A., Kotková, L., et al. 2010, A&A, 509, A18 [CrossRef] [EDP Sciences] [Google Scholar]
 Zahn, J. P. 1975, A&A, 41, 329 [NASA ADS] [Google Scholar]
 Zahn, J. P. 1977, A&A, 57, 383 [Google Scholar]
Appendix A: The ClairautRadau equation
A partial derivation of the ClairautRadau 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_{1}(θ, ϕ) with its spherical expansion:
where r̅ denotes the averaged radial coordinate of the considered equipotential and the spectral terms of the equipotential radii that can be obtained with
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.
Fig. A.1. Scheme of a binary system separated by a distance r. The different quantities used to model the deformations in the perturbative framework are detailed and illustrated here. The blue curve is the surface of the deformed star while the red curves are the spheres of equivalent averaged radii. As the surface is represented here, r̅_{i}(μ_{i},ϕ_{i}) = R_{i}, R_{i} being the radial coordinate of the stellar surface in the spherically symmetric case. 
In the perturbative approach, r̅ can be obtained by expanding Eq. A.1 to the first order ():
The density of an equipotential passing by (r_{1}, θ, ϕ), denoted ρ(r̅(r_{1},θ,ϕ)), can be approximated to the first order using a Taylor expansion of r̅ around r_{1}:
where ρ(r_{1}) is the density of the equipotential with r̅ = r_{1}. The densities of a sphere of radius r_{1} projected on the spherical harmonics basis is obtained through
that can be rewritten with
Using the orthogonality properties of the spherical harmonics, at the first order the projected densities are
The spectral development of the gravitational potential on a sphere of radius r_{1} is given by Eq. 24. The gravitational potential on a sphere of radius r_{1} will be noted in the future Ψ_{r}(r_{1}, θ, ϕ) while the gravitational potential on an equipotential is noted Ψr̅(r̅,θ,ϕ) and can be expressed as
that can be developed to the first order using the Taylor approximation around r_{1} = r̅:
and therefore
The spectral components of the gravitational potential (Eq. 24) on an equipotential can be rewritten using Eq. A.7 before integrating by parts:
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:
implying that for all the ℓ and m different from 0,
Eq. A.13 can be developed and multiplied by r̅^{ℓ+1}/4πG:
At first order and , By differentiating Eq. A.16, we obtain
This expression can then be divided by r̅^{2ℓ} and again differentiated by r̅ to obtain
where we used the definition of the mean density under an isobar:
To the first order meaning that therefore the last term of the previous equation is cancelling out during the differentiation. Finally by defining
and developing the second derivative of :
we obtain the well known ClairautRadau equation expressed as
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:
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
The integral term is negligible at the surface. For an arbitrary chosen primary with a radius R_{1} the surface gravitational potential is now expressed as
Eq. A.17 is then evaluated at the stellar surface and to the first order to obtain
that we inject in Eq. B.3 to express the perturbation of the spectral surface gravitational potential
With the definition of the Love numbers given in Eq. 28, we have finally:
Appendix C: First terms of the centrifugal, gravitational, and tidal potentials
C.1. Centrifugal potential in unsynchronised systems
In a nonsynchronised system, the centrifugal potential expressed with the Cartesian coordinates centred on the primary is given by
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^{2} + y^{2} can be written when passing in spherical coordinates as with z = r_{1} cos θ that can be expressed with the spherical harmonics
The centrifugal potential can therefore be expressed as
and decomposed as
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
where r is the instantaneous orbital separation of the stars and P_{ℓ}(sin θ_{1} cos ϕ_{1}) is the classical Legendre polynomial that can be decomposed using spherical harmonics as
where are constant coefficients linked to the normalisation of the spherical harmonics. In the following we define and use
Using this notation:
yield . Similarly,
yields . Therefore,
Summing Eqs. (C.3), (C.11) and using the definition of ñ (Eq. 5), we have:
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
and
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:
that can be rewritten, using Eq. C.6 as
where P_{2} 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
Following Kopal (1959), different coordinates related quantities from the secondary star can be approximated with series developments in the coordinates of the primary star:
and
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:
that can be rewritten using the third Kepler law including the periods of the system as follows
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_{1}/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
when treating the primary as point like (r_{1} = 0). In this case, only the gradient of the ℓ = 1 component is nonnull at the centre of the star (r_{1} = 0) as
Therefore, by only taking the ℓ = 1 component, the force can be expressed as
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
Beyond the surface of the secondary, at a distance r_{2} from its centre, the potential perturbation from the ℓ = 3 perturbed potential of the secondary is expressed as
Based on Kopal (1959), λ_{2} can be expended as
meaning that
and
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:
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
Following the Sect. 3.2, the acceleration perturbation undergone by the system due to the tidal perturbation is expressed as
when defining a as the semimajor axis of the system. Following the same approach as in Sect. 3.3, the time orbital average of the perturbation is formally expressed as
which can be developed into polynomial functions of the eccentricity using Hansen’s functions:
with
The apsidal motion caused by the ℓ = 3 acceleration perturbation can finally be obtained by inserting Eq. (E.12) in Eq. 19 and simplifying:
with,
All Tables
Decomposition of the apsidal motion of the modelled observed systems including the ℓ = 3 term of the perturbative approach.
All Figures
Fig. 1. View of a 1.0 M_{⊙} MS primary star in a binary system, cut in the plane including all the rotations axes of the system. The companion to this star is a 0.2 M_{⊙} MS secondary star with an orbital separation of r = 2R_{1}. The black curve corresponds to the Roche lobe of the primary star. The black lines respectively correspond to the reference direction for θ and to the reference direction in witch our model is given by the 1D spherical symmetric model (θ = μ_{crit}). The blue lines are equipotentials inside the star with the colourbar representing the total potential. On each equipotential the pressure, density, and temperature, are constant. L_{1} is the Lagrangian point and R_{1} is the radius of the initial 1D model in the reference direction. 

In the text 
Fig. 2. Illustration of a binary system configuration with nonsynchronised rotations. The system is composed of 1 M_{⊙} MS primary star and a 0.2 M_{⊙} MS companion separated by r = 2R_{1}. The two grey planes are the planes of symmetry, each rotation axis is marked with a black line. The orbital rotation is originating from the system centre of mass, located inside the primary. The black curve is the secondary orbital path, a black axis is joining the centre of each star with and emphasis on the Lagrangian point L1. Finally, the colour code corresponds to the surface effective gravity, the local flux emitted by a star being a function of this quantity. 

In the text 
Fig. 3. Orbital trajectory of a star in a binary system with an eccentricity of 0.2. The orbital plane is in red, the grey plane is a the reference plane (i.e. the plane of the sky), i is the inclination of the system compared to the reference plane, Ω is the longitude of the ascending node, Ω_{0} is the ascending node, ω is the argument of the periastron, and φ is the true anomaly. The reds dots on the trajectory represent the different points at which the snapshots are taken to compute the apsidal motion. 

In the text 
Fig. 4. Surface potential discrepancy between perturbative and nonperturbative modelling approaches as a function of the orbital separation normalised by the stellar radii for the different binary components presented in Table 1. The upper panel corresponds to the ℓ = 2, m = 0 component of the potential while the lower panel is the ℓ = 2, m = 2 term. 

In the text 
Fig. 5. Perturbed acceleration discrepancy as a function of the separation normalised by the stellar radii for the twin binary systems composed of the stars presented in Table 1. The colour code denotes the different twin binary systems. 

In the text 
Fig. 6. Evolution of the apsidal motion relative difference as a function of the orbital separation normalised by the stellar radii. The eccentricity of the orbits were set to e = 0.1 and the colour code denotes the different twin binary systems. 

In the text 
Fig. 7. Apsidal motion relative different for a twin binary system composed of two 20 M_{⊙} stars presented in Table 1. Each colour correspond to the variation of discrepancy as a function of the orbital separation scaled by the eccentricity and stellar radii for different values of the orbital eccentricity. The black lines are lines of constant model relative difference, corresponding to 5% and 2% discrepancy. 

In the text 
Fig. 8. Mapping of the acceleration perturbation discrepancy for our grid of lowmass stars and intermediatehigh mass stars. The colour code denotes the difference of perturbed acceleration between the perturbative and our nonpeturbative approach as a function of the models k_{2} and orbital separation. 

In the text 
Fig. 9. Apsidal motion discrepancy as a function of the orbital separation scaled with the q parameter. Each system is composed of a primary 1.5 M_{⊙} RGB star presented in Table 1, and a secondary star with a lower mass. The stars composing the systems have the same age (2 Byr) and initial chemical composition (X_{0} = 0.72, Z_{0} = 0.010). Each colour corresponds to a secondary with a different mass. 

In the text 
Fig. 10. Apsidal motion discrepancy as a function of the orbital separation scaled with the q parameter. Each system is composed of a primary 1.5 M_{⊙} MS star, and a secondary star with a lower or equal mass. The stars composing the systems have the same age (1.5 Byr) and initial chemical composition (X_{0} = 0.72, Z_{0} = 0.010). Each colour corresponds to a secondary with a different mass. 

In the text 
Fig. 11. Decomposition of the perturbed acceleration for the IMPer system with the perturbative and nonperturbative models as a function of the system orbital separation normalised by the stellar radii. Each colour represents one spherical order: the lighter curves correspond to the perturbative model while the darker curves correspond to the with nonperturbative approach. The red curves are the total of all the spherical orders. 

In the text 
Fig. 12. Repartition of observed binary systems from our catalogue affected by our methodology. Each point correspond to an observed system with a given orbital separation and a relativistic contribution to the apsidal motion. The blue dots correspond to massive binaries, the yellow points are intermediate mass binary stars, and the red dots are lowmass stars. For this diagram, the primary was chosen as the most massive star in the system. 

In the text 
Fig. A.1. Scheme of a binary system separated by a distance r. The different quantities used to model the deformations in the perturbative framework are detailed and illustrated here. The blue curve is the surface of the deformed star while the red curves are the spheres of equivalent averaged radii. As the surface is represented here, r̅_{i}(μ_{i},ϕ_{i}) = R_{i}, R_{i} being the radial coordinate of the stellar surface in the spherically symmetric case. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.