Issue 
A&A
Volume 676, August 2023



Article Number  A22  
Number of page(s)  14  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202346175  
Published online  28 July 2023 
MoBiDICT: New 3D static models of close synchronised binaries in hydrostatic equilibrium
STAR Institute, University of Liège, 19C Allée du 6 Août, 4000 Liège, Belgium
email: lfellay@uliege.be
Received:
17
February
2023
Accepted:
25
May
2023
Context. In close binary systems, tidal interactions and rotational effects can strongly influence stellar evolution as a result of mass transfer and common envelope phases. These aspects can only be treated following improvements of theoretical models to take into account the breaking of spherical symmetry occurring in close binaries. Current models of binary stars rely on either the socalled Roche model or the perturbative approach, both of which result in several assumptions concerning the gravitational, tidal, and centrifugal potentials.
Aims. Our aim is to develop a precise 3D model of stellar deformations and to study the robustness of the Roche and perturbative models in different deformation regimes.
Methods. We developed a new nonperturbative method to compute the precise structural deformation of a binary system in three dimensions that is valid even in the most distorted cases. We then compared our new method to the Roche and perturbative models for different orbital separations and binary components.
Results. We found that in the most distorted cases, both the Roche and perturbative models significantly underestimate the deformation of binaries. The effective gravity and the overall structural deformations are also noticeably different in the most distorted cases, leading to modifications of the usual gravity darkening generally obtained through the Roche model when interpreting the observations. Moreover, we found that the dipolar term of the gravitational potential, usually neglected by the perturbative theory, has the same order of magnitude as the leading tidal term in the most distorted cases.
Conclusions. We developed a new method that is capable of precisely computing the deformations of a binary system composed of any type of stars, even compact objects. For all the stars we studied, the differences in deformation with respect to the Roche or perturbative models are significant in the most distorted cases, which impacts both the interpretation of observations and the theoretical structural depiction of these distorted bodies. In the weaker deformation regimes, the Roche model is a viable option for studying the surface properties of binaries, while the perturbative model is strongly favoured when evaluating structural deformations.
Key words: binaries: close / binaries: general / stars: interiors / stars: evolution / Celestial mechanics
© The Authors 2023
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
The recent work of Sana et al. (2012, 2013, 2014) showed that a significant fraction of massive binary stars (> 70%) are impacted during their lifetime by interactions with their companions. While massive stars have a high occurrence rate in a binary system (Sana & Evans 2011; Sana et al. 2012, 2013, 2014), both evolved and main sequence (MS) lowintermediate mass stars are also commonly found in binary systems (Slawson et al. 2011; PriceWhelan et al. 2018; Gaia Collaboration 2022). The recent Gaia Data Release 3 (Gaia Collaboration 2021) offers the possibility to catalogue nonsingle stars, as more than 800 000 of them were discovered (Gaia Collaboration 2022), allowing us to find and study the most peculiar systems.
In a binary system, the presence of a nearby companion alters the classical evolutionary path of each binary component (Hurley et al. 2002; Sana et al. 2012), leading in the most extreme cases to Xray binaries (Savonije 1978; Patterson 1984) or gamma ray bursts (Aharonian et al. 2005), for example. For any type of close binaries, the Rochelobe overflow is a common phenomenon (Sana & Evans 2011; Sana et al. 2012, 2013, 2014). It is characterised by a material and angular momentum transfer that can even strip stars from their envelope in some systems, and it is believed to be at the origin of hot subdwarfs (Han et al. 2002, 2003; Heber 2009). Moreover, phases of common envelope can allow the binary components to exchange angular momentum or chemicals, for example. Before reaching such stages, stars experience deformations induced by their mutual interactions dominated by the gravitational, tidal, and centrifugal forces.
A correct modelling of these deformations is crucial in order to determine an accurate stellar structure of each component of the binary system during these phases and to study their evolution. Such models are also necessary as starting points to study the physics behind mass transfer in more detail, which currently relies on numerous free parameters and assumptions. Finally, these models can contribute to precise computations of system apsidal motions that impose direct observable constraints on the stellar structure in a binary system usually theoretically obtained with the perturbative model (Rosu et al. 2020).
The modelling of deformations is currently based on two independent types of models: the Roche model and the perturbative model. The principle of the Roche model is to treat binary components as pointmass bodies using a simplified expression of the gravitational and tidal potentials. On the other hand, the principle of the perturbative model is to treat the centrifugal and tidal forces as a first order perturbation to the spherical symmetry, only accounting for the leading orders of gravitational potential (ℓ = 0, 2 Kopal 1959, 1978) and neglectingthe dipolar term (ℓ = 1) in particular. Both models are first order approximations that consider spherically symmetric stellar density distributions. In the least distorted cases, the different assumptions of the classical modelling procedure are justified and allow for the study of the stellar surface deformations, effective gravity, and even structure (Siess et al. 2013; Palate et al. 2013; Prša 2018; Fabry et al. 2022). However, in the most distorted cases, the redistribution of mass caused by the deformation, as well as the different assumptions made by the classical models, can cause inaccuracies in the depiction and study of the stellar deformations, particularly in the lowmass regime (Landin et al. 2009).
To study the limits of classical models and get an accurate computation of the stellar surface deformations, even in the most distorted cases, we developed a new nonperturbative approach in three dimensions based on a method presented in Roxburgh (2004, 2006) and recently used by Manchon (2021), Houdayer & Reese (2023) to model rotating stars in two dimensions. This tool, called Modelling Binary Deformation Induced by Centrifugal and Tidal forces (MoBiDICT), was developed in Fortran 95 for computational efficiency. The principle of MoBiDICT is to iteratively solve Poisson’s equation, including the centrifugal and tidal forces in a nonperturbative way. More than obtaining accurate surface quantities, such nonperturbative models give access to the entire, precise 3D structure of each binary system component, allowing accurate computation of gravity darkening or the apsidal motion constant, even in the most distorted cases. Finally, the impact of deformations on the time evolution of the stellar internal structure becomes accessible with such detailed models.
In this work, we start by giving a technical and physical description of our new modelling method in Sect. 2. In Sect. 3, we present the posttreatment of our code that computes the effective gravity and gives the physical insights to couple MoBiDICT to classical 1D stellar evolution codes. Section 4 is dedicated to the presentation of the Roche model and its comparison with our new models. The presentation and comparison with the perturbative model is done in Sect. 5. Finally, the conclusion, the possible improvements, and applications of our new types of models are discussed in Sect. 6.
2. Stellar models and properties
In this section, we describe and give the technical details of our novel modelling procedure to compute 3D nonperturbative models of close synchronised binaries in hydrostatic equilibrium. Our method is an iterative procedure designed to include the impact of stellar deformations on the stellar structure of the whole system.
In Sect. 2.1, we present the general properties of the problem that we solve as well as the general principle of MoBiDICT. Then in Sects. 2.2 and 2.3, we detail the modelling procedure by describing the different steps of each iteration. Finally, in Sect. 2.4, we present our convergence criterion and the different readjustments made at the end of each modelling iteration.
2.1. General modelling procedure
In the case of a binary system, the equations governing the hydrostatic structure of a given rotating component i are given by
and the Poisson equation
where ρ(r, μ, ϕ) and P(r, μ, ϕ) respectively denote the density and the pressure in the spherical polar coordinates (r, μ, ϕ), with μ = cos(θ), θ being the polar angle, ϕ being the azimuthal angle, and r being the radial coordinate. The gravitational potential of the arbitrary selected primary star of the system is denoted by Ψ_{1}(r, μ, ϕ); Ψ_{2}(r, μ, ϕ) is the gravitational potential of its companion; and Ψ_{centri} is the centrifugal potential later defined in Eq. (13). Equation (2) is a second order nonlinear differential equation in three dimensions with Ψ_{i}(r, μ, ϕ)→0 at the infinity as the boundary conditions.
Assuming stars with a solid body rotation, all the forces here considered were derived from a potential. Thus the problem is conservative. By combining Eq. (1) and its curl,
we deduced that the structure of each star is barotropic, which means that the density and pressure are constant on a given equipotential. At a given chemical composition, the structure of stars becomes thermotropic, implying that the temperature is also constant on the equipotentials. To take advantage of these properties, we assumed that the density of each star is given by a one dimensional, spherically symmetric input model along a given arbitrary direction (μ_{crit}, ϕ_{crit}). Therefore, in this direction, we can associate a total potential to a density distribution, and by computing the total potential on the whole star, we deduce the 3D density profile of each component of the system. The choice of the direction (μ_{crit}, ϕ_{crit}) is detailed in Sect. 3.1.
By assuming that stars have their rotation axis aligned to the orbital rotational axis, two planes of symmetry appeared. The first one is the orbital plane, and the second one is in the plane including the two stellar rotation axes and the orbital rotation axis. With these two symmetries, a description of one quarter of each star was sufficient to model the entire system in 3D.
The principle of our modelling method is to iteratively solve Poisson’s equation on a spherical harmonic basis for each star in order to compute, in a nonperturbative way, the impact of tidal and centrifugal forces on the stellar hydrostatic structure. This method is a generalisation of the method of Roxburgh (2004, 2006) in 3D and for binary systems. Each iteration of our modelling technique proceeded as follows:

As an initial model, we took an input density profile from a spherically symmetric 1D model (i.e., ρ_{i}(r, μ, ϕ) = ρ_{i, 1D}(r)) obtained through a stellar evolution code for each star i composing the system.

We solved Poisson’s equation (Eq. (2)) knowing ρ_{i}(r, μ, ϕ), and we computed Ψ_{tot} = Ψ_{1} + Ψ_{2} + Ψ_{centri} for each star.

We deduced ρ_{i}(r, μ, ϕ) from Ψ_{tot}(r, μ, ϕ) using the barotropic property of our problem. More precisely, along the direction (μ_{crit}, ϕ_{crit}), we computed Ψ_{tot}(r, μ_{crit}, ϕ_{crit}) and imposed that ρ_{i}(r, μ_{crit}, ϕ_{crit}) = ρ_{i, 1D}(r). Since ρ_{i} and Ψ_{tot} are two monotonic functions of r, we obtained the function ρ_{i}(Ψ_{tot}). Finally, we simply had ρ_{i}(r, μ, ϕ) = ρ_{i}(Ψ_{tot}(r, μ, ϕ)).

We estimated the quantities δρ_{i}(r, μ, ϕ) and δΨ_{i}(r, μ, ϕ) that we used as a convergence indicator and readjusted the model. If the model did not converge, we started again from step two, using ρ_{i}(r, μ, ϕ) obtained in step three as an input.
2.2. Solution for Ψ_{tot} given ρ_{i}
Given ρ_{i}(r, μ, ϕ) on a grid (r, μ, ϕ) for each star i, we solved Poisson’s equation (Eq. (2)) in order to obtain the potential Ψ_{i}(r, μ, ϕ) and deduce Ψ_{tot} = Ψ_{1} + Ψ_{2} + Ψ_{centri}. We decided to adopt an individual spherical coordinates grid (r_{i}, μ, ϕ) centred on each star i with r_{i} pointing toward the centre of mass of the system when ϕ = 0 and μ = 0. The colatitude μ and ϕ share a common mesh for each star, while r_{i} depends on the mesh of the input 1D models. Taking advantage of the two symmetrical plans of the problem, we could limit the angular domain to μ ∈ [0, 1] and ϕ ∈ [0, π].
Generalising the method presented by Roxburgh (2006) in 3D and for binaries, we expressed ρ and Ψ as a finite sum of spherical harmonics with
and
where p = 0 if ℓ is even, p = 1 if otherwise, and m = 2k + p. The free parameter of our model defining the maximum degree of the spherical harmonics to account for is as denoted L. In our work, the spherical harmonics were normalised as
where δ_{ij} is Kronecker’s delta, only ≠0 and = 1 when i = j. To obtain , we solved Poisson’s equation on the spherical harmonics basis, which is expressed as
with the boundary conditions as if ℓ ≠ 0 and
where R_{i, 0} denotes a given radius at the exterior of the star i. To project the density on the spherical harmonics basis, we used the integral relation
As mentioned by Roxburgh (2006), the computation of through Eq. (9) using a detailed mesh can become a real problem in terms of computation time and the precision required. To reduce this problem, we decided to use a mesh of μ and ϕ following the GaussLegendre quadrature method. This method allows for an exact solution to an integral of a polynomial of degree 2n − 1, where n is the number of points in the integration domain. Drastically reducing the number of points in grids while increasing the integral precision, the GaussLegendre quadrature method is particularly efficient for computing integrals of functions, such as our spherical harmonics.
An integral representation of the solution of Eq. (7) with the boundary conditions presented in Eq. (9) can be written as
An equivalent representation of Eq. (10) used by Kopal (1959, 1978), Fitzpatrick (2012), for example, is expressed as
Both equivalent solutions were implemented in our code, and the only difference between the two is related to the numerical treatment of the integrals in the limit r_{i} → 0. To avoid further complications in this region, we exploited the linearity of close to the singularity with our trapezoidal integration method to increase the integration precision. With both solutions, the same results were obtained. For all the work presented in the rest of this article, we used the solution given by Eq. (10).
Inside MoBiDICT, all the equations solved and quantities used are dimensionless in order to facilitate the computations. To nondimensionalise the system, we used constants of a binary system. The different dimensionless variables of the system are expressed as
where a is the separation of the binaries, M_{tot} is the total mass of the binary system, Ω_{⋆} is the orbital rotation rate, and G is the gravitational constant. From a numerical point of view, having dimensionless variables is useful. In our case, the quantities used to nondimensionalise were chosen to simplify the 4πG constant in Poisson’s equation (Eq. (2)) and the integral representation of its solution (Eqs. (10) and (11)). In the rest of this article, the equations and quantities presented are with dimensions to maintain the physical sense of the problem.
Using the presented method, the following steps are performed to compute the total potential of each star:

With the density provided for each star ρ_{i}(r_{i}, μ, ϕ), we determined their spectral expansion through Eq. (9).

We numerically computed Poisson’s equation (Eq. (7)) solution given by Eq. (10) for each star in order to obtain their gravitational potential projected on the spherical harmonics basis: and .

Using the proper coordinates conversion and the previously determined spectral quantities, namely and , we determined each gravitational potential on the grid of each star, obtaining Ψ_{1}(r_{1}, μ, ϕ), Ψ_{1}(r_{2}, μ, ϕ), Ψ_{2}(r_{1}, μ, ϕ), Ψ_{2}(r_{2}, μ, ϕ). Finally, we computed the total potential Ψ_{tot} = Ψ_{1} + Ψ_{2} + Ψ_{centri} on the grid of each star.
The centrifugal potential of each star, assuming a synchronised solid body rotation aligned to the orbital rotation axis, can be expressed as
where ϖ(r, μ, ϕ) is the distance, on the orbital plane from the orbital rotation axis of the element considered, Ω_{⋆} is the stellar and orbital rotation rate, x_{CM} is the distance of the star from the centre of mass of the system.
Our formalism is not limited to solid body rotation with aligned rotational axis. By modifying the centrifugal potential, we can extend the method to model the cases of nonsynchronised, nonaligned solid body rotations. These specific cases could be important to model since we expect nonsynchronised stellar rotations close to the critical velocity as a result of mass transfer within binary systems. This generalisation of our formalism will be presented in a forthcoming article.
2.3. Solution for ρ given Ψ_{tot}
Given Ψ_{tot}, which we determined in the Sect. 2.2, we wanted to estimate the corresponding density for of each star. We exploited the fact that the structure of each star is barotropic, meaning that the density is constant on the equipotentials.
Assuming that the density of a star i along a given direction μ_{crit} and ϕ_{crit} is the density of the corresponding averaged 1D input model, we could determine the functions ρ_{i}(Ψ_{i, tot}(r_{i}, μ_{crit}, ϕ_{crit})) = ρ_{i, 1D}(r_{i}) for each star. Then, by simple function composition and interpolation, we could estimate the density of the entire star with the total potential on each point of the grid.
2.4. Readjustment and criterion for the convergence of the method
A readjustment of the density was necessary to conserve the total mass of each star and the system. By the simple integration of we obtained the total mass of each star with our new density profiles. We then corrected the entire stellar density profiles by the ratio of the mass from the average 1D input models to the new mass obtained.
At this modelling step, several readjustments were made in order to verify different basic laws impacted by the nonsphericity of the models. In particular, we solved a system of equations verifying the balance of forces at the centre of each star in order to properly recompute the centre of mass and the orbital period of the system and maintain the balance of the forces. The system of equations solved at this modelling step is presented in Appendix A.
The convergence of our method was guided by the variations of the density and gravitational potential between two successive iterations. These variations quantified the contribution of each iteration to the deformation of each point of a star and thus gave a direct control of the convergence of our modelling. The point at which deformation was expected to be the most important is at the point closest to the Lagrangian point L1 in each star i: (R_{s, i}, μ = 0, ϕ = 0). Thus, if the variation of potential and density at this point was negligible, compared to the previous iteration, then the entire model had converged. In the rest of our work, R_{s, i} denotes the radial coordinate of the stellar surface in a given direction(μ, ϕ). Our normalised convergence criterion for a model at the iteration j was
or
Given the normalised precision of our 1D average input models of about 10^{−8} on the density, we used and as the criterion for convergence for each star. This value can be adjusted in order to reduce the computation time of our models if highprecision modelling is not required.
3. Posttreatment of the models
Once we obtained a final model with the desired precision, the models went into to a postprocessing phase. During this phase, the different properties of the models were extracted and used to potentially include our new 3D models in stellar evolutionary models using the method of Kippenhahn & Thomas (1970). Later revisited by Meynet & Maeder (1997), this method has been extensively used in the literature (Roxburgh 2004; Landin et al. 2009; Siess et al. 2013; Palate et al. 2013; Fabry et al. 2022) to couple deformations to 1D stellar evolutionary models.
Following this method, the stellar structural equations were averaged over isobars of our 3D models to give a one dimensional description of the stellardistorted structure usable by stellar evolution codes.
3.1. Mass conservation equation
Next, we introduce m_{p}, denoting the mass encompassed under an isobar of pressure p and the averaged radius r_{p}, such as , which is the volume under this isobar. With this notation, the mass conservation equation is unmodified and given by
By construction, the critical direction (μ_{crit}, ϕ_{crit}) was chosen to verify Eq. (16) at each layer of the star, implying that r_{p, i} = r_{i} and m_{p, i} = m_{i} in the whole star.
3.2. Hydrostatic equilibrium
Following Eq. (1), the infinitesimal distance dn between the equipotentials Ψ_{tot} and Ψ_{tot} + dΨ_{tot} is given by
With this formalism, the volume between two isobars is obtained with
where the integrals on Ψ_{tot} denote integrals taken over isobars of constant Ψ_{tot} and dσ is an isobar surface element expressed as
In Eq. (19), g_{eff} denotes the norm of the effective gravity presented in Eq. (1), and r(μ, ϕ) is the radial position of the chosen isobar in the direction (μ, ϕ). The analytical expression of the effective gravity in our new models is developed in Appendix B. Equation (18) can be integrated and combined to Eqs. (17) and (19) to obtain the volume under an isobar:
where Ψ_{tot, c} is the total potential at the centre of the star. Finally, the mass between two isobars is expressed as
Therefore, the equation of hydrostatic equilibrium can be expressed as
which can be written as
where,
3.3. Transport of energy
In this section, we assume that the chemical composition of each star is constant on each isobar. Therefore, the isobars become isotherms.
In the radiative zone, the norm of the local flux is given by
where T is the temperature of a given isobar, κ is its opacity, a is the radiation constant, and c is the speed of light in the vacuum. Equation (25) can be rewritten using Eq. (17), and the mass conservation equation (Eq. (16)) can be rewritten as
which is the classical formulation of the gravity darkening.
By integrating on an entire equipotential, we obtained the luminosity:
The temperature gradient could thus be expressed as
where we introduced the corrective factor
Finally, the correction of the radiative gradient is
To find the boundaries of convective zones, this corrective factor had to be applied to the radiative gradient. We note that the computed 3D deformations are mainly significant in the upper regions of stars, which means that f_{T}/f_{p} ≈ 1 close to the convective boundaries for most of the stars.
3.4. Conservation of energy
The energetic balance of a shell of mass m_{p} is trivial and gives the classical equation
where ϵ denotes the net energy production rate and s is the entropy. We note that to couple MoBiDICT to 1D stellar evolution codes, the computation of the effective gravity is required. The analytical expression of this effective gravity in our new models is given in Appendix B. Other main requirements included the detection of the equipotential and the interpolation of different quantities on each equipotential, assuming that the potential of each star in the direction (μ_{crit}, ϕ_{crit}) corresponds to the reference 1D stellar input model. In this work, f_{T} and f_{p} are also employed as structural deformation indicators in order to study the differences between models, as these quantities are isobar averages of the stellar deformation over the entire stellar structure.
4. Comparison with the Roche model
One classical way of describing the total potential in a binary system is the Roche model. While this model makes strong assumptions regarding stellar structure, it also allows for simple analytical expressions of the total potential at a given point of a system.
In this section, we explore the Roche model, looking in particular for the differences between the Roche model and our more accurate models. Section 4.1 is dedicated to a theoretical presentation of the Roche model and the different assumptions made by this model. In Sect. 4.2, we present the different stellar models and binary systems that we model and compare in later sections.
Section 4.3, we study the difference in surface deformation between the Roche model and our models for different types of stars. Section 4.4 is focussed on the differences of surface gravity that are essential to not only the computation of the gravity darkening but also the interpretation of binary system observations. Finally, in Sect. 4.5, we explore the differences in structural modelling of the Roche model and our modelling for different types of stars composing the binaries we selected.
4.1. Theoretical formulation of the Roche model
The Roche potential describes the total potential on a given point of the system that is produced by two pointlike gravitationally interacting bodies and the consequent centrifugal potential created. The main advantage of the Roche model is the possibility to have simple analytical expression of the total potential. To obtain such expressions, several assumptions are made: The orbits of the stars are assumed to be circular, and the stellar rotation is assumed to be synchronised to the orbital rotation. The stars are assumed to behave in a pointlike manner, which means that the mass of their envelope is neglected to compute the tidal potential. The stars are also assumed to have a spherically symmetric density profile, even in the most distorted cases.
The Roche model is overly used in the literature to compute the gravitational potential of distorted binaries at their surface and outside of stars (Siess et al. 2013; Palate et al. 2013; Fabry et al. 2022). In the most extreme cases, this potential is even used to compute structural deformation of the components of binary systems through the coefficients f_{p} and f_{T}, as presented Sects. 3.2 and 3.3.
In its classical formulation, the total Roche potential of a given point in the grid of the arbitrarily chosen primary of the system is given by
where Ψ_{1} and Ψ_{2} are, respectively, the gravitational potential of the primary and secondary stars, expressed as
depending on the author and the utilisation of the models. The same assumptions that were made in our code were made concerning the circular orbits and the synchronised rotation. Thus Ψ_{centri} is given by Eq. (13).
The gravitational potential inside a star given by the Roche model as expressed in Eq. (34) is not a solution of Poisson’s equation (this can be directly seen by inserting Eq. (34) in Eq. (2)). Therefore, using this formulation of the Roche potential to correct a stellar structure at the hydrostatic equilibrium, as done in Sect. 3, can be problematic. To solve this particular issue, we derived a “refined” version of the Roche model, guaranteeing that the gravitational potential is a solution of Poisson’s equation while keeping it untouched outside of the star and expressed as in Eq. (34). This version of the Roche potential, which we call the “refined Roche” in the rest of this article, is expressed as
In the following subsection, we focus on the differences in surface and structural deformations as well as the surface gravity obtained by the classical Roche model, our refined Roche model, and our nonperturbative modelling done with MoBiDICT.
4.2. Types of stars studied
In this work, we studied the surface deformation given by our models and the Roche model for several types of stars. We made four 1D stellar models with the Code Liegeois d’Evolution Stellaire (CLES; Scuflaire et al. 2008) in order to represent a diversity in the types of stars studied and compared in the work. The stellar properties of each model are presented in Table 1. Our method is not limited to these particular types of stars or binaries. As long as a 1D density profile of each star composing a system can be provided, the deformation can be computed. In the most extreme cases, our method is even able to model the deformations of systems with compact objects, such as white dwarfs.
Summary of the stellar properties of the 1D input models used in this work.
In this article, we focussed on twin systems with two identical stars in order to have a point of reference to compare the distortions of different stellar models, as the deformation of a binary system component depends on the stellar structure of its companion. However, MoBiDICT also works perfectly with two nonidentical components. In MoBiDICT, the type of star studied is only defined by the input density profile of each star, then associated to a total potential along the direction (μ_{crit}, ϕ_{crit}). We illustrated the different 1D density profiles normalised by the core density and the mass profiles in function of the normalised radius of each star studied in this article in Fig. 1.
Fig. 1. 1D input density profile normalized by the core density and mass profile of each component of the twin binary systems studied in this article. The upper panel corresponds to the stellar density profile while the bottom panel corresponds the mass profile normalized by the total mass of each star. 
Two groups of models can be distinguished in Fig. 1: models with high or low envelope mass fractions. While massive and solarlike stars have a nearly negligible envelope mass, for lowmass and red giant branch (RGB) stars, most of their mass is located in their envelope. In particular, the mass of about 80% of our RGB star is located in the convective envelope, with the transition from the radiative core to the convective envelope located at r/R_{1D} ∼ 0.05 for this star.
In the following sections, we compare the deformations obtained from the different types of stars presented above with the Roche models and with our new type of models. We then relate those differences directly to their 1D density and mass profiles presented in Fig. 1.
4.3. Surface deformation
The surface deformation of binary stars is directly related to its total potential. In our models, the surface of a star composing the binary system is defined by the equipotential corresponding to the photosphere of our 1D input stellar model in the direction (μ_{crit}, ϕ_{crit}), as presented in Sect. 2.3.
Both formulations of the classical Roche model and our refined version give the same deformations at the surface, as the total potentials have the same expressions outside of the star. In order to study the deformation given by the Roche model, we assumed, as done in our models, that the volumetric average radius coincides with the 1D input average radii, meaning that the volume and mass of the models are conserved. Using this property, we could define the surface given by the Roche model with the same criterion as our models. In Fig. 2, the position of the surface as given by the Roche model and our model is compared in the case of the 0.2 M_{⊙} star with an orbital separation of 2.8 R_{1D}. The R_{1D} denotes the radius of the 1D spherically symmetric model.
Fig. 2. Surface deformation of the twin binary system of our 0.2 M_{⊙} stars with a separation of 2.8 R_{1D} (corresponding to a period of 2 h and 11 min) as viewed from the side of the system. The black curve corresponds to the Roche lobe of the system, the yellow curve is the surface of each star given by the Roche model while the purple curve corresponds to the surface of each star given by the modelling with MoBiDICT. In this particular case, the filling of the Roche lobes is of 81.7% with each model. 
Figure 2 shows a visible difference in surface deformation between our models and the Roche model. To quantify this difference, we introduced a new quantity:
This quantity represents the difference in deformation at the closest point normalised by the deformation between the Roche model and the 1D input model. For the system presented in Fig. 2, this factor is ΔR = 0.368, which means that we have 36.8% more deformation with MoBiDICT than with the Roche model at the closest surface point to the Lagrangian point L1. The difference in surface deformation is a function of the orbital separation between the two components of the binary system. To draw a comparison between the different types of stars, we used the ratio of orbital separation to initial 1D stellar radius a/R_{1D} as an indicator of the orbital separation with respect to the type of stars in the system.
While we focus on the most distorted region in the rest of this deformation study, Fig. 2 also shows that surface deformation discrepancies appeared at the poles, where our model is more contracted than the Roche model. To study the significance of the modifications to the surface deformation obtained with MoBiDICT for the different types of stars, we computed a grid of deformed models of twin binary systems composed of each stellar type previously presented. For each system, we varied the orbital separation to explore the dependency of the deformation discrepancy to the orbital separation of the components. The results obtained with this grid of models are presented in Fig. 3.
Fig. 3. Difference in deformation at the closest point ΔR between Roche and our models as a function of the orbital separation of the components normalized by the input 1D stellar radius a/R_{1D}. 
Figure 3 shows that the differences of the deformations are extremely dependent on the type of star studied. In the case of a massive and solarlike MS star, Fig. 3 shows minor discrepancies of surface deformation, even when the stars are close to filling their Roche lobes. At most, the deformation discrepancy reaches 7% before a common envelope or mass transfer phase occurs.
As showed in Fig. 3, RGB stars are significantly more deformed by our modelling (up to ∼35%) than with the Roche model when a/R_{1D} ≲ 3. Moreover, this discrepancy remains significant even at large a/R_{1D} with an asymptotic value of ΔR = 0.10.
Finally, we observed that the corrections are the most significant for lowmass stars. As can be seen in Fig. 3, when a/R_{1D} ≲ 3.5, the Roche model totally underestimates the deformations up to ∼80% just before reaching a common envelope or mass transfer phase.
One would expect that for RGB stars, our new modelling does not significantly change the surface deformation obtained with the Roche model, as RGB stars evolve with a high density contrast between the radiative core and convective envelope, justifying the pointlike assumption of the Roche model. In practice with our modelling, we saw that by comparing Figs. 1 and 3, deformation discrepancies are directly related to the 1D mass profiles of each star. Stars with the largest mass fraction in their upper layers are the most impacted by our modelling. In particular, in our RGB model, the low density of the envelope is totally compensated by its significant volume, leading to an envelope mass of about 0.8 M_{tot} (see Fig. 1). MoBiDICT is necessary to model deformations of close binaries composed of stars with a high envelope mass fraction, such as RGB or lowmass stars.
All the differences previously presented can have a significant impact on the theoretical evolution of binary stars, as mass transfer occurs at higher orbital separation than predicted by the Roche model. Therefore, mass transfers and common envelope phases can occur in earlier system evolutionary stages than with the Roche model, impacting the future evolutionary path of such systems. In that respect, the significant deformation discrepancy found for RGB stars is crucial, as most of the mass transfers are expected to occur during late evolutionary stages, due to the rapid radius inflation resulting from stellar evolution. For example, in the system presented in Fig. 2, a common envelope phase is expected to occur when the orbital period reaches 2 h and 05 min with our modelling, while with the Roche model an orbital period of 2 h is required.
In Sect. 2.4, we mentioned and presented the different readjustments made at the end of each modelling iteration. One of these readjustments, presented in Appendix A, concerns modifying the computation of the system orbital period so that the nonspherical shapes of deformed binaries are taken into account. In the case of lowmass binaries, we saw an approximately 1% increase in the orbital period for a given orbital separation. Compared to the usual uncertainties in the determination of orbital periods, this effect dominates the uncertainties when linking the orbital period of a system to its orbital separation.
4.4. Effective surface gravity
On the observational side, more than having significantly larger deformations than the Roche model, with MoBiDICT the surface effective gravity is directly impacted, as the topology of the total potential is modified. The local brightness of a star as given, for example, in Eq. (26) is proportional to the surface local effective gravity. Therefore, a modification of the total potential topology can directly be seen in the light curves of eclipsing binary systems and impact the derivation of the orbital parameters, for example. As we focus on the surface effective gravity and both our refined Roche model and the classical Roche model give the same result, only one of them is considered in the rest of this section.
In our models, the effective gravity has an integral expression, given in Appendix B, and for the Roche model, the effective gravity is given by g_{eff}=∇Ψ_{tot, Roche}, where Ψ_{tot, Roche} is given in Sect. 4.1. For the Roche model and our models, Fig. 4 compares the surface effective gravity along the meridian passing by the most distorted point and the poles of the stars for a twin system composed of two 0.2 M_{⊙} stars with an orbital separation of 2.8 R_{1D}.
Fig. 4. Surface effective gravity of the meridian passing by the most distorted point of the star, its opposite, and the pole. The orange curve corresponds to the surface effective gravity given by the Roche model while the purple curve corresponds to the effective gravity obtained with our modelling. The system modelled here is a twin system composed of two 0.2 M_{⊙} stars with an orbital separation of 2.8 R_{1D}. The absence of points near θ = 0 is explained by the GaussLegendre quadrature method that is not creating any points in this region of our grids. θ > 0 corresponds the direction where the stars are facing each other, while the direction θ < 0 corresponds to the back of the system. 
Figure 4 shows that a difference of surface effective gravity exists between the Roche model and our modelling. The most significant g_{eff} discrepancies are located at the region closest to the L1 point, where we have a normalised the g_{eff} difference of 12% with respect to the Roche model.
The discrepancies found can be explained by several factors. First, as our models are noticeably more distorted than the Roche model, the surface effective gravity is consequently significantly lower at the surface. To a lesser extent, not correcting the Roche model for the redistribution of mass in the star can impact the surface g_{eff}.
To summarise the differences obtained, we computed the effective gravity of the grid of models presented in Sect. 4.3 and quantified the differences of surface effective gravity between our models with the quantity Δg_{eff, surface}, defined as
In Fig. 5, we summarise all the differences of surface effective gravity for the twin binary systems presented in Sect. 4.3 and for a wide diversity of orbital separations. Figure 5 shows the same overall results as Fig. 3. The models with most surface deformation discrepancies compared to the Roche model are the same models exhibiting the stronger surface effective gravity differences. In particular, Fig. 5 shows that higher discrepancies of g_{eff} are found for the RGB and lowmass stars when a/R_{1D} ≲ 3.2, as the other MS stars have limited corrections from our modelling.
Fig. 5. Difference of surface effective gravity of the most distorted point as a function of orbital separation for the four types of twin binaries studied. 
For the interpretation of observations, the modifications of g_{eff} can, in theory, strongly impact the light curves of close binary systems, as the local effective temperature of a star is directly proportional to the surface effective gravity. With our modelling, we expected that a star would appear colder than with the Roche model. In particular, the equator is the region that would appear noticeably colder and thus redder. To confirm these considerations, simulations of binary light curves using models from MoBiDICT have to be developed.
4.5. Structural deformation
The last point of comparison between the Roche model and our models is the structural deformation of the bodies composing the binary system. For this study, we compared the classical Roche model, our refined Roche model, and the model obtained through MoBiDICT. The structural discrepancy can be first seen in the total potential Ψ_{tot} and the effective gravity g_{eff}, as these quantities define the structural properties of our models. Figure 6 illustrates the evolution of Ψ_{tot} and g_{eff} in a particular direction, the axis linking the centre of both stars, and the system centre of mass for the 0.2 M_{⊙} stars with an orbital separation of a = 2.8R_{1D}.
Fig. 6. Normalized total potential Ψ_{tot} and effective gravity g_{eff} in the stellar interior along the axis joining the centre of the two stars. The system studied here is a twin binary system composed of 0.2 M_{⊙} stars with an orbital separation of a = 2.8R_{1D}. The blue and green curves are corresponding to the classical Roche model while the orange and violet curves are respectively our refined Roche model and our new types of models. 
Figure 6 illustrates the limitations of Roche models to study the stellar structure of deformed binaries. Classical ways of computing the Roche potential (Eqs. (33) and (34)) are inaccurate when modelling the stellar interior, while the surface properties are well reproduced. In addition, Fig. 6 shows that our refined Roche model is particularity accurate in obtaining both Ψ_{tot} and g_{eff}, even if discrepancies can be seen in the entire star.
In the work of Fabry et al. (2022), the potential used to compute the f_{p} and f_{T} factors, Eqs. (24) and (30), corresponds to the green curve in Fig. 6. However, with assumptions on its derivative, the corresponding effective gravity used is the expression given by our refined Roche model. Even if the effective gravity is correct and able to obtain some f_{p} and f_{T}, the potential used to identify the equipotentials on which the average quantities are computed is not representative of the stellar deformed structure. Corrections can therefore be expected on the averaged internal quantities f_{p} and f_{T} in Fabry et al. (2022), for example.
To quantify the structural differences on the entire stellar structure, we used the coefficients f_{p} and f_{T} presented in Sect. 3. As our work considers the entire deformed structure as a function of the radius, we limited our study to twin binary systems composed of our four types of stars presented in Sect. 4.2 with an orbital separation of 2.7 R_{1D}. The resulting f_{p} and f_{T} obtained with the different models are respectively illustrated in Figs. 7 and 8. As classical expressions of the Roche model lead to unphysical results on the majority of the stellar structure, f_{p} and f_{T} were not studied for these models.
Fig. 7. Evolution of f_{p} as a function of the radius of stars in twin binary systems for the different models compared. From top to bottom, each panel respectively represents f_{p} for the 0.2, 1.0, 1.5 and 20.0 M_{⊙} stars in twin binary systems. The orange curves are the results with our refined Roche model while the purple curves are the results from MoBiDICT. 
Fig. 8. Evolution of f_{T} as a function of the radius of stars in twin binary systems for the different models compared. From top to bottom, each panel respectively represents f_{T} for the 0.2, 1.0, 1.5 and 20.0 M_{⊙} stars in twin binary systems. The orange curves are the results with our refined Roche model while the purple curves are the results from MoBiDICT. 
Figures 7 and 8 show that models with the most deformation discrepancies with respect to the refined Roche model are also the models with the highest f_{p} and f_{T} differences. This structural difference can mainly be seen in the upper layers of the 0.2 and 1.5 M_{⊙} stars, as these bodies have the largest deformation discrepancies as compared to the Roche model. Moreover, in the case of the 0.2 M_{⊙} binary system, a discrepancy in the entire structure can be seen. By predicting a lower f_{p} and f_{T}, our modelling impacts the structure of binary stars. Combined with the discrepancies seen in deformations, our modelling can significantly alter the evolutionary path of binary systems. Finally, as previously seen, for MS massive and solarlike stars, our modelling does not significantly change their structural depiction.
5. Comparison with perturbative models
The principle of perturbative models is to treat the centrifugal and tidal forces as perturbations of the spherical symmetry by only accounting for the leading terms. The potentials and densities of the models are also projected on a basis of spherical harmonics by only taking the leading order, namely, the ℓ = 0 and ℓ = 2 terms. The objectives of this section are to determine the precision limits to only account for the leading order terms in the spherical projections and to see whether the treatment of deformations by a perturbative approach is justified in the most distorted cases.
In Sect. 5.1, we present the theoretical formulation of the perturbative models. Section 5.2 is dedicated to the exploration of the limitations of the selection of spherical harmonic projections of the gravitational potential, while in Sect. 5.3 we study the limit of the perturbative approach in the treatment of the tidal and centrifugal forces.
5.1. Theoretical formulation of perturbative model
As explained previously, the perturbative models are a class of models treating the centrifugal and tidal forces as a first order perturbation to the spherical symmetry. The gravitational potential is obtained, as in MoBiDICT, through Eq. (11) but only taking in account the first order terms: the ℓ = 0 and ℓ = 2. With this perturbative modelling, the feedback process where the modification of the potential leads to a redistribution of the masses (therefore impacting the potential that requires updating in an iterative process) is not present. Without going through a detailed demonstration of the equations used in the perturbative model, we mention that the assumptions made allowed us to simplify the expression of the deformations that are related to the structural coefficient η_{ℓ}. These coefficients are linked to the deformations of a star through their definition
where is the average radius of an isobar that is defined as
and that is different from r_{p} corresponding to the radius of the sphere with a volume identical to the one under an isobar. The radii of the equipotentials r(μ, ϕ) projected on the spherical harmonics basis are noted as and expressed as done in Eq. (9):
The overall deformation of a given point in a model can thus be obtained by projecting back these on the spherical coordinates basis as follows:
with p = 0 if ℓ is even, p = 1 otherwise, and m = 2k + p.
With the coefficient η_{ℓ} and using a simple integration, can be obtained, and thus the deformation of a model can be computed. Physically, these coefficient η_{ℓ} correspond to the structural answer of a body if a perturbative potential is applied. Using simplifications performed in the perturbative approach, η_{ℓ} become independent from the deformation undergone by a star and are only related to its unperturbed structure. In the perturbative approach, η_{ℓ} can be obtained by solving the wellknown ClairautRadau second order differential equation expressed as
using the boundary conditions
To compute the perturbative models, we developed a small feature in MoBiDICT to solve the ClairautRadau equation (Eq. (42)) with an order two RungKutta method and using the 1D input stellar models given in this tool.
For comparison with the model obtained with MoBiDICT, we decided to directly compare the η_{ℓ} of the two models. To obtain these coefficients with our models, we computed the deformations as described in Sect. 2. Then we obtained and projected the deformations on a basis of spherical harmonics (Eq. (40)) to have the that were used in Eq. (38) to compute the coefficients directly given by MoBiDICT. We mention that when using these equations, none of the assumptions of the perturbative approach were made.
5.2. Gravitational potential
One main approximation of perturbative models is to only account for the leading order terms (ℓ = 0 and ℓ = 2) on the projections of the potential and densities. Verifying these assumptions with MoBiDICT is straightforward, as our modelling directly gives the spherical expansion of the gravitational potential of each deformed star as an output of the code. Figure 9 illustrates the gravitational potential projected on the spherical harmonics basis of the twin binary system composed of the 0.2 M_{⊙} stars at an orbital distance of a = 2.8R_{1D}. This is the exact system configuration illustrated in Fig. 2.
Fig. 9. Normalized gravitational potential projected on the spherical harmonics of a star of 0.2 M_{⊙} in a twin binary system with an orbital separation of a = 2.8R_{1D}. The dotted lines are corresponding to the surface of each star of the system in the direction including the most distorted point of each star and the Lagragian point L1 (μ = 0, ϕ = 0). 
Figure 9 illustrates the amplitude hierarchy of gravitational potential spherical terms. As assumed by the perturbative model, the leading order terms are the ℓ = 0 and ℓ = 2 terms. The dipolar term (ℓ = 1) is the next most important term, especially outside of the star, as each spectral term is proportional to r^{−ℓ}. In these regions, the dipolar term is included in the computation of the tidal potential applied to the opposite star. Neglected by the perturbative model, the ℓ = 1 term is less than one order of magnitude smaller than the ℓ = 2 term, meaning that the tidal force is significantly greater with our models and the assumption to neglect the ℓ = 1 is not justified for the perturbative models in the most distorted cases. The next order terms (ℓ = 3, 4) have amplitudes similar to the ℓ = 1 term inside the star, while they are one order of magnitude smaller than this same dipolar term outside the star.
5.3. Validity of the ClairautRadau equation
In this section, we verify the limits of the perturbative approach for close binaries, in particular when the stars undergo high deformations. To compare the different models, as explained in Sect. 5.1, we directly looked at the coefficient η_{ℓ} obtained through Eq. (42) for perturbative models. With our models, these quantities are computed with a projection on the spherical harmonics of the deformations. More details are given in Sect. 5.1. We also compared these models to our refined Roche model through the same coefficients η_{ℓ} obtained in the same way as with our models. We focussed our work on the ℓ = 2 term, as this coefficient is assumed by the perturbative model to be the leading contribution to the apsidal motion of a binary system, and it is used, for example, in Rosu et al. (2020, 2022a, b) to constrain the stellar structure. The evolution of η_{2} as a function of the average radius is illustrated in Fig. 10 for the extensively studied system of this article, the twin binary system composed of 0.2 M_{⊙} stars at an orbital distance of 2.8R_{1D}.
Fig. 10. Evolution of η_{2} as a function of the average radius for the different models studied in this work. The system studied here is a twin binary system composed of two 0.2 M_{⊙} stars with an orbital distance of 2.8R_{1D}. 
Figure 10 illustrates the different significant results from our modelling. First, comparing the Roche model to other models, we observed that both perturbative and MoBiDICT models predict an important difference of η_{2} in all the stars. These results can be interpreted as difficulties of the Roche model to precisely depict the stellar deformations. These limitations originate from the assumptions on the stellar gravitational potential.
The second major result illustrated in Fig. 10 is the comparison between our new models and the perturbative models. First, in the less distorted regions of the star (when ), the same results were obtained with the models, indicating that our modelling gives a good depiction of stellar deformation in these regions, as we expected to be in the limit of validity of the perturbative modelling. Then, after , nonnegligible differences between the two models appeared, with discrepancies increasing with , as the upper layers of the star are more deformed. These differences can be explained by the limitations of the perturbative approach in these regimes of deformation, while our model is designed to study cases with high stellar deformations.
To study the apsidal motion of a binary system, one would require the apsidal motion constant k_{2} that is usually obtained through η_{2} evaluated at the surface of the primary star. To study how our new modelling procedure impacts this apsidal motion constant in regard to the twin binary system studied and orbital separation, we introduced a new quantity Δη_{2} defined as the difference of η_{2} at the stellar surface between our model and the perturbative model. The quantity Δη_{2} is expressed as
The comparison of Δη_{2} from the grid of the twin binary models presented in Sect. 4.2, as a function of their orbital separation, is illustrated in Fig. 11.
Fig. 11. Evolution of the surface Δη_{2} for different types of stars composing the binaries and with different orbital separation. 
Figure 11 illustrates that the differences of η_{2} between the perturbative and MoBiDICT models are almost independent from the type of star composing the binary system and only dependent on the orbital separation of the components. Significant differences of surface η_{2} can only be seen in the most distorted cases, when a < 5R_{1D}. Otherwise, the perturbative treatment of the tidal and centrifugal forces is justified to compute η_{2}. We expect the differences we found to directly impact the apsidal motion computation, resulting in corrections to apsidal motion constant k_{2} and adding k_{1} constant to account for neglected dipolar terms. The impact of our improved modelling on apsidal motion will be more deeply investigated in a forthcoming article.
6. Conclusion
In this work, we developed a new tool called MoBiDICT to compute 3D static models of close synchronised binaries in hydrostatic equilibrium in a nonperturbative way. Sections 2 and 3 were respectively dedicated to giving a technical description of this code and presented the posttreatment done in our models to couple, in the future, MoBiDICT to classical 1D evolution codes. We then compared our new type of modelling to classical models, namely, the Roche and perturbative models.
First, for the Roche model, our study showed in Sect. 4 that the differences of deformation mainly arise in very close binaries (when a < 5R_{1D}). The discrepancies are significant for lowmass and RGB stars in twin binary systems where the deformation differences can reach, at most, 80% for our 0.2 M_{⊙} stars. Massive and solarlike MS stars are less distorted, with a difference of 7% at most. The envelope mass is the key parameter controlling these deformation discrepancies. Stellar models with the largest envelope masses compared to the total mass were observed to be the most modified by MoBiDICT. All the quantities studied were impacted by this difference. In particular, significant changes of the surface effective gravity were seen, which directly affected gravity darkening and thus the interpretation of the observations of such stars. In the lower deformation regime, the Roche model provides a good approximation of the stellar surface properties, although we would not favour using Roche to study the structural properties of stars.
The discrepancies with respect to the perturbative approach were discussed in Sect. 5. Our work showed that the assumption of the perturbative model to only account for the leading terms of the projected gravitational potential (ℓ = 0 and ℓ = 2) is not justified in the highdistortion cases (when a < 5R_{1D}). In particular, the dipolar term (ℓ = 1) represented 31% of the quadrupolar term (ℓ = 2) outside of the studied star and 14% at its surface. In addition, we showed that in such regimes of distortion, the perturbative approach reaches its limit and significant deformation discrepancies can be seen, particularly in the upper layers of the stars (when r > 0.5R_{1D}). These differences of deformation were seen through the quantity η_{2} that was evaluated at the surface to estimate the apsidal motion of binary systems. In such a distortion regime, the theoretical apsidal motion from the quadrupolar term is generally considerably underestimated, and an additional term from the dipolar component is likely to arise. In the least distorted regime (a > 5R_{1D}), the perturbative approach remains valid, and accurate η_{2} values were obtained.
A few conclusions can be drawn from our work. First, in least distorted regimes, the classical modelling methods for binaries are valid, at the surface for the Roche model and in the entire structure for perturbative models. In the most distorted regimes, however, both the perturbative and Roche models fail to describe the structural and surface properties of stars with high envelope masses. In particular, significant deformation discrepancies were observed for lowmass stars, in agreement with the results of Landin et al. (2009). Our work highlights the necessity of nonperturbative modelling to study close binary interactions and evolution, as mass transfer is expected to occur earlier in the lifetime of a system under our modelling approach. In that regard, deformation discrepancies found for RGB stars significantly facilitate mass transfer occurrence during a stellar evolutionary phase already believed to be a source of mass transfer.
Our method, while effective in modelling highly distorted cases, has demonstrated limitations when deformations are extremely weak. Specifically, we observed numerical noise arising from higher spherical orders of the potential and densities that impact the derived η_{ℓ}. To address these issues, we plan to implement a coupling of MoBiDICT, in the least distorted regions, to a refined perturbative approach accounting for all the spherical orders and for terms of the tidal and centrifugal potential that are usually neglected. The advantage of this new form of the perturbative model is the possibility to study solidbody nonsynchronised systems that are also easily generalised in MoBiDICT. All these modifications will be the subject of a forthcoming article.
Acknowledgments
The authors are thanking the anonymous referee for his/her constructive comments. L.F was supported by the Fonds de la Recherche Scientifique FNRS as a Reaserch Fellow. The authors would like to thank Martin Farnir from “Psychedelic Kiss Designs” for his help in finding our code name and designing MoBiDICT’s logo.
References
 Aharonian, F., Akhperjanian, A. G., Aye, K. M., et al. 2005, A&A, 442, 1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fabry, M., Marchant, P., & Sana, H. 2022, A&A, 661, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fitzpatrick, R. 2012, An Introduction to Celestial Mechanics (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
 Gaia Collaboration 2022, VizieR Online Data Catalog: I/357 [Google Scholar]
 Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Han, Z., Podsiadlowski, P., Maxted, P. F. L., Marsh, T. R., & Ivanova, N. 2002, MNRAS, 336, 449 [Google Scholar]
 Han, Z., Podsiadlowski, P., Maxted, P. F. L., & Marsh, T. R. 2003, MNRAS, 341, 669 [NASA ADS] [CrossRef] [Google Scholar]
 Heber, U. 2009, ARA&A, 47, 211 [Google Scholar]
 Houdayer, P. S., & Reese, D. R. 2023, A&A, 675, A181 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [Google Scholar]
 Kippenhahn, R., & Thomas, H. C. 1970, in Stellar Rotation, IAU Colloq., 4, 20 [NASA ADS] [CrossRef] [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]
 Landin, N. R., Mendes, L. T. S., & Vaz, L. P. R. 2009, A&A, 494, 209 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Manchon, L. 2021, Ph.D. Thesis, Université ParisSaclay, France [Google Scholar]
 Meynet, G., & Maeder, A. 1997, A&A, 321, 465 [NASA ADS] [Google Scholar]
 Palate, M., Rauw, G., & Mahy, L. 2013, Cent. Eur. Astrophys. Bull., 37, 311 [Google Scholar]
 Patterson, J. 1984, ApJS, 54, 443 [NASA ADS] [CrossRef] [Google Scholar]
 PriceWhelan, A. M., Hogg, D. W., Rix, H.W., et al. 2018, AJ, 156, 18 [Google Scholar]
 Prša, A. 2018, Modeling and Analysis of Eclipsing Binary Stars; The Theory and Design Principles of PHOEBE (Bristol: IOP Publishing) [Google Scholar]
 Rosu, S., Noels, A., Dupret, M. A., et al. 2020, A&A, 642, A221 [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]
 Roxburgh, I. W. 2004, A&A, 428, 171 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Roxburgh, I. W. 2006, A&A, 454, 883 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sana, H., & Evans, C. J. 2011, IAU Symp., 272, 474 [Google Scholar]
 Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444 [Google Scholar]
 Sana, H., de Koter, A., de Mink, S. E., et al. 2013, A&A, 550, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sana, H., Le Bouquin, J. B., Lacour, S., et al. 2014, ApJS, 215, 15 [Google Scholar]
 Savonije, G. J. 1978, A&A, 62, 317 [NASA ADS] [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]
 Slawson, R. W., Prša, A., Welsh, W. F., et al. 2011, AJ, 142, 160 [Google Scholar]
Appendix A: Verification of Kepler’s third law
In this section, we show the system of equations solved at the end of each modelling iteration in order to readjust the centre of mass and the orbital period of our models. We consider an arbitrary synchronised binary system with an orbital separation a, an orbital rotation rate Ω_{⋆}, and a centre of mass noted CM. This system can be illustrated as in Fig. A.1.
Fig. A.1. Schema of an arbitrary chosen binary system composed of two stars with different masses. a_{1} and a_{2} are the orbital separation of respectively the primary and secondary star from the centre of mass. and are the unit vectors of the x_{1} and x_{2} direction considered in this problem. 
The orbital separation of the system a is a constant of the model. Thus, the centre of each star is at the dynamical equilibrium, meaning that the sum of the gravitational, tidal, and centrifugal forces is null in their centre. Mathematically this condition is expressed as
where x_{CM} is defined as a_{1}/a.
The system given in Eq. A.1 can be summed and rewritten as
The second equation of this system can be solved with MoBiDICT to obtain the necessary to maintain the orbital separation of the binary system. Having , the first equation can be solved to obtain the new position of the centre of mass of the system while accounting for the redistribution of mass in the stars.
Appendix B: Expression of the effective gravity
In this section, we develop the expression of the effective gravity used in MoBiDICT.
The effective gravity is a vector defined by
with each component of a given primary star noted as ‘star 1’ expressed as
In this case, all the necessary quantities are partial derivatives of solutions to differential equations that have integral representations. Consequently, numerical methods are not required to perform the derivatives of the potentials. The contribution to effective gravity from gravitational potential of the primary star is given by
The contribution from centrifugal force is expressed as
Finally, the gravitational contribution originating from the companion is written as
where J^{−1} is the inverse of the Jacobian matrix taken between the system of coordinates of our primary and secondary star, and expressed as
The norm of the effective gravity is simply given by
All Tables
All Figures
Fig. 1. 1D input density profile normalized by the core density and mass profile of each component of the twin binary systems studied in this article. The upper panel corresponds to the stellar density profile while the bottom panel corresponds the mass profile normalized by the total mass of each star. 

In the text 
Fig. 2. Surface deformation of the twin binary system of our 0.2 M_{⊙} stars with a separation of 2.8 R_{1D} (corresponding to a period of 2 h and 11 min) as viewed from the side of the system. The black curve corresponds to the Roche lobe of the system, the yellow curve is the surface of each star given by the Roche model while the purple curve corresponds to the surface of each star given by the modelling with MoBiDICT. In this particular case, the filling of the Roche lobes is of 81.7% with each model. 

In the text 
Fig. 3. Difference in deformation at the closest point ΔR between Roche and our models as a function of the orbital separation of the components normalized by the input 1D stellar radius a/R_{1D}. 

In the text 
Fig. 4. Surface effective gravity of the meridian passing by the most distorted point of the star, its opposite, and the pole. The orange curve corresponds to the surface effective gravity given by the Roche model while the purple curve corresponds to the effective gravity obtained with our modelling. The system modelled here is a twin system composed of two 0.2 M_{⊙} stars with an orbital separation of 2.8 R_{1D}. The absence of points near θ = 0 is explained by the GaussLegendre quadrature method that is not creating any points in this region of our grids. θ > 0 corresponds the direction where the stars are facing each other, while the direction θ < 0 corresponds to the back of the system. 

In the text 
Fig. 5. Difference of surface effective gravity of the most distorted point as a function of orbital separation for the four types of twin binaries studied. 

In the text 
Fig. 6. Normalized total potential Ψ_{tot} and effective gravity g_{eff} in the stellar interior along the axis joining the centre of the two stars. The system studied here is a twin binary system composed of 0.2 M_{⊙} stars with an orbital separation of a = 2.8R_{1D}. The blue and green curves are corresponding to the classical Roche model while the orange and violet curves are respectively our refined Roche model and our new types of models. 

In the text 
Fig. 7. Evolution of f_{p} as a function of the radius of stars in twin binary systems for the different models compared. From top to bottom, each panel respectively represents f_{p} for the 0.2, 1.0, 1.5 and 20.0 M_{⊙} stars in twin binary systems. The orange curves are the results with our refined Roche model while the purple curves are the results from MoBiDICT. 

In the text 
Fig. 8. Evolution of f_{T} as a function of the radius of stars in twin binary systems for the different models compared. From top to bottom, each panel respectively represents f_{T} for the 0.2, 1.0, 1.5 and 20.0 M_{⊙} stars in twin binary systems. The orange curves are the results with our refined Roche model while the purple curves are the results from MoBiDICT. 

In the text 
Fig. 9. Normalized gravitational potential projected on the spherical harmonics of a star of 0.2 M_{⊙} in a twin binary system with an orbital separation of a = 2.8R_{1D}. The dotted lines are corresponding to the surface of each star of the system in the direction including the most distorted point of each star and the Lagragian point L1 (μ = 0, ϕ = 0). 

In the text 
Fig. 10. Evolution of η_{2} as a function of the average radius for the different models studied in this work. The system studied here is a twin binary system composed of two 0.2 M_{⊙} stars with an orbital distance of 2.8R_{1D}. 

In the text 
Fig. 11. Evolution of the surface Δη_{2} for different types of stars composing the binaries and with different orbital separation. 

In the text 
Fig. A.1. Schema of an arbitrary chosen binary system composed of two stars with different masses. a_{1} and a_{2} are the orbital separation of respectively the primary and secondary star from the centre of mass. and are the unit vectors of the x_{1} and x_{2} direction considered in this problem. 

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.