© ESO 2022

1. Introduction

Massive stars (M_initial  ≳  8 M_⊙) give rise to a multitude of astrophysical phenomena on the stellar to the galactic scale. Binarity is responsible for several of these and it is expected that the vast majority of massive stars are in binary systems that interact during their lifetime (Vanbeveren et al. 1998; Sana et al. 2012), these effects cannot be ignored. Furthermore, population synthesis calculations predict that about a quarter of all O-type stars will merge with its companion (Sana et al. 2012). Stellar mergers are thought to be responsible for strong magnetism in massive stars (Ferrario et al. 2009; Wickramasinghe et al. 2014; Schneider et al. 2016, 2019) and are a potential source of the high rotation needed in models of long-duration gamma-ray bursts (Yoon & Langer 2005; Yoon et al. 2006; Aguilera-Dena et al. 2018) and superluminous supernovae (Justham et al. 2014; Aguilera-Dena et al. 2020). As mergers are preceded by an overcontact phase, understanding the properties of an overcontact model can give crucial insights into the merger process itself and the transients associated with them, or, in a broader sense, into the question of whether the merger happens at all.

In the past, a significant effort has been undertaken to observe and model the low-mass overcontact binaries known as W Ursae Majoris (W UMa) binaries. These systems have short periods of less than a day and total masses on the order of 1 M_⊙. From the observational side, surveys from Eggen (1967) and Binnendijk (1970) characterized the masses and mass ratio of the first identified W UMa systems and, most notably, they found that the majority of them have unequal mass components. Later, over a thousand low-mass, near-contact or overcontact binaries were found based on the light curve analysis of OGLE-I survey data (Szymanski et al. 2001). Extensive studies have aimed to characterize the population of W UMa binaries and identify systems that are expected to merge (Gazeas et al. 2021). Recently, more and more massive overcontact binaries have been identified, the most massive of which, VFTS 352 (Almeida et al. 2015), has components over 30 M_⊙.

The fact that a large number of unequal overcontact binaries were observed challenged the argument from Kuiper (1941), who stated that overcontact stars cannot be stable; however, it was shown by Lucy (1968) that for the convective envelopes of W UMa binaries, this argument does not apply and they constructed an approximate model of a low-mass overcontact binary. Shu et al. (1976) expanded on these models and argued there must be significant energy transport in the innermost overcontact layers between both components. However, the theory of overcontact models of massive stars with radiative envelopes has not yet been well developed.

In the context of gravitational wave astrophysics, massive overcontact binaries have been proposed as progenitors of merging binary black holes (BBH) (Marchant et al. 2016; Mandel & de Mink 2016), as tidal synchronization leads to rapid rotation and enhanced rotational mixing (de Mink et al. 2009). Recently, new analysis techniques of observations have allowed for an improved characterization of the components of a massive overcontact binary (Palate et al. 2013; Abdul-Masih et al. 2020, 2021). Such observations give critical constraints to evolutionary models that try to explain the evolution of overcontact systems on their way to BBH mergers. However, to model these strongly interactive systems, important approximations need to be made (Marchant et al. 2016; Menon et al. 2021). Even though orbital evolution, tidal torques and mass transfer by Roche lobe (RL) overflow is usually taken care of, whereas the tidal deformation of the stellar structure from a binary companion or energy transfer in overcontact layers are ignored.

Since centrifugal forces introduce a latitudinal dependence and tidal forces an additional azimuthal one, in order to properly compute stellar models for single rotating or binary stars, the stellar structure equations would need to be recast and solved in 2D or 3D, respectively. From a computational standpoint, doing full 3D stellar evolution is not feasible with current resources, so simplifying assumptions need to be made in order to remain in 1D. Kippenhahn & Thomas (1970) introduced a methodology to correct the spherical stellar structure equations to account for deformation from the spherical geometry by rotation, which was later extended by Endal & Sofia (1976). This method is built on the premise that hydrostatic layers of stars fall on equipotential surfaces of the system. For single star rotation, this method has been implemented in the stellar evolution code Methods for Experiments in Stellar Astrophysics (MESA, Paxton et al. 2011, 2013, 2015, 2018, 2019). Currently, these single star corrections are also called upon when using the code to evolve binary stars. This means that while the centrifugal deformation from stellar rotation is accounted for, tidal deformation as a result of a nearby companion is not.

In this first paper of a three-part series, we apply the method of Kippenhahn & Thomas (1970) and Endal & Sofia (1976) to synchronized binary stars, as a first step in consistently modeling tidally deformed stars. The problem of energy transfer in overcontact layers will be studied in the next paper of this series. Still, the methods developed here are applicable to model twin overcontact binaries (binaries of mass ratio unity), where no energy is transferred between components, as well as for very close detached and semi-detached systems that are tidally synchronized. In Sect. 2, we calculate the modified 1D stellar structure equations in a conservative potential and compute the structure correction factors f_P and f_T for a Roche binary. In Sects. 3 and 4, we consider the Eddington limit and boundary conditions as appropriate to the modified equations, respectively. Section 5 shows stellar models computed using the MESA code comparing currently implemented physics with the new methods in representative case studies. Lastly, Sect. 6 presents our final remarks. All the input files needed to reproduce the results and associated data products in this paper can be downloaded from the Zenodo data repository¹.

2. Methods

2.1. Modifications to the spherical stellar structure equations

In non-rotating stars, the equations of stellar structure and evolution are solved in a 1D space of the radius, r, or mass coordinate, M, and assume spherically symmetric stellar shells. As such, the equation of mass continuity takes the following form:

(1)

with ρ as the density at radius, r, and M = M(r) as the total mass interior to r. The equation of hydrostatic equilibrium is written as follows:

(2)

where P is the total pressure and G the gravitational constant. Finally, in radiative regions of a star, the energy transport is described by:

(3)

with L as the luminosity, T the gas temperature, κ the local opacity, a the radiation constant, and c the speed of light. However, several physical effects can break the spherical symmetry that is assumed when calculating 1D stellar evolution models. Centrifugal forces introduce a latitudinal dependence and a binary companion introduces an additional azimuthal dependence. Yet, if the total potential Ψ, which is the sum of all gravitational and rotational contributions, is conservative, Kippenhahn & Thomas (1970) showed that such effects can be accounted for in 1D equations. The only assumption that remains is the one of shellularity: We assume that stellar isobars exactly coincide with equipotential surfaces of the given geometry, therefore, the total pressure, P, is constant at constant Ψ. von Zeipel (1924) showed that in hydrostatic equilibrium, a solidly rotating stellar shell has all its thermodynamical quantities, including T and ρ, remaining constant along isobars. However, he also showed that in this case, the stellar shell cannot be in radiative thermal equilibrium. Therefore one expects a deviation from shellularity coupled with large scale meridionial flows (Eddington 1925; Sweet 1950). Still, large horizontal flows are also present, keeping the rotation of isobars near to solid (Zahn 1992), so that the star is kept approximately shellular.

Instead of spherically symmetric stellar shells, we write the equations of stellar structure in terms of these equipotentials surfaces, which can be characterized by a single independent variable and thus be appropriate for 1D stellar evolution. Here, we summarize the model of Endal & Sofia (1976, henceforth ES) that cast the equations of stellar structure in 1D given a conservative potential Ψ. This potential can, for example, be the Roche potential of a single rotating star (Sect. 2.2) or of a binary star (Sect. 2.3).

We start by defining r_Ψ as the equivalent radius of a sphere enclosing the same volume V_Ψ enclosed by an equipotential surface:

(4)

Hence, as the spherical relation from radius to volume is retained and the mass density, ρ, is assumed constant over an equipotential, the equation of continuity is unaffected:

(5)

where M_Ψ is the total mass enclosed by the equipotential. The variable r_Ψ (or M_Ψ in a Lagrangian prescription) goes on to act as the new independent variable to describe stellar models.

Next, we denote with g the magnitude of the effective gravity at each point in space:

(6)

with dn being the distance between potential surfaces Ψ and Ψ + dΨ. Considering then a change in equipotential volume, we have:

(7)

(8)

where S_Ψ is the surface area of the equipotential, and we defined for convenience, the average of a quantity f over an equipotential as:

(9)

Continuing from Eq. (8), we have

(10)

which, along with the statement of hydrostatic equilibrium ∇P = −ρ∇Ψ gives the modified momentum equation:

(11)

with the correction to the standard spherical equation contained in the factor

(12)

This factor thus measures the strength of the pressure gradient acting on a deformed shell (as a result of the non-spherical potential) with respect to a spherical shell enclosing the same volume and mass.

Similarly to the mass continuity equation, the equation of energy conservation retains its spherical form as mass, radius, and volume relate to each other spherically, hence,

(13)

where L_Ψ is the energy flowing through the equipotential per unit time, ϵ is the (specific) nuclear energy generation rate, U the specific internal energy, and t is time.

Considering lastly the equation of energy transport through radiative layers of the star, we have for the radiative flux:

(14)

by the chain rule. Equation (14) is to the von Zeipel theorem (von Zeipel 1924), stating that the radiative flux ℱ is proportional to the local effective gravity g = |∇Ψ|. Taking then the integral over an equipotential surface:

(15)

using the previous results, namely, Eqs. (6), (9), and (10). Rewriting this using (11) and (12) transforms this into the usual spherical-like form of radiative energy transport:

(16)

with the second correction factor to the spherical equations being

(17)

Similarly to f_P (Eq. (12)), the correction f_T/f_P measures the strength of the radiative gradient of a deformed stellar shell relative to that of a spherical shell of the same volume and enclosed mass.

2.2. Single star rotation

A stellar shell rotating with angular velocity Ω exterior to a mass M has a total potential defined by

(18)

where θ is the polar angle and Φ is the gravitational potential due to the mass interior to the shell. Since stars, especially massive ones with radiative envelopes, are heavily condensed toward their centers, significant deformation from spherical symmetry is only expected in the outermost layers, justifying that the gravitational potential Φ of the star is taken as of a point mass, with M as the total mass of the star.

Paxton et al. (2019) calculated analytical approximations to the corrections, f_P and f_T, and provided polynomial fits for these as a function of the variable , with r_e the equatorial radius of the shell. In Appendix A, we provide new, more accurate fits without increasing the degree of the polynomial expansions.

2.3. Synchronized Roche binaries

We expect that tidal torques in binary systems are very sensitive to the separation of components (Zahn 1975), causing stars that significantly fill their RL to be synchronized with their orbits. Therefore, the total potential of these systems is the well-known Roche potential and may be written in the corotating frame as follows:

(19)

where Φ is the gravitational contribution due to the two component stars, a their separation, while the second term represents the centrifugal contribution of the synchronous rotation around the center of mass. We employ the same reasoning as in Sect. 2.2 to approximate the gravitational potential Φ by point mass potentials of mass M₁ and M₂, that is, the total masses of the component stars. Finally, defining the mass-ratio q = M₂/M₁, normalizing all distances r → r/a and shifting , we arrive at the dimensionless form of the potential in spherical coordinates:

(20)

where r² = x² + y² + z², r sin θ cos ϕ = x, r sin θ sin ϕ = y and r cos θ = z relate to the Cartesian frame, and . In Fig. 1, we show an illustration of the Roche geometry in the equatorial plane and highlight the colinear Lagrangian point equipotentials.

Fig. 1. Illustration of the Roche geometry in the equatorial plane. It shows the L1, L2, and L3 equipotentials (gray), as well as the splitting surfaces (blue dashed). The splitting surfaces have their normals always perpendicular to the potential gradient.

2.4. Overcontact shells

In the context of 1D evolution models, we have to distinguish the two stellar components within the geometry of a Roche binary. For equipotential shells lying within their respective RLs, this is trivial, as such layers are physically separated, but shells in overcontact are shared between the two stars. We therefore construct three “splitting surfaces”, one through each of the colinear Lagrangian points L₁, L₂, and L₃, separating the Roche geometry in three main parts. This will ensure distinction between the two components, as well as either from regions beyond the outer Lagrangian points L₂ and L₃, from which outflows are expected if any component overflows them. We construct them in such a way that the condition:

(21)

holds on all points of these splitting surfaces (see Fig. 1, in blue dashed lines), such that the surface normal is always perpendicular to the gradient of the potential (20). We note that by our construction of the splitting surfaces, the scalar product dn ⋅ dS is zero. This ensures that Eq. (7) still holds, and if the overcontact layers are shellular, there is no net flux of radiative energy between the components following Eq. (14) However, since the condition of shellularity cannot be held exactly, there will exist a flow L_Ψ, trans = ∫_dΨ[BOLD]ℱ ⋅ dS_split across equipotentials, representing energy transfer from one component to the other. In this work, we assume this contribution is zero and defer the modeling of energy transfer to later work. We note also that in the case q ≠ 1, the L₁ splitting surface is not a vertical plane through L₁, which differentiates our results of equipotential volumes and surface areas from those of Mochnacki (1984) and Marchant et al. (2021).

2.5. Numerical calculations

In order to use the modified stellar structure equations in a stellar evolution code, we require for each shell in the geometry the quantities V_Ψ, S_Ψ, ⟨g⟩ and ⟨g⁻¹⟩ which specify r_Ψ, f_P and f_T of Eqs. (4), (12), and (17). Additionally, the specific moment of inertia, i_rot, is needed to calculate the rotational velocity of the equipotential shell from the angular momentum: Ω_Ψ = j_rot/i_rot. We approach this problem using numerical methods. First of all we note that the problem is symmetric under y → −y and z → −z, so we need only explicitly do calculations in the quadrant of positive y and z.

Given a set of mass ratios, q, and target equipotential values, Ψ₀, we calculate, in a spherical coordinate system on rays of constant (θ, ϕ), the points r = r(θ, ϕ) that lie on the corresponding equipotential surface, such that:

(22)

When solving Eq. (22), care has to be taken in order to avoid roots beyond the splitting surfaces (Sect. 2.4), as these belong to a part the geometry not associated to the considered component star, and would instead belong to the other component or to the region beyond the outer Lagrangian points. We therefore restrict the solver to values r smaller than the value r_split(θ, ϕ) on the splitting surface that the ray (θ, ϕ) crosses. If it finds no roots in the interval (0, r_split), r_split is recorded instead as the splitting surfaces mark the maximal extent of the components.

Next, for these points r(θ, ϕ), we numerically integrate the following quantities:

(23)

(24)

(25)

(26)

(27)

where the and we used the notation dm_Ψ = ρdndS to signify a small mass element within the equipotential shell, such that dM_Ψ = ∫_Ψdm_Ψ. From these, f_P and f_T can be computed for each equipotential shell with Eqs. (12) and (17).

2.6. Integration results

We consider the component of a Roche binary at the origin with mass M₁ and a companion with mass M₂, and we computed all the necessary integrals of Eqs. (23)–(27) in the two-dimensional parameter space of mass ratio q ≡ M₂/M₁ and Roche equipotentials, Ψ₀. We emphasize here that we calculate solely the properties of the component with mass M₁ and that the properties of the companion M₂ can be obtained by inverting the mass ratio. We sampled 280 equally spaced mass ratios in logarithmic space in the interval log q ∈ [ − 7, 7], as well as the equal mass ratio log q = 0 case while in potential space, we sampled 158 values in several logarithmic ranges from Ψ₀ = 50Ψ_L1 to of different densities. With L_out we denote the outer Lagrangian point of the considered component of the binary (always having x_Lout < 0, which is L₃ if q ≤ 1 and L₂ if q ≥ 1). As quantities exhibit important variations around the overflow potentials Ψ_L1 and Ψ_Lout, those will be the regions we sample more densely. We note that Ψ_L1 < Ψ_L2 ≤ Ψ_L3 < Ψ_L4 = Ψ_L5 as in our convention Ψ is strictly negative.

Figure 2 shows the integration results of f_P for various mass ratios as a function of their radius relative to the RL radius r_Ψ/r_RL, where we defined r_ΨL1 ≡ r_RL. We observe discontinuities at the Lagrangian overflow points, which is a result of the splitting surfaces limiting the considered volume to one component (Fig. 1). For a comparison with the single rotating star deformation model, we overplot the equivalent f_P values as if the binary star was interpreted as a single rotating star with the same volume and rotational velocity (see Appendix B for details).

Fig. 2. Structure correction factor fP as a function of fractional RL radius rΨ/rRL for various mass ratios. The discontinuities in all graphs correspond to the crossing of the Lagrangian points L1 (at rΨ/rRL = 1, by construction) and Lout at larger radii. For reference, the equivalent single rotating star inferred fP is overplotted in dashed lines.

2.7. Comparison to literature

As verification of our integrations, we compare our results to calculations performed by Mochnacki (1984). As this study considered vertical splitting surfaces through L₁ at all mass ratio’s, for overcontact layers we are limited to comparing the mass ratio unity case. Figure 3 shows excellent agreement between our f_P values and those computed via Eq. (12) using the results from Mochnacki (1984), with the relative difference being less than about 1 part in 10³. A resolution convergence test of our integrations show our calculations converge to the 10⁻⁵ level at eight times our original default resolution. The final computations of all integrals for use in stellar evolution instruments are done with this eight-fold resolution increase, and corresponds to dividing the interval ϕ = [0, π] in n_ϕ = 5280 parts and cos θ = [0, 1] in n_θ = n_ϕ/2 parts so that Δϕ = π/n_ϕ and Δ cos θ = 2/n_ϕ.

Fig. 3. Comparison of the calculations of the structure correction factor fP from Mochnacki (1984) and those in this work, as a function of the fractional RL radius rΨ/rRL. We performed a resolution convergence test with double, quadruple and octuple the default spacial resolution, and conclude our integrations have well converged at an eight-fold resolution increase.

3. The ΩΓ limit

Thanks to the radiation pressure of the photon flux generated in the core, a star gains radiative support against its self-gravity. In the situation where these effects are canceled out exactly, the star has reached the Eddington limit, which is given as (for a single, non-rotating and thus spherical star):

(28)

Defining then the Eddington factor Γ as:

(29)

with L_rad the radiative luminosity, equal to L in a radiative envelope, Γ = 1 is then an equivalent statement of the Eddington limit. If Γ > 1, the sum of radiative and gravitational acceleration is directed outward. In the interior of the star this can in principle be compensated by an inversion of the gas pressure (Joss et al. 1973), but if a star approaches Γ = 1 at the surface, strong outflows are expected to develop (Gräfener et al. 2011).

After adding a centrifugal contribution to the force balance, we have:

(30)

where we write g = ∇Ψ as the effective gravity. Langer (1997) then considered the Eddington limit on the equator of a rotating (but non-deformed) star, and, assuming that the radiation field is isotropic, the radiative flux is ℱ = L/4πR², so that the force balance reduces to

(31)

where we have defined ω = Ω/Ω_{crit, class} as the fractional classical critical rotation rate . It is important to make the distinction here as the formal critical velocity is lowered to:

(32)

Therefore, the closer a star is to the Eddington limit, the lower the critical velocity will be.

Maeder & Meynet (2000) later noted that this treatment does not take von Zeipel’s theorem into account, namely that the radiative flux is dependent on local effective gravity, Eq. (14). Expressing this for a rotating star to lowest order in rotational velocity Ω gives

(33)

where ϑ signifies a polar dependence and Γ_Ω ≡ Γ × f(Ω, ϑ) is now a rotation-dependent Eddington factor. Since the strength of the radiative flux is a function of the effective gravity through the von Zeipel effect, Maeder & Meynet (2000) showed that, neglecting the ϑ dependence, Eq. (33) bifurcates when Γ = 0.639. Below this, g = 0 is the only solution, so that, independent of luminosity, the break up velocity is the classical velocity , with r_e the equatorial radius of the rotating star. Above the bifurcation point, both g = 0 and, importantly,

(34)

determine the critical velocities, the second of which is lower than the classical critical velocity. The lowest of these two velocities then determines the physical break-up velocity of a rotating star.

We further expand on the notion of break-up limits by considering the case of deformed stars due to a conservative potential. As a starting point, we consider the same criterion of zero net acceleration, g_tot = 0, as the stability boundary. The contributions to this acceleration are the effective gravity, g, as a result of the potential, and the radiative acceleration. Using von Zeipel’s expression for radiative flux (Eq. (14)), and combining with the equation of hydrostatic equilibrium results in:

(35)

which, upon using Eq. (16) as the radiative gradient in the non-spherical geometry, gives the stability criterion:

(36)

We recognize however that von Zeipel’s law of radiation can break down close to critical velocity due to significant baroclinicity, that is, departure from shellularity (see for example Espinosa Lara & Rieutord 2011). Continuing on our assumption of shellularity, Eq. (36) is of a very similar form as that of Maeder & Meynet (2000), namely, g(1 − Γ_Ψ) = 0, where Γ_Ψ = Γ × f is again a product of the classical Eddington factor with a function dependent on the geometry. For the single rotating star deformation, f = f(ω), dependent on the fractional critical rotation rate, while in the tidal deformation case of f = f(r_ψ/r_RL, q), dependent on the mass ratio and degree of Roche filling. In Eq. (36), the classical solution g = 0 is present of course, where the effective gravity vanishes as a result of rotational support only. As above, in this case there is no corresponding Eddington luminosity as it is independent of any radiative acceleration due to the Von Zeipel effect. The other solution is given by:

(37)

which translates to a modified expression of the Eddington luminosity:

(38)

or the Eddington factor (see also Sanyal et al. 2015):

(39)

These expressions, as in the case of the equations derived in Sect. 2.1, are applicable in any conservative potential Ψ.

As the criteria for instability, Eqs. (31), (34), and (37) all express a departure of the Eddington limit of a star from the classical Eddington limit of Γ = 1. In the model developed here, it amounts to a reduction of the maximal Eddington factor by the ratio f_P/f_T, which is smaller than one in both the single rotating star or synchronized binary case. In Fig. 4, we show the comparison of this reduction from the results of Langer (1997) to this work in the case of a single rotating star. Moreover, consistent to Maeder & Meynet (2000), we find a bifurcation in the critical rotation velocity at at Γ = 0.639. Below this number, only the classical critical velocity ω = 1 will make the star unstable, while above it, this velocity is reduced. The case of a synchronized Roche binary is shown in Fig. 5. For example, an equal mass ratio binary overflowing to its outer Lagrangian point has its Eddington limit reduced to about 84% its classical value.

Fig. 4. Maximal Eddington factor Γmax that a rotating star can have as a function of fractional classical critical rotation rate of ω = Ω/Ωcrit, class. Our model of shellular stars has a maximal reduction to Γmax = 0.639.

Fig. 5. Reduction, fP/fT, of the maximal Eddington factor as a function of mass ratio, q, and fractional overflow radius, rΨ/rRL. The equivalent radius of a star filling up to its outer Lagrangian point in marked with the green line.

4. Atmospheric boundary conditions

The equations of stellar structure must be supplied with appropriate boundary conditions (BCs) at the center and atmosphere of the star. In the center, the conditions are such that r = M = L = 0, which transform to r_Ψ = M_Ψ = L_Ψ = 0 in our equipotential shell model. At the surface, an atmosphere model is integrated to supply the pressure P_surf and temperature T_surf. A common choice is the so called gray Eddington atmosphere, which is a plane parallel model integrated from τ = 0 to some predetermined atmospheric optical depth τ_surf, and is presented in, for example, Cox & Giuli (1968, henceforth CG).

In spherical stars, the surface optical depth τ_surf is reached at the same physical perpendicular depth, z, at all points in the star. For deformed shellular stars, however, the optical depth τ across an equipotential depends on local effective gravity, as follows:

(40)

with κ and ρ as the opacity and density of the equipotential, respectively. Consequently, equipotential surfaces do not coincide with surfaces of constant optical depth. For example, in single rotating stars, a fixed τ_surf is reached at a deeper equipotential on the poles compared to on the equator, resulting in a different observed effective temperature between those regions. The same is true for highly deformed binary stars like overcontact systems. The atmosphere properties output by stellar evolution codes should therefore also take into account this effect of position dependent atmospheric depth.

Within the treatment of deformed stars in a conservative potential, we proceeded as follows. We start by defining the global effective temperature T_eff, Ψ of an equipotential by the law of Stefan-Boltzmann:

(41)

that is, the global effective temperature of an equipotential shell is that of a black body radiating a luminosity L_Ψ over an area S_Ψ, which need not be spherical. Next, using the equation of radiative flux Eq. (14), we can define a local effective temperature T_eff, l that varies across the equipotential:

(42)

We now consider a gray, plane parallel, Eddington-approximated atmosphere to construct appropriate BCs. At each point on the surface equipotential, the temperature profile can be written as (CG):

(43)

This expression is derived using two standard assumptions that can be similarly applied to the case of deformed stars. First, it is assumed that the radiation pressure at all optical depths can be written as . Second, the intensity of radiation on top of the atmosphere is assumed to be isotropic in the outward directions, meaning that we ignore limb darkening, so that the pressure at the τ = 0 surface is computed to be (CG):

(44)

The first BC we consider constrains the temperature T_surf of the outermost boundary by requiring that T_surf = T_eff, Ψ. To find the second BC on the surface pressure P_surf, we consider the equation of hydrostatic equilibrium:

(45)

Taking a point on the surface equipotential where T_eff, l = T_eff, Ψ, then by Eq. (42), we have that g = ⟨g⟩. As we also require that T_surf = T_eff, Ψ, Eq. (43) implies that the surface optical depth is τ_surf = 2/3. The surface pressure can then be integrated at this point of the equipotential surface by assuming a thin atmosphere such that g and κ are constants.

(46)

where we substituted ⟨g⟩ with the definitions of f_P and f_T in Eqs. (12) and (17), respectively.

In summary, the atmospheric BCs are represented by the following system of equations:

Tsurf = [4]LS We emphasize that this surface temperature T_surf of the stellar models is not the local effective temperature across the whole surface equipotential, the variations of which can be observed. However, it is possible to relate the surface temperature output by stellar models using these BCs to an observed local effective temperature at any point on the surface equipotential by applying Eq. (42). Finally, we mention that going forward, whenever the effective temperature is referenced, we mean the global effective temperature as defined by Eq. (41).

5. Stellar model comparison

We assess the effect of including tidal distortion on the structure of the star as well as the modified boundary conditions by computing stellar evolution models using the 1D MESA code, version 15140, with the mesasdk-x86_64-macos-21.2.1 SDK. Starting at the zero age main sequence (ZAMS), we compare the evolution of stars in three cases: (1) a 5 M_⊙ single rotating star, (2) a 5 M_⊙ detached star in orbit with a compact object, and (3) a 32 M_⊙ twin binary in overcontact.

The next subsection lists the physical assumptions made when computing stellar models. The following subsections detail and discuss the results of the respective case studies.

5.1. Physical ingredients

5.1.1. Mixing, microphysics, and stellar winds

In all stellar models, we used the Ledoux criterion (Ledoux 1947) for determining convective regions. Within such regions, convection itself is modeled via the mixing length theory of Böhm-Vitense (1958) as described by CG, with the mixing length parameter α = 2. We used efficient semiconvective mixing (Langer et al. 1983), with α_sc = 100 as proposed by Schootemeijer et al. (2019). Thermohaline mixing is applied as prescribed in Kippenhahn et al. (1980), with an efficiency of α_th = 1. During the main sequence, the convective core is allowed to overshoot into the radiative envelope. We use a step overshoot as calibrated by Brott et al. (2011) with f = 0.345, f₀ = 0.01, meaning that the diffusion coefficient is taken as constant from f₀ pressure scale heights into the convective boundary up to f − f₀ pressure scale heights out of the boundary. All other extra mixing processes, including rotational mixing from, for example, Eddington-Sweet type circulations are ignored in these sample models so as to isolate the effects of including the different forms of rotational deformation.

In terms of microphysics, we use the basic nuclear network consisting of ¹H, ³He, ⁴He, ¹²C, ¹⁴N, ¹⁶O, ²⁰Ne and ²⁴Mg, which is sufficient to follow the main sequence evolution (Paxton et al. 2011) and the reaction rates are taken from the JINA library (Cyburt et al. 2010). MESA uses a blend of different equation of states (EOS) from Saumon et al. (1995), Timmes & Swesty (2000), Rogers & Nayfonov (2002) and Potekhin & Chabrier (2010), as described in Paxton et al. (2019). Radiative opacities are similarly determined from a blend of opacity tables from Iglesias & Rogers (1996) and Ferguson et al. (2005). The solar metallicity Z_⊙ = 0.0142 as well as the relative metal fractions from are taken from Asplund et al. (2009).

Wind mass loss of individual stars is included following the prescription of Brott et al. (2011). As we focus on main sequence evolution, where the surface hydrogen fraction is X > 0.7, this method takes rates from Vink et al. (2001) when the effective temperature is higher than the iron bi-stability jump (as calibrated by Vink et al. 2001). For effective temperatures below the bi-stability, the mass loss rate is computed as the maximum of the rates of Vink et al. (2001) and Nieuwenhuijzen & de Jager (1995), where the latter is scaled with the same metallicity factor (Z/Z_⊙)^0.85 as with the rates of Vink et al. (2001).

5.1.2. Rotation and binary physics

At all times during the evolution of all models, angular momentum is continuously diffused throughout the whole star so as to enforce rigid rotation during the evolution. This, as with the omission of rotational mixing, is a way to reduce differences between the different rotational distortion models. When modeling single star rotation, the structure correction factors f_P and f_T that we use are the ones presented in Appendix A. In binary models, the angular velocity is always synchronized to the orbital period, so that the Roche potential (20) describes the geometry of the system, and the timescale of this synchronization is the orbital period itself. The structure corrections f_P and f_T are determined by interpolating the results of Sect. 2.6 on a grid of mass ratio and fractional RL radius (q = M₂/M₁, r_Ψ/r_RL), where the RL sizes are determined according to the Eggleton approximation (Eggleton 1983)². We do not model mass nor energy transfer in these models, as this is not needed for detached or equal mass binaries.

5.2. Main sequence evolution of a single rotating star

We compare the main sequence evolution of a 5 M_⊙ star rotating initially at around 70% critical velocity at Solar metallicity with and without the modified atmospheric BCs developed in Sect. 4. In MESA, the default BCs are integrated from a standard spherically symmetric, plane parallel, gray Eddington atmosphere (Eddington 1926). Our new BCs, however, take into account the modified pressure and area of the outer surface of the star, as represented by the correction factors, f_P, f_T, and the non-spherical area, S_Ψ, in Eqs. (47) and (48).

Figure 6 shows the evolutionary track of the single rotating stars in a Hertzsprung-Russell diagram. We see that the change in BC has little impact on the overall evolution, with the largest difference of the surface properties occurring near hydrogen depletion (X_c = 0.001), where at constant L, the corresponding T_eff varies by about 0.7%. Since the spherical BC calculates T_eff from , in comparing with Eq. (48) at constant luminosity and radius results in the temperatures differing by a factor . We plot the quantity as a function of stellar age in Fig. 7. As the star loses thermodynamic equilibrium nearing hydrogen exhaustion, it contracts and tends toward critical rotation at the surface (recall we are enforcing solid body rotation). At this limit, the area ratio reaches its minimum of about 0.945, which, per Eq. (41) and at constant luminosity, leads to an effective temperature ratio of 0.986, meaning a 2.4% difference. In the bulk of the main sequence, however, we have an area ratio of around 0.99, which corresponds to just a 0.25% difference in effective temperature.

Fig. 6. Hertzsprung-Russell diagram of the main sequence evolution of 5 M⊙ rotating stars in MESA using the default and new BCs. The inset zooms in on the blue hook near the terminal age main sequence.

Fig. 7. Fraction of equivalent spherical area to true area at the surface as a function of age of the single rotating 5 M⊙ model. A lower fraction indicates a further departure from spherical symmetry due to rotation. Top panel: fractional critical rotation rate for comparison.

5.3. Evolution of a detached B-type star

We modeled the evolution of a 5 M_⊙ solar metallicity B-type star in a binary with a point mass companion of 3 M_⊙ (q = 0.6) and with an initial period of three days until the star reaches its RL. With this setup we mimic the situation of intermediate mass stars found in short-period eclipsing binaries. We varied the rotation distortion model between (1) no rotation, (2) single star rotation where only centrifugal forces are taken into account, and (3) tidal deformation that account for both centrifugal and tidal forces. For the latter two, we also varied the atmospheric BCs between the new atmospheric BCs (Sect. 4) and the default spherical ones. For the non-rotating model, the new BCs reduce to the default ones, as there is no departure from spherical symmetry.

Figure 8 shows the radius evolution of the stars as a function of the central hydrogen fraction, X_c. At the onset of evolution, where the star is well within its RL, we see negligible difference between the various assumptions on BCs or rotation. However, as the star approaches its RL size, at constant X_c, the family of rotating tracks differ from the non-rotating model by up to 5%. This is a significant difference, and is resolvable by very precise observations of eclipsing binaries, whose errors on measured radii can be better than 3% (Torres et al. 2010; Graczyk et al. 2018; Pavlovski et al. 2018; Tkachenko et al. 2020, and references therein). Therefore, modeling stars in eclipsing binaries as non-rotating in evolution codes could induce a systematic mismatch when fitting measured radii to inferred radii from such models. Finally, we note that the difference in radii between the tidal and single rotating star models, which is at most 1%, is only resolvable by the most accurate of observations, if at all. Furthermore, the difference in radii between considered BC is on the order of 0.1% and is not perceivable on the figure.

Fig. 8. Radius evolution of a 5 M⊙ star starting on a 3-day orbit around a 3 M⊙ point mass, as a function of central hydrogen fraction Xc. The difference between the rotation models is small while the one resulting from the different BCs is imperceivable.

5.4. Evolution of a twin overcontact binary

We evolved a 32 M_⊙ star with a metallicity of Z = Z_⊙/2 in a binary system with its twin, meaning all properties of the stars are identical. In particular, this means that the mass ratio is unity, and we can ignore binary interaction effects like mass or energy transfer. We start the evolution at the ZAMS with a 1.5 day orbit and stop the simulation once the stars reach the outer Lagrangian points. This setup is the closest we can currently and consistently model massive overcontact binaries such as VFTS 352 in the Large Magellanic Cloud (Almeida et al. 2015).

In Fig. 9, we show the effective temperature evolution of the 32 M_⊙ model under the various rotation model and BC assumptions. We notice again, like in the detached case, that the observable stellar parameter (in this plot T_eff) can vary up to 5% in the late stages of the overcontact configuration. The difference between the rotation models are minor, on the order of 1%, and the tidally distorted model with the new BCs of Sect. 4 shows a rather complex evolution. It kinks at the RL filling point and subsequently crosses the other rotating tracks as the degree of overcontact increases. This behavior may be related to the crossing of the structure correction factors, f_P and f_T, as r_Ψ keeps rising (see Fig. 2). Moreover, there is a ∼25% age discrepancy between the rotating and non-rotating models. This difference can be understood in the context of spin orbit coupling. In the non-rotating model, the star can store no angular momentum, so that tidal synchronization has no effect (at least with the current implementation in MESA) and the RL size slowly increases due to wind mass loss (see top panel Fig. 9). However, when the star is modeled with angular momentum, wind mass loss causes an efficient spindown of the star. Tidal synchronization then keeps the star rotating close to the orbital velocity at the cost of decreasing the orbital separation, thus reducing the RL size and accelerating the points of contact and L₂ overflow. This effect is also visible in the detached binary of Sect. 5.3.

Fig. 9. Effective temperature evolution of the 32 M⊙ contact model with the different assumptions on BC and rotation model. We notice that the models including rotation are up to 5% cooler during overcontact. Between the different rotation models and BCs, we find only minor differences. Top panel: RL evolution, which explains the difference in age at onset of overflow of the Lagrangian points through spin orbit coupling.

6. Discussion and conclusions

We developed the methodology needed to account for tidal deformation in binary stars in 1D stellar evolution codes. They are represented in the structure-correction factors of f_P and f_T that multiply the equations of hydrostatic equilibrium and radiative energy transport, respectively. Additionally, modified expressions of the Eddington limit and atmospheric boundary conditions as a result of the deformation have been calculated.

We show that the radii predicted by stellar evolution models can shift by ∼5% if centrifugal deformation is included, while also including tidal deformation leads to a smaller ∼1% effect. This means comparing observed stellar parameters to non-rotating models could result in a mismatch and this could have implications for high-precision astrophysics. In a study where observed surface properties with uncertainties on the order of 5% are used, centrifugal deformation has to be included in the stellar models. In binaries, if the precision of measurement is better than 1%, then the tidal deformation has to be taken into account. We note that the effect observed here is only due to geometrical considerations, that is, the equator bulges out due to rapid rotation or the presence of a companion. Mixing as a result of rapid rotation was not included and will further affect the radius evolution. Generally, the effect of rotational mixing is to keep the stars compact as fresh fuel is introduced into the core, so this may somewhat cancel the geometrical effect. However, the exact interplay and outcome of these competing effects is far from established.

The methods developed here have a wide range of applicability as they can be used for any (semi-)detached binary of arbitrary mass ratio. This can be applied, for instance, in studying the effect of tidal forces in mass transfer stability, which can have a potential impact on the formation of gravitational wave sources through stable mass transfer (van den Heuvel et al. 2017; Marchant et al. 2021). There are, however, limitations to our model, most notably that the structure correction factors computed from the binary Roche potential only apply to fully tidally locked systems. The general case, where the rotation of the components is free, cannot be described with a conservative potential, and is thus not suitable to the treatment of Kippenhahn & Thomas (1970) and Endal & Sofia (1976) used here. Developing a global model of free rotation including tides is a complex problem, but a potential workaround would be to consider free rotation for stars well within their RLs and use the centrifugal deformation corrections of single stars, and then implementing a switch to the synchronized binary corrections once a star fills an appreciable amount of its RL.

For contact configurations, currently only the mass ratio of unity can be considered consistently, as the energy transfer in the overcontact layers, is zero. When we include the processes of mass and energy transfer in overcontact layers, which is the subject of the following paper in this series, the methods we developed in this work enable the consistent modeling of all overcontact binary systems.

Acknowledgments

M.F. thanks the Flemish research foundation (FWO, Fonds voor Wetenshappelijk Onderzoek) PhD fellowship No. 11H2421N for its support. P.M. acknowledges support from the FWO junior postdoctoral fellowship No. 12ZY520N. The research leading to these results has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement numbers 772225: MULTIPLES).

References

Appendix A: Polynomial fits in the rotating shell potential

In the case of the rotating star potential (18), integral expressions can be written down to compute all relevant quantities required to use the modified stellar structure equations presented in Sect. 2.1. Paxton et al. (2019) numerically integrated these and provided polynomial fits in the variable ω = Ω/Ω_crit, where is the classical critical velocity of equatorial radius r_e.

To improve on the errors introduced by such fits, we recalculated the integrals and used the Levenberg-Marquardt non-linear least squares algorithm from scipy (Virtanen et al. 2020) to find more accurate polynomial fits.

First, for the equipotential volumes V_Ψ, we have the fit:

(A.1)

with a maximal error of 0.20% in ω ∈ [0, 1]. We note that since V_Ψ(1) is an analytical result (Kopal 1959) and the leading order term in the expansion is set by requiring that a slow rotating shell is ellipsoidal, this fit has only one free parameter. For the volume equivalent radius, r_Ψ, and equatorial radius, r_e, have the following fits:

(A.2)

with a maximal error of 0.11% in ω ∈ [0, 1], and,

(A.3)

with a maximal error of 0.15% in ω ∈ [0, 1]. Next, the surface area S_Ψ and average effective gravity ⟨g⟩ are approximated by:

(A.4)

(A.5)

with a maximal error of 0.07% and 0.25% in ω ∈ [0, 1], respectively. Since the inverse gravity diverges when ω → 1, we require a non-polynomial term to represent it appropriately in that limit. Therefore, constructing the auxiliary functions:

(A.6)

(A.7)

we have then for the average inverse effective gravity ⟨g⁻¹⟩, and the specific moment of inertia i_rot, the following fits:

(A.8)

(A.9)

with a maximal error of 0.85% and 0.68% in ω ∈ [0, 0.9999], respectively.

For the correction factors, f_P and f_T, we construct the following fits:

(A.10)

(A.11)

with a maximal error of 0.56% and 0.58% in ω ∈ [0, 0.9999], respectively.

Finally, with the auxiliary function:

(A.12)

the implicit equation from which ω can be obtained from j_rot is:

(A.13)

with a maximal error of 0.54% in ω ∈ [0, 0.9999].

With these fits, we improve all the fits of Paxton et al. (2019), especially the maximal error of f_T, namely, taking it from 1.6% to 0.58%.

Appendix B: Matching the tidal and single rotating star distortion models

In Fig. 2, we make the comparison of the structure correction factor, f_P, of the tidal model as well as the single rotating star model in one plot. This is meant to illustrate the corrections that would be used in a simulation with the MESA code if rotation is included in a model for a binary system without including tidal deformation. In practice, a tidally distorted star does not have a well defined equatorial radius with which to define ω, but nevertheless the current implementation of MESA makes use of the equipotential radius r_Ψ, the mass contained within the shell M_Ψ and the angular frequency of the shell Ω_Ψ to derive the corrections as if these properties corresponded to a single rotating star with cylindrical symmetry.

For a tidally synchronized binary, the angular velocity is derived from Kepler’s third law:

(B.1)

with q = M₂/M₁ as the mass ratio and a as the separation of the binary. Next, we consider the equivalent single rotating star. Given the mass M₁ and equatorial radius r_e, it has a critical rotation rate:

(B.2)

Dividing these equations, and setting ω = Ω/Ω_crit as the dimensionless rotation rate, we arrive at an implicit equation for ω:

(B.3)

where f_re relates the volume equivalent radius to the equatorial radius as a function of the dimensionless rotation rate in a single rotating star. For a given value of r_Ψ/a and q, this implicit equation gives the value of ω that a MESA simulation would determine to compute the corrections to the stellar structure equations, despite there not being a well defined equatorial radius. The values shown in Fig. 2 for the single star equivalent correspond to evaluating the results of Appendix A for f_P(ω), with ω derived from r_Ψ/a and q with the implicit Eq. (B.3).

All Figures

	Fig. 1. Illustration of the Roche geometry in the equatorial plane. It shows the L1, L2, and L3 equipotentials (gray), as well as the splitting surfaces (blue dashed). The splitting surfaces have their normals always perpendicular to the potential gradient.
In the text

	Fig. 2. Structure correction factor fP as a function of fractional RL radius rΨ/rRL for various mass ratios. The discontinuities in all graphs correspond to the crossing of the Lagrangian points L1 (at rΨ/rRL = 1, by construction) and Lout at larger radii. For reference, the equivalent single rotating star inferred fP is overplotted in dashed lines.
In the text

	Fig. 3. Comparison of the calculations of the structure correction factor fP from Mochnacki (1984) and those in this work, as a function of the fractional RL radius rΨ/rRL. We performed a resolution convergence test with double, quadruple and octuple the default spacial resolution, and conclude our integrations have well converged at an eight-fold resolution increase.
In the text

	Fig. 4. Maximal Eddington factor Γmax that a rotating star can have as a function of fractional classical critical rotation rate of ω = Ω/Ωcrit, class. Our model of shellular stars has a maximal reduction to Γmax = 0.639.
In the text

	Fig. 5. Reduction, fP/fT, of the maximal Eddington factor as a function of mass ratio, q, and fractional overflow radius, rΨ/rRL. The equivalent radius of a star filling up to its outer Lagrangian point in marked with the green line.
In the text

	Fig. 6. Hertzsprung-Russell diagram of the main sequence evolution of 5 M⊙ rotating stars in MESA using the default and new BCs. The inset zooms in on the blue hook near the terminal age main sequence.
In the text

	Fig. 7. Fraction of equivalent spherical area to true area at the surface as a function of age of the single rotating 5 M⊙ model. A lower fraction indicates a further departure from spherical symmetry due to rotation. Top panel: fractional critical rotation rate for comparison.
In the text

	Fig. 8. Radius evolution of a 5 M⊙ star starting on a 3-day orbit around a 3 M⊙ point mass, as a function of central hydrogen fraction Xc. The difference between the rotation models is small while the one resulting from the different BCs is imperceivable.
In the text

Fig. 9. Effective temperature evolution of the 32 M⊙ contact model with the different assumptions on BC and rotation model. We notice that the models including rotation are up to 5% cooler during overcontact. Between the different rotation models and BCs, we find only minor differences. Top panel: RL evolution, which explains the difference in age at onset of overflow of the Lagrangian points through spin orbit coupling.

In the text