Underestimation of the tidal force and apsidal motion in close binary systems by the perturbative approach: Comparisons with non-perturbative models

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. Our objective is to assert the reliability of the perturbative approach when applied to close and strongly deformed binary systems. We have developed a non-perturbative 3D modelling method designed to account for high stellar deformations to explore the limitations of the perturbative models. 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 first-order treatment in the perturbative approach, and cannot be rectified using straightforward mathematical corrections due to the strong non-linearity 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)/R1<6.5. The perturbative approach underestimates tidal interactions between bodies up to ~40% for close low-mass 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.


Introduction
Tidal forces in binary systems result from the gravitational interaction between non-point-like bodies and can be categorised into two components: the equilibrium tides, corresponding to 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 corrections1 (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(Rosu et al. , 2022a(Rosu et al. ,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 1959Kopal , 1978)).This approach allows one to obtain simplified structural deformations from the unperturbed spherically symmetric structure through the well-known stellar structure constant, the k 2 coefficient that is obtained from solving the Clairaut-Radau equation (Kopal 1959).With the perturbative formalism, the dipolar components of the tidal and centrifugal forces cancel out (Kopal 1959;Fitzpatrick 2012), implying that no asymmetry exists in the stellar models, in contradiction with the fundamental Roche geometry.Moreover, our previous work (Fellay & Dupret 2023) already demonstrated that the perturbative approach underestimates the deformations when compared to our non-perturbative models.Therefore, it becomes crucial to verify the consequences of the assumptions inherent to the perturbative approach.Finally, it is worth mentioning that the tidal distortion and dissipation formalism (Zahn 1977;Hut 1981) employed by widely recognised and publicly available codes such as MESA (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018(Paxton et al. , 2019) ) are based on the perturbative approach.
In this study, we expand on the work done in Fellay & Dupret (2023) by further investigating the limitations imposed by the assumptions of the perturbative approach, using our non-perturbative method.Our tool, called Modelling Binary Deformations Induced by the Centrifugal and Tidal forces (MoBiDICT), allows us to determine the deformed structure of close binaries in a non-perturbative manner, accounting for the redistribution of mass in the bodies.By utilising the deformed structure obtained from MoBiDICT, we could calculate the instantaneous tidal acceleration perturbation and its consequence on the apsidal motion of binary systems.Finally, MoBiDICT's formalism enabled us to locate and identify the specific spherical order responsible for the discrepancies seen.
In this article, we start, in Sect.2, by introducing the physics behind MoBiDICT and the methodology developed to compute the apsidal motion of binaries in a non-perturbative way.In Sect.3, we explain how we developed the perturbative approach formalism usually used to model the apsidal motion of binaries.Sect. 4 is dedicated to exploring the discrepancies in models for different types of theoretical stars.In Sect. 5 we study binary systems for which the apsidal motion has been measured and provided in the literature.In Sect.6, we discuss the implication of our findings on the stellar modelling methods used for binary stars.Finally, in the Appendices, we develop all the equations of the perturbative approach starting from a single star's gravitational potential to arrive at the apsidal motion in binaries.

Set-up of the problem and first models
In this article, we investigate the deformations in binary systems resulting from three primary sources: the gravitational force, the tidal force arising from a companion, and the centrifugal force originating from the orbital rotation of the system and the individual stellar rotations.Accurate modelling of these phenomena necessitates a non-perturbative treatment of deformations in three dimensions (3D).
We considered a binary system consisting of two distinct stars separated by a distance r and with an orbital rotation rate n.Assuming the stars (i = 1, 2) within the system are in hydrostatic equilibrium, their structures are governed by the equation of hydrostatic equilibrium: 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.
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 non-perturbative approach to treat the deformations.Our method iteratively solves the Poisson equation on a 3D structure, allowing to obtain the gravitational and tidal potential accounting for the redistribution of mass within the bodies.MoBiDICT takes advantage of the conservative nature of all the forces considered, resulting in a barotropic model structure (i.e. the densities are constant on equipotentials).By assuming that in a given direction (µ crit , φ crit ) each star is described by a 1D spherically symmetric input model, we can recompose the entire 3D density profile from the total potential (Ψ tot = Ψ 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 non-perturbative modelling method.In this case, we applied our tool to a binary system composed of a primary 1 M main-sequence (MS) star and a 0.2 M MS companion separated by r = 2R 1 , R 1 being the radius of the primary given by the 1D spherically symmetric model.
In the scenario depicted in Fig. 1, we assumed that the rotational velocities of both stars are identical and equal to the orbital A210, page 2 of 20 rotation rate, in other terms the system is tidally locked.However, when studying eccentric binaries, pseudo-synchronisation becomes the end result of tidal interactions, making impossible to model the system as entirely synchronised.Furthermore, mass transfers between binary components often result in significant transfer of angular momentum (Packet 1981), leaving the accretor star close to its critical rotation while the donor slowly rotates.Generalising our model to non-synchronised binary systems is necessary to model eccentric binaries and binary evolution.

Non-perturbative modelling of non-synchronised binary systems
We considered a scenario where each star i composing a system has its individual rotation rate Ω ,i , independent from the orbital rotation rate n.In this situation, assuming alignment of all rotational axes, two planes of symmetry emerge and are exploited by our method.The first plane of symmetry is the orbital plane, the second is the plane including all the rotation axes.A 3D representation of this configuration is shown in Fig. 2 for the same system as in Fig. 1.With our non-perturbative method, unsynchronised binary systems can be directly modelled by modifying, and generalising the centrifugal potential.The previous expression of the centrifugal potential, which relied on the assumption of solid body synchronised rotation and the circular orbit was and becomes in the non-synchronised 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 major-axis of the orbit) pointing towards the companion, x CM denotes the distance between the system's centre of mass and the star's centre, and y is the coordinate along the perpendicular axis to x included in the orbital plane.
As expected from the comparison of Eqs. ( 3) and ( 4), for a synchronised system (Ω ,i = n), these expressions are equivalent up to a constant factor.The effective gravity resulting from the centrifugal force can be obtained by differentiating Eq. (4).
Our model can be generalised to eccentric systems by decomposing the motion of each element as a Keplerian motion around the system's centre of mass plus a cylindrical rotation of the star.To obtain the centrifugal potential at a given point of the orbit where the instantaneous distance between the centre of the two binary components is r, n 2 has to be replaced, in Eq. ( 4), by ñ2 defined as so that ñ2 r 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 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.
is the gradient of the centrifugal potential that is expressed for the primary star as where the bold symbols denote vectorial quantities, êx 1 is the unit vector linked to the x-axis of the primary and ês 1 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.

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, êr 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 A210, page 3 of 20 secondary on its primary can be decomposed as In Eq. ( 9), F R is the radial perturbed acceleration resulting from considering stars as non-spherical extended interacting bodies.This perturbation is responsible for all the equilibrium tidalinduced dynamical effects of binary evolution, including orbital migration (i.e.variation of systems semi-major axis), eccentricity dissipation (i.e.reduction of orbital eccentricty), and apsidal motion (i.e.time variation of the argument of the periastron).
In the framework of our modelling, the force exerted by the secondary on the primary can be expressed as 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 non-spherical perturbations, respectively ρ 1 and F 2 , 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 non-symmetric terms.In our formalism the perturbed acceleration can thus be obtained with Finally, each force term projected on the spherical harmonic basis F m 2, can be decomposed and re-expressed 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 ∂Ψ m 2, Y ,m /∂r 1 , ∂Ψ m 2, Y ,m /∂θ 1 , and ∂Ψ m 2, Y ,m /∂φ 1 terms can be found in Fellay & Dupret (2023).More than having a non-perturbative description of the resulting force, the main advantage of our formalism is the possibility of naturally compare the contributions of the different spherical harmonic orders by selecting the desired orders in the sums.

Non-perturbative apsidal motion computation
In this section, we developed the formalism used to model the apsidal motion with our non-perturbative method.Figure 3 illustrates the setup of an eccentric binary system, the apsidal motion of a system corresponding to dω/dt.In eccentric binary systems, aspidal motion is a direct observational constraint on the deformations of distorted bodies.Recent investigations (Rosu et al. 2020a(Rosu et al. , 2022a,b) ,b) highlight its utility in constraining physical processes occurring during stellar evolution, such as internal mixing.
The apsidal motion is induced by the equilibrium tidal perturbed acceleration F R , the general relativity and dynamical tidal perturbed acceleration.The general relativistic component of the apsidal motion can be expressed as where a is the system's semi-major axis, e its eccentricity, P orb its orbital period, c is the celerity, and G the gravitational constant.
In this article, we studied the equilibrium tides component of the apsidal motion and neglected the dynamical tides component, the latter will be the focus of a forthcoming article.
From the Gauss planetary equations for a perturbed acceleration along the apsidal line (see e.g.Kopal 1959), the apsidal motion of a binary system is classically given by ωN 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 non-perturbative modelling is the time average perturbed acceleration.To achieve this, we took advantage of the fast computing capabilities of our method, we captured snapshots of stellar deformations at various points along the binary orbits, as depicted in Fig. 3. Using these snapshots, we can derive the timeaveraged perturbed acceleration by reformulating its expression with the second Kepler's law: A210, page 4 of 20 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.
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 non-perturbative 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.

The Clairaut-Radau equation
Perturbative models are obtained by assuming that the deformations of the stellar models are small (more precisely, it is assumed that |r m | 1 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 Clairaut-Radau first-order differential equation (Kopal (1959), and derived in Appendix A) expressed as using the boundary conditions In the Clairaut-Radau equation, r is the average radius of a chosen isobar while ρ is the average density under the chosen isobar.
The relationship between the deformed structure and the η coefficients is explained and derived in Appendix A.

Perturbation of the gravitational potential and force by a companion
We considered the perturbation induced by a point-mass secondary on the primary.As developed in Appendix B, the potential perturbation caused by the secondary at the surface of the primary is given by where Ψ m 2, (R 1 ) 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 Ψ m c, .The first terms of each of these quantities can be found in Appendix C. The coefficients k ,1 also known as Love numbers, characterise the surface non-spheroidal deformation response of a body subjected to a perturbative potential.Mathematically, these coefficients are defined as: Equation ( 27) provides a direct comparison between the models, where the quantity Ψ m 1, (R 1 ) is equivalent to Ψ m 1, (R 1 ) of our non-perturbative iterative method.One particular application to Eq. ( 27) concerns the dominant contribution to deformations and perturbations, the quadrupolar terms ( = 2 terms).For those terms, the sum of the centrifugal and tidal potentials generated at the surface of the primary by the secondary treated as a pointmass body is expressed, in the non-synchronised case, as and In both synchronised and non-synchronised systems, the dipolar component of the total potential is null (Ψ 1 c,1 (R 1 ) + Ψ 1 2,1 (R 1 )) when the secondary is treated as a point-mass, as gravitational and centrifugal forces cancel out according to the third Kepler law.As highlighted earlier, when computing the apsidal motion of a system, it is imperative to account for the tidally induced acceleration perturbation F R .By only considering the spherical A210, page 5 of 20 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 .

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 re-expressed with the zero-order Hansen coefficients X n,m 0 (e) 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, (1 − e 2 ) 5 , (41) and, (42)

General model comparisons
This section is dedicated to a theoretical comparison of different types of stars encompassing a range of potential binary systems across different mass regimes.Here, we present the input stellar models used to compare the perturbative and non-perturbative modelling approaches.We chose four models to represent a broad spectrum of stars and structural characteristics.The physical ingredients of the chosen models are detailed in Table 1, those models are the same models compared in Fellay & Dupret (2023).Table 1 in Fellay & Dupret (2023) has only been extended to include the k 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 low-mass and red giant branch (RGB) stars are the most distorted by a companion while the massive and solar type stars are expected to be less deformed and impacted by our nonperturbative method for a given a/R 1 ratio.In the following sections, we compared the surface deformations and apsidal motion A210, page 6 of 20 obtained with our non-perturbative modelling method and the perturbative modelling for twin binary systems (q = 1) composed of the stars presented in Table 1.Twin binary systems being systems composed of two identical stars (i.e.stars having the same mass, radius, effective temperature, and rotational period).

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 Ψ m i (R i ) with > 0. In the perturbative method, this quantity is given by Eq. ( 27), which is applied to the = 2 terms in Eqs. ( 29) and (30).Within our non-perturbative method, Ψ m i (R i ) 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∆Ψ m (R 1 ) and defined as: where Ψ m ,MoBi is the spectral surface potential obtained with our non-perturbative method and Ψ m ,pert is the same quantity obtained with the perturbative approach.In Fig. 4, we show this surface potential discrepancy as a function of the orbital separation for twin binary systems composed of the stars presented in Table 1.Our focus in this analysis is on the dominant terms, the = 2 terms.The surface potential discrepancy, varies significantly depending on the stellar type and the orbital separation.For both = 2 components, stars with the higher k 2 exhibit more pronounced effects from our non-perturbative treatment, resulting in differences of up to ∼40% in the = 2, m = 2 term for low mass stars.This implies that the higher the deformations, the more substantial the underestimation of the surface deformations by the perturbative approach.Furthermore, orbital separation plays a critical role, with closer star configurations yielding stronger deformations and more notable disparities in the surface potential.
By comparing the two panels in Fig. 4, we observe a significantly smaller difference in surface potential for the = 2, m = 0 term, regardless of the model.The = 2, m = 0 compo-0.000.02 0.04 0.06 0.08 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.
nent arises from both centrifugal and tidal deformations, while the = 2, m = 2 term is solely originating from the tidal deformation.Both our modelling and perturbative modelling are sharing the same assumptions on the centrifugal potential, therefore, the difference originating from tidal contribution in the = 2, m = 0 term is diluted by the centrifugal contribution.The major difference between the perturbative and non-perturbative modelling approaches is thus appearing in the treatment of the tidal deformation.

Apsidal motion
In Sect.2.4 we showed that the apsidal motion of an eccentric binary system is proportional to the time average of the perturbed acceleration.As detailed in Sect.4.1, the perturbed acceleration denotes the tidal acceleration originating from a nearby companion.A modification of the surface = 2 potential as illustrated in Fig. 4 introduces variations in the perturbed acceleration.In this section, we investigated the perturbed acceleration discrepancy resulting from the potential differences seen and the consequent impact on the apsidal motion of eccentric binaries.
In the classical perturbative model, the perturbed acceleration F R is solely originating from the quadrupolar component of the deformation.The detailed expression of F R is given in Eq. ( 33) and derived in Appendix D. In opposition, our methodology incorporates a non-perturbative treatment of all spherical orders spanning = 1, 2, ...L.The exact expression of F R is given by Eq. ( 16).To quantify the disparities between perturbative and non-perturbative perturbed accelerations, we define their difference ∆F R as follows: The difference of perturbed accelerations as a function of the separation between the twin binary systems components presented in Table 1 is illustrated in Fig. 5.The perturbed acceleration discrepancy observed in Fig. 5 is in agreement with the results showed in Fig. 4. Models A210, page 7 of 20  1.The colour code denotes the different twin binary systems.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.Solar-like stars and massive stars are comparatively less affected, even if in close orbit F R discrepancies can still reach up to, respectively, 20% and 12%.
To assess the discrepancy of apsidal motion resulting from the perturbed acceleration differences, we need to consider binaries with eccentric orbits.To limit the free parameters of the system we compare binary systems with synchronised rotations and a fixed eccentricity, e = 0.1.The impact of the eccentricity on the apsidal motion is explored for systems in Sect.4.3.We introduce the apsidal motion difference: 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 non-perturbative method we required to capture snapshots of the binary system along half of an orbit to estimate the perturbed acceleration time average and hence the apsidal motion (see Sect. 2.4).For each binary system composed of the stars presented in Table 1, Fig. 6 illustrates ∆ ω as the function of the binaries orbital separation normalised by the stellar radii.
The apsidal motion discrepancy is a direct result of the perturbed acceleration illustrated in Fig. 5. Binary systems with greater F R discrepancis exhibit the highest apsidal motion differences.For the low mass and RGB stars, the apsidal motion discrepancy can reach up to respectively 70% and 45% when the stars are close to contact, in their periastron.Similarly, for the solar type and massive stars we found an apsidal motion discrepancy up to respectively 30% and 15%.These discrepancies vary depending on the exact structure of the stars or the architecture and the orbital parameters of the system.

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.
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.

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 .
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 intermediate-high 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 non-linear behaviours of the models greatly reduces its feasibility.We attempted an MCMC analysis using linear combinations of power laws of the orbital separation, for both the intermediate-high mass stars and low-mass stars.However, the results were unsatisfactory to propose an empirical parametrisation.As detailed in Sect.5.2, the modelling discrepancies indeed arise from a combination of discrepancies in the = 1, 2, 3 terms.Non-perturbative modelling is necessary to recover the discrepancies seen in this section.However, Fig. 8 can be used to get an order of magnitude of the underestimation of tidal forces introduced by the perturbative modelling.The updated grid of k 2 values provided by Claret (2023) could for example be used to easily locate observed stellar systems in this diagram.

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 non-pertubative 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 sub-product of the observations of binaries.To study and modify the q parameter in our modelling, different stellar models with the same age and initial chemical properties have to be used for the secondary star, impacting 4 of the parameters of our modelling for a given orbital separation: (n, (R 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.
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 x-axis of Fig. 9 to be q −1/5 R 1 /a to maintain the scale between (R 1 /a) 5 and q.
We can see that all the systems are superposed in Fig. 9 which indicates a clear proportionality of the apsidal motion discrepancy on q.When R 5 1 R 5 2 , 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 A210, page 9 of 20 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 q 1/5 a/R 1 0.0 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.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.
proportional to q for a given orbital separation, the apsidal motion discrepancy is also proportional to q.
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.
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 mod-ified 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 R 5 1 R 5 2 however, as showed before, the dependency on q still dominates.

Modelling of observed binaries and quantification of the extra mixing
In this section, we apply our non-perturbative modelling technique to existing twin binary systems, comparing the stellar parameters obtained with our approach to those from the perturbative theory.As we study twin binary systems, modelling one star is sufficient to describe the entire system.We integrated MoBiDICT to a minimisation method following a Levenberg-Marquardt algorithm to compute on the fly the non-perturbative apsidal motion resulting from twin binaries which unperturbed structure are given by our 1D single star stellar evolution code, CLES.By using this methodology we followed the procedure of Rosu et al. (2020aRosu et al. ( , 2022aRosu et al. ( ,b, 2023)), neglecting the interactions of the binaries during their evolution.To avoid such assumptions, our modelling method would have to be directly coupled to binary stellar evolution codes, this work will be conducted in the future with the code BINSTAR (Siess et al. 2013) and will be presented in a forthcoming article.We selected four well-known close observed twin binary systems from the literature: PV Cas (Torres et al. 2010;Claret et al. 2021;Marcussen & Albrecht 2022), IM Per (Lacy et al. 2015;Claret et al. 2021), Y Cyg (Gimenez et al. 1987;Harmanec et al. 2014;Claret et al. 2021;Marcussen & Albrecht 2022), and HD152248 (Rosu et al. 2020a,b).The first binary system is composed of MS intermediate mass stars, the second one of sub-giant stars, and the last two ones of massive stars.The properties of these systems are given in the second and fifth line blocks of Table 2.For each of these binary systems, an apsidal motion was determined in the literature and is used in this section as a constraint on the stellar structure.The systems were selected to be twin systems with accurate determination of the orbital and stellar parameters, and an apsidal motion primarily driven by tidal forces.
For this modelling we were inspired by the findings of Rosu et al. (2020aRosu et al. ( , 2022aRosu et al. ( ,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. 2020aRosu et al. , 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 (τ) A210, page 10 of 20 relation from Model-C of Vernazza et al. (1981) for the atmosphere, the mixing length theory of convection implemented as in Cox & Giuli (1968), and the nuclear reaction rates of Adelberger et al. (2011).For the massive stars, the mass loss rate was computed, accounting for the metallicity of the models, with the prescription of Vink et al. (2001) assuming that the scaling factor for the mass loss rate equals ξ = 1.The equivalent averaged mass loss along the entire evolution is given in the fourth line block of Table 2.
In our modelling procedure, we started by finding models reproducing the observations with the perturbative approach.We then fixed the initial mass and hydrogen mass fraction, only using the constraints on the radii and apsidal motion to model the systems with MoBiDICT.This approach aims to directly evaluate the impact of our non-perturbative modelling on the required additional mixing and age.The objective of our modelling is not to provide an accurate estimation of the stellar properties of each system but rather to compare the stellar properties obtained when modelling the apsidal motion with the perturbative and nonperturbative approaches.The results from our modelling procedure are given both for the non-perturbative and perturbative modelling approaches in the third line block of Table 2.
For the modelling of each system, the parameters presented in the fourth line block of Table 2 were empirically selected to facilitate the convergence of our Levenberg-Marquardt method.For the first system, Z 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 non-perturbative approaches can accurately reproduce all observed system constraints.For the four systems, we find that extremely large extra-mixing is required to reproduce the observational properties of stars.Comparing the stellar parameters obtained with MoBiDICT and the perturbative approach, our non-perturbative technique necessitates even more extra mixing for the first three systems.While for the three first systems, the leading contribution to the χ 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 non-perturbative 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).

Apsidal motion contribution from the different sources
In this section, we focus on exploring the origin of the apsidal motion discrepancy seen for the observed targets presented in Sect. 5. We reused the structures from the modelling done in Table 2 and looked at the different contributions to the apsidal motion for each binary system.Table 3 provides a decomposition of the apsidal motion for each observed binary system using the models obtained with our non-perturbative method.In the right part of Table 3, we also decomposed the difference between the perturbative and non-perturbative approaches based on the Fig. 11.Decomposition of the perturbed acceleration for the IM-Per system with the perturbative and non-perturbative 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 non-perturbative approach.The red curves are the total of all the spherical orders.considered spherical orders, revealing the key contributors to the apsidal motion discrepancies.
The analysis of Table 3 reveals that, for all chosen targets, the contributions from our non-perturbative modelling approach surpass the uncertainties in the observational apsidal motion determinations.Generally, the corrections provided by MoBiDICT are comparable in magnitude to the general relativistic corrections.Then, going into more details, the perturbed acceleration discrepancy between the non-perturbative and perturbative methods can be attributed to several main components: -= 1: Constituting ∼23% of the modelling discrepancy, this component cannot be recovered by the perturbative model due to its inherent formalism and assumptions; -= 2: Representing ∼40% of the modelling discrepancy, this component cannot either be recovered by the perturbative approach.The differences arise from the first-order assumption of the perturbative approach.We tried taking into account the effect from the quadrupolar deformation of the primary star on the quadrupolar component of its companion and the reverse impact of this modification on the quadrupolar deformation of the primary star.However this effect is negligible as each of such additional effect considered has an increase dependency to (R/a) 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 non-perturbative models.

Targets
PV Cas (1,2,3)  IM Per (2,4)  Y Cyg (5,2,3) HD152248 (6,7)   Pert.Notes.The first line block corresponds to the observational constraints used in our modelling, the second line block gives the values of the free parameters obtained, the third line block gives the assumed parameters of the models, and the fourth, provides the assumed parameters of the system originating from observations.All the models were compared using χ 2 defined as the sum of the squared relatives differences between the observed constraints and the models. References.

Higher order perturbative method
An advantage of the perturbative approach lies in the possibility to extend the modelling to higher spherical orders.In the previous section we showed that the = 3 component of the acceleration perturbation represents ∼35% of the discrepancies seen between our non-perturbative and the perturbative models.
In Appendix E, we developed the equations necessary to extend the perturbative approach to the = 3 term.In this section, we verify the relevance of adding the = 3 term in the computation of the apsidal motion in the perturbative approach.As shown in Appendix E, the apsidal motion arising from the = 3 component of the perturbed acceleration is given by (1 − e 2 ) 7 .
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.
References.By comparing Tables 3 and 4 we conclude that by including the = 3 component of the apsidal motion we corrected ∼20% of the discrepancies seen in Sect.5.2.However, the = 3 term is only partially corrected when included in the perturbative approach.Indeed, we found that for all the systems modelled this term stills represent ∼10% of the observed model discrepancies.Due to the reduction of the contribution from the = 3 term, the quadrupolar term now represents ∼50% of the observed discrepancies while the dipolar term contributes to ∼30% of the differences.The exact repartition of these contributions can vary depending on the orbital separation and the stellar type, in particular, we saw a highly non-linear comportment of the dipolar term.Nonetheless, in the future, the = 3 contribution to the apsidal motion have to be included when modelling the apsidal motion of close binaries with the perturbative approach.

Discussion
In the previous sections, we showed that non-pertubative modelling is necessary for a precise interpretation of observations in close binary systems.Aside from the orbital evolution considerations, the detection of tight systems with apsidal motion is favoured due to the (R/a) 5 dependency of the tidal component.The combination of high-precision observations and strong tidal effects on apsidal motion motivates their investigation as stellar laboratories and, consequently, the need to model them comprehensively.Our modelling has conclusively shown that, for close binary systems, a non-perturbative treatment of deformations is necessary for accurately characterising stellar global properties and stellar structure.
To assess the impact of our modelling method on observed binary systems with apsidal motion, we constructed a comprehensive catalogue of such systems.In total, our catalogue includes 61 systems, encompassing all systems with constraints on their stellar parameters and observed apsidal motion found in the literature.These systems were obtained by combining two recent catalogues of observed binary systems with apsidal motion (Claret et al. 2021;Marcussen & Albrecht 2022) to systems from various literature sources where stellar parameter determinations were available (Gimenez et al. 1987;Benvenuto et al. 2002;Wolf et al. 2006Wolf et al. , 2008Wolf et al. , 2010;;Torres et al. 2010;Pablo et al. 2015;Baroch et al. 2021Baroch et al. , 2022;;Rosu et al. 2020bRosu et al. , 2022aRosu et al. ,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.
The observed binary systems in our catalogue can be broadly categorised into two equal-sized groups: close binaries with most of the apsidal motion driven by tidal effects and wider binaries where apsidal motion is mostly due to relativistic corrections.The boundary between these two groups is located around q −1/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 non-perturbative method.On the other hand, for binary systems where apsidal motion is dominated by tidal effects (close systems), which is the group preferred for drawing constraints on stellar structure with apsidal motion, our non-perturbative models show discrepancies compared to perturbative models.In Sect.4.4 we showed that, for a given orbital separation, intermediate-high mass stars are the least impacted by our 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 intermediate-mass binaries.In these catalogues and in general, low-mass binaries are under-represented, and not commonly observed with apsidal motion probably due to observational biases.
The discrepancies observed in the apsidal motion are a consequence of the underlying tidal force discrepancy.However, various aspects and phenomena of binary systems' evolutions are impacted by an underestimation of the tidal force.One key quantity affected by this underestimation is tidal dissipation, which plays a crucial role in the evolution of binary systems through orbital circularisation, transfer of angular momentum between bodies, and orbital migrations.Current models of the tidal dissipation rely on the perturbative approach (Hut 1981) combined to a free parameter that can be adjusted to control the dissipation.Even if this parameter can be tuned to reproduce the dissipation obtained by our non-pertrubative models with enhanced tidal interactions at a given point of the evolution, it cannot capture the non-linearities that arise during the evolution of binary systems.Consequently, our methodology predicts that tidal dissipation models used by the majority of binary stellar evolution codes A210, page 13 of 20 Fellay, L., et al.: A&A, 683, A210 (2024) 10 3 4 5 6 7 8 9 20 30 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 low-mass stars.For this diagram, the primary was chosen as the most massive star in the system.
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 non-perturbative methodology.To precisely capture the interior deformations, mass redistributions, and structural properties of deformed stars, a transition to 3D modelling becomes imperative.Having highly accurate structural models is crucial for the precise study and modelling of the pulsation properties of deformed celestial bodies.Among the oscillations propagating in deformed stars, gravito-inertial waves, called dynamical tides, significantly contribute to the orbital evolution of binary systems.However, the current dynamical tides models rely on structural modelling techniques based on the perturbative approach combined with 1D non-adiabatic spherically symmetric oscillations codes to treat this highly sensitive and non-linear problem.As such, the dynamical tide models, as obtained in Ma & Fuller (2023), for example, are only estimations of the impact of dynamical tides on close binary systems.Without having appropriate 3D non-perturbative deformed structures and the associated 3D non-perturbative non-adiabatic oscillation codes this problem cannot be precisely solved.
Our methodology, similar to the perturbative approach, does have a notable limitation: it assumes solid body rotation within stars.For single stars, this assumption is not fully justified because they often exhibit strong differential rotation, which varies over the course of stellar evolution.However, as we showed in Sect.5.1, interactions between binary components induce significant extra mixing within deformed stars.This additional mixing reduces the differential rotation in interacting bodies, satisfying our assumption.The subject of angular momentum transport processes in binary systems is a complex and extensive area of research that warrants further investigation.MoBiDICT represents an important initial step in this direction.It allows for the first precise 3D computation of stellar interactions and the resulting forces at each point within the deformed stellar structure.
Finally, the modelling of mass transfer and phases of common envelopes in binary systems is still associated with significant uncertainties and often relies on the use of multiple free parameters.In theory, MoBiDICT has the potential to be adapted for the static modelling of common envelope phases or, at the very least, to provide initial 3D deformed structures that can be used in conjunction with more sophisticated codes to simulate common envelope evolution.Furthermore, our 3D non-perturbative models can be valuable in modelling mass transfer processes in binaries or providing initial 3D deformed structures for further studies.In the near future, it would be useful to investigate the oscilation properties of binaries before they fill their Roche lobes as determining their configurations and properties before these phases play a crucial role in determining their ultimate fate.

Conclusion
In this work, we investigated the impact of non-perturbative modelling on the theoretical computation of the apsidal motion and tidal force in eccentric binaries.We developed formalisms for apsidal motion in our non-perturbative approach and the perturbative approach in respectively Sects. 2 and 3.Then, we compared the perturbed acceleration and apsidal motion for theoretical models, in Sect.4, and observed binary systems, in Sect. 5. Finally, we discussed the implications of our results, in Sect.6, on the observed population of binaries and the general problem of binary modelling.
All theoretical stellar models were significantly impacted at low orbital separation by our methodology.For instance, we observed a maximum increase in tidal force of about 40% for our low-mass stars, 28% for the RGB, 20% for the solar types, and 12% for massive stars.Consequently, the apsidal motion of these systems is impacted.For an orbital eccentricity of e = 0.1, these discrepancies led to a respectively increase in apsidal motion of 70%, 45%, 30%, and 15% for the low-mass, RGB, solar type, and massive twin binary systems.The discrepancies seen may vary depending on the exact structure of the stars or the architecture and the orbital parameters of the system.These results are consistent with Fellay & Dupret (2023): the models with higher envelope mass experience greater impacts from our methodology.From this study we noted that models with higher k 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 A210, page 14 of 20 stars have similar radii.In the latter case, the properties of the secondary are modifying this relationship and the exact properties of the system have to be accounted for.
In our analysis of observed systems, we explored in more details the origins of these discrepancies.Specifically, we modelled four observed twin binary systems, PV Cas, IM Per, Y Cyg, and HD152248, composed of respectively, intermediate mass, sub-giants, and massive stars for the two latests.We started by modelling these systems using the pertubative model to reproduce the observed values of the constraints parameters.When best fitting models were found, we started from this solution to model stars using our non-perturbative method to obtain the theoretical apsidal motion.In agreement with the results of Rosu et al. (2020aRosu et al. ( , 2022aRosu et al. ( ,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. (2020aRosu et al. ( , 2022aRosu et al. ( ,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 non-negligible 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 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 non-perturbative non-addiabatic 3D astreroseimic modelling is required to study the propagation of gravito-inertial waves in the stellar interior and their impact on the dynamics of binary systems.In the future, our aim is to develop such a tool.In addition, in the future, the inclination of binary systems have to be considered when modelling the tidal contribution to the apsidal motion to account for all the parameters of the orbits.An extension of MoBiDICT for non-aligned systems can be developed; however, all the system's symmetries exploited in this work will disappear, consequently our methodology will require significantly more resources to run.
Importantly, the underestimations seen extends beyond the apsidal motion to encompass tidal interactions between the binary components, thereby impacting tidal dissipation pro-cesses.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.be written when passing in spherical coordinates as r 2 1 − z 2 with z = r 1 cos θ that can be expressed with the spherical harmonics The centrifugal potential can therefore be expressed as and decomposed as Ψ m c, = 0 for all other combinaisons of andm.

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 (C.13) 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: Ψ 2,2 (R 2 ) = Ψ 0 2,2 (R 2 )Y 2 0 (θ 2 , φ 2 ) + Ψ 2 2,2 (R 2 )Y 2 2 (θ 2 , φ 2 ) (D.1) that can be rewritten, using Eq.C.6 as

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 Table1.The colour code denotes the different twin binary systems.

Fig. 6 .
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.

Fig. 7 .
Fig.7.Apsidal motion relative different for a twin binary system composed of two 20 M stars presented in Table1.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.

Fig. 8 .
Fig. 8. Mapping of the acceleration perturbation discrepancy for our grid of low-mass stars and intermediate-high mass stars.The colour code denotes the difference of perturbed acceleration between the perturbative and our non-peturbative approach as a function of the models k 2 and orbital separation.
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 Table1, 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.
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.

Table 1 .
Summary of the stellar properties of the 1D input models used in this work.

Table 2 .
Observed and model properties of the systems studied.

Table 3 .
Decomposition of the apsidal motion of the modelled observed systems.

Table 4 .
Decomposition of the apsidal motion of the modelled observed systems including the = 3 term of the perturbative approach.