Issue 
A&A
Volume 672, April 2023



Article Number  A175  
Number of page(s)  11  
Section  Stellar structure and evolution  
DOI  https://doi.org/10.1051/00046361/202346277  
Published online  19 April 2023 
Modeling contact binaries
II. Effects of energy transfer
^{1}
Institute of Astronomy (IvS), KU Leuven, Celestijnenlaan 200D, 3001 Leuven, Belgium
email: matthias.fabry@kuleuven.be
^{2}
ArgelanderInstitut für Astronomie, Universität Bonn, Auf dem Hügel 71, 53121 Bonn, Germany
^{3}
MaxPlanckInstitut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
Received:
28
February
2023
Accepted:
15
March
2023
Context. It is common for massive stars to engage in binary interactions. In close binaries, the components can enter a contact phase, when both stars simultaneously overflow their respective Roche lobes. While observational constraints on the stellar properties of such systems exist, the most detailed stellar evolution models that feature a contact phase are not fully reconcilable with those measurements.
Aims. We aim to consistently model the contact phases of binary stars in a 1D stellar evolution code. To this end, we have developed a methodology to account for energy transfer in the common contact layers.
Methods. We implemented an approximative model for energy transfer between the components of a contact binary based on the von Zeipel theorem in the stellar evolution code MESA. We compared structure and evolution models both with and without this transfer. We then analyzed the implications for the observable properties of the contact phase.
Results. Implementing energy transfer helps in eliminating baroclinicity in the common envelope between the components of a contact binary, which (if present) would drive strong thermal flows. We find that accounting for energy transfer in massive contact binaries significantly alters the mass ratio evolution and can extend the lifetime of an unequal mass ratio contact system.
Key words: stars: evolution / binaries: close
© 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
Massive stars (M_{initial} ≳ 8 M_{⊙}) are predominantly found in binaries that interact over the course of their lifetimes (Sana et al. 2012). This means that the impact of binarity on the evolution of massive stars must be taken into account in stellar structure and evolution models, as evidenced by population synthesis calculations (Vanbeveren et al. 1998; de Mink et al. 2014). An extreme form of binary interaction happens when the stellar companions enter a contact configuration, where both stars overflow their respective Roche lobes (RLs). The stars connect at the socalled neck and form a peanutlike shape, as determined by the geometry of the Roche potential.
Understanding the evolution of contact systems is of great importance in the study of stellar mergers (among others), as these are all preluded by an (unstable) contact phase. Mergers could be the potential source of strong magnetic fields observed in about 7% of massive stars (Ferrario et al. 2009; Wickramasinghe et al. 2014; Fossati et al. 2015; Schneider et al. 2016, 2019; Grunhut et al. 2017; Frost et al. 2023). Furthermore, close binaries are expected to rotate rapidly due to the high Keplerian velocity and quick tidal synchronization. This might lead to the physical conditions necessary to create longduration gammaray bursts and magnetardriven superluminous supernovae (Yoon & Langer 2005; Yoon et al. 2006; AguileraDena et al. 2018, 2020). The high rotational velocity could also lead to chemically homogeneous evolution, which is a proposed channel to produce gravitational wave progenitors (Marchant et al. 2016; Mandel & de Mink 2016).
In the past, lowmass contact systems, known as W Ursae Majoris (W UMa) systems, have been extensively studied. These systems have an orbital period of less than a day with a total mass on the order of 1 M_{⊙}. Observational campaigns from Eggen (1967) and Binnendijk (1970) identified the first of these stars and found that a significant fraction had mass ratios away from unity. At present, thousands of W UMa systems have been identified thanks to large surveys such as the OGLE survey (Szymanski et al. 2001). Observations of massive contact systems however are much rarer, but fortunately, two dozen or so systems have been identified from the VLT Flames Tarantula Survey (VFTS; Evans et al. 2011), MACHO survey (Alcock et al. 1997), and OGLE in the Magellanic clouds, as well as studies of galactic targets (e.g., Lorenzo et al. 2014, 2017; Yang et al. 2019).
Also on the modeling side, most work in the literature is dedicated to understanding the W UMa systems as opposed to massive contact systems. Kuiper (1941) initiated theoretical considerations of contact binaries based on the argument that due to the stark difference between RL radii and the massradius relation of main sequence stars, contact systems with mass ratios that are away from unity cannot be stable. The Kuiper argument eventually transformed into what is now referred to as Kuiper’s paradox once a multitude of contact systems were observed to have unequal mass components. However, it should be immediately clear that this paradox is only apparent as the massradius relation of single main sequence stars need not apply for stars in a contact configuration.
Lucy (1968) first considered energy transfer (ET) in common convective envelopes of stellar components, and computed the first approximated structure models for W UMa stars. Later, Lucy (1976), Flannery (1976), Hazlehurst (1985) and Kähler (1989) calculated models that are out of thermal equilibrium and show cyclic behavior. Shu et al. (1976, 1979) and Lubow & Shu (1977, 1979) (collectively SLA) constructed models of contact binaries where they dropped the requirement of a continuous structure at the layer coinciding with the equipotential surface of the first Lagrangian point (L_{1}). These models resolve Kuiper’s paradox regardless of the thermal structure of the envelopes, namely, whether they are radiative or convective. Despite much criticism (see e.g., Hazlehurst 1993; Kähler 1989), this is the simplest model of ET in radiative envelopes available in the literature.
Accurate modeling of massive, longlived contact systems has been attempted in the past. Marchant et al. (2016) computed detailed evolution model grids of massive binaries with initial periods down to 0.5 d, which includes the regime of contact binaries at the zero age main sequence. Those models however did not include energy transfer between contact components, as the models concerned binaries of mass ratio close to unity M_{2}/M_{1} = 0.8 − 1 and the effect was thought to be minimal.
Sen et al. (2022) used the models of Marchant (2016) to study the semidetached Algol systems, although these models also included contact phases. Menon et al. (2021) extended these grids by computing models with initial mass ratios down to M_{2}/M_{1} = 0.6 with the study of massive, longlived contact binaries in mind. They found, also without including energy transfer in contact phases, a strong correlation between observed mass ratio and period in contact systems that is broadly in agreement with observations. However, the mass ratio distribution they derive is heavily skewed toward values close to unity, which is not supported observationally. These authors suggested that including energy transfer in contact phases of unequal mass components could alleviate this discrepancy.
In this work, we apply the ET model of SLA in common stellar layers to modern stellar structure models, to further advance our evolutionary modeling of massive contact binaries. In Sect. 2, we describe and discuss the theory of ET. Section 3 describes the physical setup of the stellar evolution code, along with our ET implementation. In Sect. 4, we compare models computed with and without ET in contact layers and discuss the difference in the observable properties of the models. Lastly, in Sect. 5, we provide our concluding remarks.
2. Theory of energy transfer
2.1. Simple considerations
Theoretical modeling of contact stars started with Kuiper (1941), who stated that stable contact systems of uniform composition with unequal masses cannot exist as a result of differing massradius relationships. For the galactic zero age main sequence (ZAMS) models of Brott et al. (2011), we derive the mass radius relation of single stars to be approximately:
where we defined q ≡ M_{2}/M_{1}. However, the condition of contact in a binary following the Roche geometry constrains the surface of the stars to same equipotential, which leads to:
for stars not overflowing their RL to excess. Clearly, Eqs. (1) and (2) cannot be satisfied simultaneously unless M_{1} = M_{2}. However, while the contact condition needs to be satisfied from dynamical arguments, the stellar structure does not, a priori, need to follow a single star model. Even though the dense stellar core will be largely unperturbed due to the companion, the outer layers of contact components are highly distorted, causing variations in the total radii with respect to spherical models, as seen explicitly in Fabry et al. (2022, henceforth Paper I). This result had not yet taken the possibility of energy transfer into account and we expect further changes under the consideration that a hotter gas under similar pressure takes up more volume. Therefore, the ZAMS massradius relation of Eq. (1) is not expected to be satisfied under general contact conditions, and so Kuiper’s paradox can be resolved by providing alternative stellar models that have a massradius relation closer to the contact condition of Eq. (2).
Following the Roche lobe geometry, Lucy (1968) finds the ratio of the surface areas S_{2}/S_{1} of two stars in contact to be proportional to (M_{2}/M_{1})^{β}, with β = 0.96. Using the approximation β ≃ 1, combined with von Zeipel’s theorem of gravity darkening, von Zeipel (1924), and the StefanBoltzmann law, , results in the simple expectation that in contact binaries, the luminosity ratio follows the mass ratio (Lucy 1968; Tassoul 2000):
Single main sequence stars, on the other hand, follow the well known massluminosity relation:
with α ≃ 2 − 3 for the upper main sequence (Köhler et al. 2015; Gräfener et al. 2011). We use the symbol ≃ to denote these relations are approximations to simple power laws. We see that this leads to a difference between the luminosity of a single star L_{s, 1} and the luminosity L_{1} of a star of the same mass in a contact binary of
and for the companion
since L_{2} ≃ qL_{1} and we require ΔL_{1} + ΔL_{2} = 0 to conserve energy. This defines f as
Therefore, the two stars in a contact binary can fulfill the single star massluminosity relation of Eq. (4) in their cores and the contact binary massluminosity condition of Eq. (3) at their surfaces if the amount of energy per time given by Eqs. (5) and (6) is transferred from the more massive to the less massive star in their common envelope.
2.2. Models of energy transfer
The general solution to Kuiper’s paradox is to consider detailed stellar models with the inclusion of energy transfer (ET) between the binary components. Several models of ET are given in the literature.
Lucy (1968) and Biermann & Thomas (1972) provided a first solution by adjusting the adiabatic constants of convective envelopes in contact components. However, this is an unsatisfactory solution, since this setup requires the stars to be burning hydrogen through different nuclear chains or cycles in the case of Lucy (1968) or that the models exhibit inaccurate light curves, as in Biermann & Thomas (1972). Other models, such as those of Lucy (1976) or Flannery (1976), relaxed the requirement of thermal equilibrium and constructed models of W UMa stars that exhibited thermal cycles. Kähler (1989) presented a detailed model that required turbulent motions in the common envelope to explain early type (radiative) W UMa binaries.
Meanwhile, SLA presented the contact discontinuity model of contact binaries, by relaxing the requirement of continuous structural quantities across the RL. This is the only model that treats the common envelope as a single volume of the binary structure, at the price of hiding a heat engine in a very thin region around the RL. One peculiar feature is that this model necessitated a temperature inversion at the RL layer in one of the components as otherwise they would not be able to construct thermally stable contact models of uniform composition (in order to model binaries at zero age). This feature has received criticism in that the proposed heat engine violates the second law of thermodynamics and cannot be stable over thermal timescales (see Hazlehurst 1993; Kähler 1989, and references therein). Eventually, Kähler (2004) drew the following conclusion based on all collected theoretical arguments: internal circulation currents must exist in the less luminous component to reduce the radiative temperature gradient, since the luminosity carried by radiation is reduced by the circulation luminosity.
Given the complexity of the theoretical problem of the structure of contact binaries, especially with radiative envelopes, it is beyond the scope of this work to further develop the analytic theory. Instead, we apply an ET model in modern stellar structure calculations. Using shellularity as our base assumption of the stellar structure (see Sect. 3.2 for the precise definition), the model of SLA is a natural choice, although we recognize that this comes with the apparent thermodynamical problems stated above. However, we believe we avoid the most fundamental one, as we do not explicitly require a temperature inversion in our models. We only used the model of SLA to compute the amount of energy transferred (see Sect. 2.3).
2.3. Energy transfer in the Roche geometry
The work of SLA provides a general model of ET by introducing the notion of an energy flow at the base of the common envelope. If the contact layers are shellular, and they satisfy von Zeipel’s gravity darkening, conservation of energy at the RL implies:
Here, S is the surface area of the RL, ⟨g⟩ is the surface averaged effective gravity at the RL, and the primed quantities specify the state just above the ET layer, while unprimed those just below. This equation specifies that the fraction of the total luminosity that each component radiates is proportional to S⟨g⟩. The transferred luminosity then equals
Comparing Eqs. (8a) and (8b) against Eqs. (5) and (6), we find for the fraction f:
which is consistent with Eq. (7) as and .
3. Methods
To investigate the effect of ET on the evolution of contact binaries, we computed binary evolution models using the stellar evolution code MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023), version r22.11.1. These models are the first massive binary evolution models that include ET in contact layers. We followed the evolution of binaries from the ZAMS until the least massive component overflows the second Lagrangian point.
3.1. Physical assumptions in MESA
The microphysical setup of our stellar evolution models remains mostly the same as in Paper I, with some additions from the newer MESA version. In summary, the nuclear net contains the eight isotopes of ^{1}H, ^{3}He, ^{4}He, ^{12}C, ^{14}N, ^{16}O, ^{20}Ne, and ^{24}Mg, sufficient for the main sequence and nuclear burning rates are taken from the JINA (Cyburt et al. 2010) and NACRE libraries, with weak interaction rates from Fuller et al. (1985), Oda et al. (1994) and Langanke & MartínezPinedo (2000). Plasma screening is included following Chugunov et al. (2007) and thermal neutrino losses are computed in Itoh et al. (1996). For the equation of state, a blend is used between different tables by Saumon et al. (1995), Timmes & Swesty (2000), Rogers & Nayfonov (2002), Potekhin & Chabrier (2010) and Jermyn et al. (2021), specified in Jermyn et al. (2023). Radiative opacities are also blended from COenhanced tables of OPAL opacities (Iglesias & Rogers 1993, 1996) and Ferguson et al. (2005) for lower temperatures. At high temperatures, Compton scattering opacity is from Poutanen (2017), and electron conduction opacities are taken from Cassisi et al. (2007) and Blouin et al. (2020). In all simulations in this work, we set the metallicity of stars to the solar value, Z_{⊙}, where Z_{⊙} = 0.0142, with metal fractions as determined from Asplund et al. (2009).
Mass loss through winds is accounted for by following the prescription of Brott et al. (2011). If the surface hydrogen fraction is X > 0.7, the mass loss rate is taken either from Vink et al. (2001), for temperatures above the iron bistability jump (calibrated also by Vink et al. 2001), or the maximum of the rates from Vink et al. (2001) and Nieuwenhuijzen & de Jager (1995), below this temperature. For X < 0.4, the wind prescription of Hamann et al. (1995) is used, albeit decreased by a factor of ten. When 0.4 < X < 0.7, the wind is linearly interpolated between the above results. We allow part of the wind launched by a star to be accreted by its companion using the BondiHoyle mechanism (Bondi & Hoyle 1944) as implemented by Hurley et al. (2002).
For mixing, we use the Ledoux criterion to determine convectively mixed regions, where the mixing length theory of BöhmVitense (1958), in the version described by Cox & Giuli (1968), is applied with a mixing length parameter of α = 2. The convective core is allowed to overshoot its boundary. Following the calibration of Brott et al. (2011), we use a step overshoot where the diffusion coefficient 0.01 pressure scale heights into the convective layer is kept constant out to 0.335 pressure scale heights beyond the boundary. Semiconvection is included following the model of Langer et al. (1983) with a high efficiency of α_{sc} = 100 as calibrated by Schootemeijer et al. (2019). We include thermohaline mixing as developed by Kippenhahn et al. (1980) with an efficiency parameter α_{th} = 1. All other mixing processes, in particular rotational mixing, are ignored in order to isolate (as much as possible) the effect of energy transfer.
At all times, the rotation of the stars is synchronized to the orbital period, as part of the requirement to use the Roche potential as the deformation geometry (Sect. 3.2). Additionally, we artificially diffuse the total angular momentum throughout the interior of the star to enforce solid body rotation. The treatment of mass transfer (MT) is explained in Sect. 3.3, while the implementation of energy transfer (ET) is shown in Sect. 3.4.
3.2. Shellularity and Roche lobe geometry
Since we deal with highly tidally deformed stars, we use the modifications to the stellar structure equations from Paper I to incorporate the RL geometry into a onedimensional (1D) stellar evolution code. The stars are therefore modeled as hydrostatic structures living in the Roche potential Ψ of a fully synchronized binary (see Eq. (19) of Paper I). Furthermore, we assume stellar layers to be shellular. Shellularity is reached when all intensive quantities, (in particular the temperature, pressure, and mass density) are constant along a stellar layer and that such a layer coincides with a unique equipotential surface.
It should be emphasized that in the context of 1D models, the full shellularity of layers of a (single) star is an assumption – not a selfconsistently modeled feature as, of course, there is no 3D structure to explicitly test the shellularity of a stellar layer. However, overflowing layers of components in a contact binary can be tested whether they are shellular with respect to each other if the temperature, density, and so on at the 1D cells are equal for equal values of Roche potential. Given the method we used to split the common envelope of the contact binary (Sect. 2.4 of Paper I), the computation of the tidal deformation corrections in Paper I based on equipotential surfaces and the usage of the corresponding outer boundary condition (Sect. 4 of Paper I), the final ingredient of ET in the common envelope developed here will ensure the shellularity of one component with the other – again, defined as having similar cell values of pressure, temperature, and so on for similar values of potential.
In contact binaries, the outer stellar layers are shared between the components and a choice must be made as to what part of the envelope belongs to each component. To this end we constructed what we call a splitting surface emanating from the Lagrangian point L_{1}, by following the gradient ∇Ψ outward in all directions (see Paper I, Sect. 2.4). With this setup, the volume of both components is uniquely defined, and in Paper I we computed integrals in the Roche potential necessary to represent distorted shells in 1D. The volume equivalent radius then acts as the new independent variable of a stellar shell, and is defined as:
Given a certain splitting strategy, and the requirement of the stellar surface being on a common equipotential, a relation of the form:
with q = M_{2}/M_{1}, constraining the radii of the overflowing components 1 and 2. This is a statement of the contact condition and it is at the basis of Kuiper’s paradox (Sect. 2). For vertical splitting surfaces through L_{1}, Marchant et al. (2016) found the following fit:
where we have denoted r_{RL} = r_{ΨL1} the volume equivalent RL radius of a component. With our setup of the splitting surface however, this approximation is no longer accurate for significantly overflowing shells of mass ratios away from unity. Therefore, in this work, we evaluate the function F of Eq. (12) by interpolating the results of the Roche integrations obtained in Paper I.
Retaining the form of the fit of Marchant et al. (2021), for the radius of a component overflowing to its outer Lagrangian point, L_{out}, which is L_{2} for the less massive and L_{3} for the more massive component, we constructed a new fit from the integration results of Paper I, which is:
This fit has an error smaller than 0.1% in the range −7 ≤ log q ≤ 7, compared to 1% if we used the fit of Marchant et al. (2021). We use this fit then to determine when the least massive component overflows to L_{2}, after which we stop the evolutionary simulation.
3.3. Mass transfer
The components of a contact binary are, per definition, in mechanical contact with each other. Thus, stellar material can be exchanged throughout the evolution of such a system. However, MT cannot occur at an arbitrary rate, since the surface of both components are constrained by the contact condition of Eq. (12). If this were not satisfied, pressure gradients would arise across the surface of the contact binary, and we expect this would induce strong, horizontal^{1} flows equilibrating the surface. Thus, in our simulations, we implicitly adjust the MT rate either to keep the donor just below its RL radius in a semidetached configuration, or so the component radii satisfy the contact condition of Eq. (12) in a contact configuration (see also Marchant et al. 2016).
3.4. Energy transfer
Not only are the contact layers of both stellar components in mechanical contact, facilitating MT, there is also thermal contact between these layers, allowing for ET. As elaborated in Sect. 2.3, we model ET in contact binaries as the transfer of luminosity of one component to the other (Eq. (9)), at the location of the RL. It should be mentioned that no further corrections to the ET need to be computed since the splitting surface we use, computed in Paper I, is parallel to the local effective gravity and thus the radiative flux. Therefore, since the regular energy transport by radiation (or by convection) in the interior of either component is vertical, no vertical energy flux crosses over toward the other component. Our splitting strategy thus nicely decouples horizontal and vertical energy flows and in this approximation, the ET computed here is not affected by the vertical energy transport present in the components.
The stability and performance of numerical simulations are severely impacted by introducing sharp discontinuities at the RL radius; in this case, that of luminosity. To remedy these issues, we used the following approach.
First, for numerical stability over evolutionary steps, we smooth the amount of ET by computing:
where is the smoothed transferred luminosity at an evolutionary step, n, L_{trans} is the luminosity to be transferred as per the Roche geometry (Eq. (9)), and p is the smoothing factor. In all simulations including ET, we use a moderate smoothing of p = 0.5.
Then, we implement the luminosity transfer as a constant extra source of specific heat ε_{RL, 1, 2} in both stars, occurring in the cells with a radius within 1.00 and 1.01 times its respective RL radius, a choice meant to reflect the estimate of Shu et al. (1979) (see Eq. (18) below and also Appendix A for further details). The magnitude of ε_{extra, 1, 2} is determined from the required luminosity:
where Δm_{1, 2} is the mass of the cells we put the extra heat in.
Since MT is modeled as removing stellar material from the top of the donor and putting it on top of the accretor, during MT, the common layers are out of thermal equilibrium due to expansion or compression. To correct for this effect and regain shellularity, especially at the surface, we transport additional heat in the common envelope according to Eq. (9), where now we evaluate S⟨g⟩ at the surface layer of both components instead of the RL. The constant heat source ε_{surf, 1, 2} is then computed similarly to Eqs. (16a) and (16b), only now Δm encompasses all shells from 1.01 times the RL up to the surface.
As a final approximation, we linearly scale the necessary ET for shallow contact systems, defined as when for either star:
We recognize that using our ET prescription in fast MT phases is the weakest element of our simulations. When the MT timescale is shorter than the thermal timescale, the star falls out of thermal equilibrium, which is characterized by strong (vertical) gradients in the luminosity profile. During such phases, the energy budget of the star is redistributed vertically in order to regain as best as possible thermal equilibrium. Adding to this process the horizontal energy transport in the layers near the RL could induce complex interactions between both energy flows. Even though we mention the horizontal and vertical energy flows to be decoupled above, this is an idealization and the approximation may break down in quickly evolving phases. Furthermore, while our computation of how much energy is transferred depends only on the geometry and the theorem of von Zeipel, which does not require thermal equilibrium of the layers, hydrodynamical effects break the validity of von Zeipel’s theorem (von Zeipel 1924), which can be significant near the Lagrangian point, L_{1}.
4. Stellar models
4.1. Detailed example of ET evolution
Here, we detail the evolution of a 25 M_{⊙} primary and a 20 M_{⊙} secondary (so q_{ini} = 0.8) in an initial orbit of 1.4 days. This initial configuration is chosen such that, at the onset of the long lasting contact phase, this model mimics the observed galactic contact binary V382 Cyg studied in AbdulMasih et al. (2021). We simulate both cases where ET, as implemented in Sect. 3.4, is included and excluded.
In Fig. 1, we show the Hertzsprung–Russel diagram (HRD) broken up into four distinct phases of the evolution, while in Fig. 2, we show the MT rate and mass ratio evolution. First, we have the preinteraction phase (leftmost column) followed by a short lived, fast case A MT phase from the primary to the secondary when the donor experiences RL overflow (second column of Fig. 1). During this MT phase, the components briefly enter into contact, as the accretor falls highly out of thermal equilibrium and swells beyond L_{1}. After thermal relaxation of the accretor, the system detaches, after which slow case A MT is initiated (third column of Fig. 1), still from the primary to the secondary, although the mass ratio has inverted by this point. This phase of nuclear timescale MT occurs in a semidetached configuration and is referred to as an Algol phase. Finally, the evolution of the now more massive secondary catches up to that of the stripped primary, contact is engaged again and MT inverts (rightmost column of Fig. 1), moving material from the secondary back to the primary (however see below). This contact configuration evolves on the nuclear timescale and ends when L_{2} overflow occurs.
Fig. 1. HRD of a 25 + 20 M_{⊙} binary with an initial period of 1.4 d. Red lines indicate a simulation with the inclusion of ET, while blue lines ignore ET. Moving left to right, the columns highlight four subsequent phases in the evolution (see Sect. 4.1 for details). Contact phases are highlighted with thick outlines. 
Qualitatively, the evolution of this system is similar with and without inclusion of ET. Up to and including the slow case A MT phase, the evolution of this system is nearly identical. One exception is the different HRD tracks during the fast MT case. In the second column of Fig. 1, the red tracks exhibit a “hook”, while the blue tracks do not. This is a direct consequence of ET in the short contact phase changing the outer structure of the components. The secondary donates energy to the primary and cools a bit, jumping to the right of the HRD (and viceversa for the primary). The point at which the primary and secondary tracks touch (at roughly log L/L_{⊙} = 4.8 and log T_{eff}/K = 4.545) corresponds to the situation where M_{1} = M_{2}. This is consistent with the theory since we expect that for M_{1} = M_{2} we have similar luminosities per Eq. (3) and also similar effective temperatures, following the law of StefanBoltzmann. When contact is disengaged, the secondary (at this point more massive than the primary) jumps back to the left as the energy sink to the primary disappears. Given the timescale of this MT phase, however, it is unlikely to observe contact binaries in this phase and thus to constrain the effects of ET.
Still, we do observe quantitative differences between the energy transferring and nonenergytransferring models during the final, inverted MT phase, which is longlived and therefore more likely to be observed. In particular, from the q(t) curves in the bottom panel of Fig. 2, we see that the ET model spends more time at mass ratios q ≳ 1.4 versus the no ET model. Furthermore, the noET model monotonically approaches a mass ratio of unity without further inverting the mass ratio or the MT direction. In contrast, the ET model experiences two inversions of the mass ratio at around τ = 3.2 and 4.0 Myr and the MT is inverted at around τ = 3.3 Myr.
Fig. 2. MT rate (top panel) and mass ratio (bottom panel) evolution of the 25 + 20 M_{⊙} binary studied in Sect. 4.1. Contact phases are highlighted with thick outlines. 
In Fig. 3 (top panel), we show the differential duration of the longlived contact phase, plotted against ^{2}, while the cumulative duration is shown in the bottom panel. The top panel shows three distinct peaks in the red distribution, two of which correspond to the periods where the ET model is inverting its MT direction (at q ≈ 1.4 and q ≈ 0.85) and the last corresponds to a mass ratio close unity at the end of the simulation. Comparing the two former peaks to the no ET curve at those shows that the ET model spends more time at these mass ratios than the no ET model. As a cumulative distribution, from the bottom panel of Fig. 3, we deduce that it is about four times more likely to observe the ET model in a configuration more extreme than than the no ET model, even though the total length of the contact phase is similar to within 0.05 Myr.
Fig. 3. Contact duration distributions as function of observed mass ratio . Top panel: differential duration (in bins of ). Bottom panel: cumulative duration in configurations more extreme than . 
When comparing the profiles of the components in both models, we note a considerable difference in the structure of the outer layers of the components. In Fig. 4, the temperature, density, and luminosity profiles of the models at the start of inverted MT phase (τ = 2.7 Myr) are shown. The independent coordinate is chosen to be , so as to map the overflowing layers to the interval [0, 1]. At this point in the evolution, in both cases of ET and no ET, the stars are in a contact configuration and the primary is stripped to M_{1} = 18.3 M_{⊙}, while the secondary accreted mass to M_{2} = 26.4 M_{⊙}. We note that these masses are close to the component masses AbdulMasih et al. (2021) found for V382 Cyg. We see that in the system without ET, the overflowing shells do not satisfy our definition of shellularity since the temperature and density profiles of the primary and secondary do not match for constant Roche potential (which approximately corresponds to constant ). The ET models however match significantly better, as seen in the left column of Fig. 4 by the joining of the temperature and density profiles of the primary and secondary. In particular, at the surface, the density of the ET components agree to within 0.2%, while the temperature to within 0.5%, whereas the surface properties of the no ET components vary more than 10% (as expected for models of differing mass).
Fig. 4. Profiles of outer layers of the binary components at the onset of inverse MT. It shows (top to bottom) the temperature, density, and luminosity profiles for both components in the ET (red, left column) and no ET (blue, right column) cases, all as function of the scaled radius coordinate . The gray vertical lines show the location of the RL. 
An significant shortcoming of our model is the assumed thickness of the ET layer. Shu et al. (1979) argued that the thickness d of the ET layer is on the order of:
with a as the binary separation, and ρ, h, and c_{s} as the density, specific enthalpy and local sound speed evaluated at the RL, respectively. Lubow & Shu (1979) computed this number to be d/a ∼ 10^{−2} for stars of masses around 4 − 8 M_{⊙} which justifies SLA in modeling the layer as a discontinuity in the stellar profile, located at the RL radius.
However, direct computation of Eq. (23) for our 25 + 20 M_{⊙} model show that this estimation breaks down for higher masses, see the red line in Fig. 5. This suggests that the energy redistribution flow, modeled as a discontinuity by SLA, is not sufficient in these higher mass stars. Except when the binary has reached considerable overflow, where d/a ≲ 10^{−2}, we find that the thickness needed can be a significant fraction of the binary separation, even surpassing it in the early stages of the contact phase.
Fig. 5. Thickness of the ET layer, d, with respect to the binary separation a as function of the radius of the primary star during the long lived contact phase. The gray line gives the physical size of the overflowing layers, and corresponds to the maximal width the ET layer can assume. 
Another way to compute an upper limit on the thickness of the ET layer is to use Bernouilli’s equation. If we consider fluid motion from far away from L_{1} on the primary star (location i) toward L_{1} (location f), we have
where we have already canceled the potential terms Ψ_{i}, Ψ_{f} since we move along an equipotential surface, and the initial velocity, v_{i}, is assumed to be negligible. Making the estimation:
and assuming that the transferred energy through the binary neck of width, b, by a mass flow of Ṁ = ρ_{i}v_{f}b^{2} is:
we compute for the minimal thickness of the ET layer:
Finally the thickness of the layer at the neck is related to the thickness far away from L_{1} via:
since the Roche potential varies quadratically near L_{1} and linearly elsewhere. Equation (22) gives a lower limit on the width of ET layer so that a balanced mass flow, Ṁ, in the contact binary can carry a to be transferred luminosity, L_{trans}. Conversely, it can be interpreted as L_{trans} being the maximal luminosity the mass flow can carry in a layer of fixed width, b.
Also, we plot the Bernouilli computed thickness of Eq. (23) in Fig. 5. Similarly, we see that the required thickness, d/a, is much larger than what the radius of the primary star allows room for. It is only at later times, when the mass ratio has equilibrated and deeper contact is engaged, that the estimated width smaller than the overflow rate of R − R_{RL} of the primary.
4.2. Mass versus luminosity ratios
During the nuclear timescale, inverted MT phase, contact is engaged so that our ET scheme acts to move luminosity from one component to the other. As mentioned in Sect. 2, we expect the luminosity ratio of contact binaries to follow the mass ratio, L ∝ M, as opposed to detached stars following a single star massluminosity relation, L ∝ M^{α}, with α ≃ 2 − 3. Figure 6 shows the evolution of the luminosity ratio as function of the mass ratio during the slow MT phase of the 25 + 20 M_{⊙} explored in Sect. 4.1. In this graph, the models evolve from the top right at q ≈ 1.4 to near equal mass ratio on the left. We see that the model not including ET follow closely a q^{2.2} relation, appropriate for single stars in the mass range of 10 − 30 M_{⊙}. As the models including ET engage into deep contact however, their massluminosity ratio changes drastically from the q^{2.2} line in near contact to the q^{1} relation in full contact.
Fig. 6. Luminosity ratio versus mass ratio during the nuclear timescale MT of the 19 + 14 M_{⊙} systems of Sect. 4.1. Measurements from observed massive (near)contact systems from AbdulMasih et al. (2021), Yang et al. (2019), and Lorenzo et al. (2014) are overplotted. The systems from Mahy et al. (2020) were classified as ‘uncertain configurations’. 
We give (overplotted on Fig. 6) the measurements of several observed massive contact binaries of AbdulMasih et al. (2021), Yang et al. (2019) and Lorenzo et al. (2014) (see also Fig. 2 of Langer 2022). Curiously, the luminosity ratios from AbdulMasih et al. (2021), while they do follow an L ∝ M trend, are offset to lower L_{2}/L_{1} than predicted. Either the luminosity of the primary is overestimated or the uncertainties are underestimated. The systems included from Mahy et al. (2020) were categorized as ‘uncertain configurations’ since the measurement of the radius was consistent with being above as well as below the RL. In context of the massluminosity relation however, we expect that the system from Mahy et al. (2020) at q ≈ 1.3 (VFTS 563) is a true contact system, while the one at q ≈ 1.2 (VFTS 217) is not (although within 1σ it could be either).
5. Conclusions
In this work, we have taken a step forward in the detailed modeling of contact binaries, by implementing a model of ET in detailed stellar structure and evolution models. From Fig. 4, we see that this relatively simple model (i.e., the inclusion of a heat source or sink in the stellar model as proxy for the ET) is capable of shellularizing the common layers of stellar components in a contact configuration. Models without such ET do not exhibit their common layers to be shellular, which, from theoretical arguments, would drive strong horizontal flows equilibrating all gradients, in particular pressure.
From Figs. 2 and 3, we see that the time spent in deep contact at mass ratios between is extended when ET was included, versus when it was ignored. This is a promising result, in that if this trend persists across the parameter space of the total mass, initial mass ratio, and initial period, ET could provide an answer to the discrepancy between the observed mass ratio distribution of massive contact systems and its predicted distribution. The computation of a full grid of models and performance of a population synthesis is the topic of future work.
Acknowledgments
M.F. thanks the Flemish research foundation (FWO, Fonds voor Wetenschappelijk Onderzoek) PhD fellowship No. 11H2421N for its support. P.M. acknowledges support from the FWO junior postdoctoral fellowship No. 12ZY520N and the senior postdoctoral fellowship No. 12ZY523N. 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
 AbdulMasih, M., Sana, H., Hawcroft, C., et al. 2021, A&A, 651, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AguileraDena, D. R., Langer, N., Moriya, T. J., & Schootemeijer, A. 2018, ApJ, 858, 115 [NASA ADS] [CrossRef] [Google Scholar]
 AguileraDena, D. R., Langer, N., Antoniadis, J., & Müller, B. 2020, ApJ, 901, 114 [NASA ADS] [CrossRef] [Google Scholar]
 Alcock, C., Allsman, R. A., Alves, D., et al. 1997, AJ, 114, 326 [NASA ADS] [CrossRef] [Google Scholar]
 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Biermann, P., & Thomas, H. C. 1972, A&A, 16, 60 [NASA ADS] [Google Scholar]
 Binnendijk, L. 1970, Vistas Astron., 12, 217 [NASA ADS] [CrossRef] [Google Scholar]
 Blouin, S., Shaffer, N. R., Saumon, D., & Starrett, C. E. 2020, ApJ, 899, 46 [NASA ADS] [CrossRef] [Google Scholar]
 BöhmVitense, E. 1958, ZAp, 46, 108 [NASA ADS] [Google Scholar]
 Bondi, H., & Hoyle, F. 1944, MNRAS, 104, 273 [Google Scholar]
 Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cassisi, S., Potekhin, A. Y., Pietrinferni, A., Catelan, M., & Salaris, M. 2007, ApJ, 661, 1094 [NASA ADS] [CrossRef] [Google Scholar]
 Chugunov, A. I., Dewitt, H. E., & Yakovlev, D. G. 2007, Phys. Rev. D, 76, 025028 [NASA ADS] [CrossRef] [Google Scholar]
 Cox, J. P., & Giuli, R. T. 1968, Principles of Stellar Structure (New York: Gordon and Breach) [Google Scholar]
 Cyburt, R. H., Amthor, A. M., Ferguson, R., et al. 2010, ApJS, 189, 240 [NASA ADS] [CrossRef] [Google Scholar]
 de Mink, S. E., Sana, H., Langer, N., Izzard, R. G., & Schneider, F. R. N. 2014, ApJ, 782, 7 [Google Scholar]
 Eggen, O. J. 1967, Mem. R. Astron. Soc., 70, 111 [NASA ADS] [Google Scholar]
 Evans, C. J., Taylor, W. D., HénaultBrunet, V., et al. 2011, A&A, 530, A108 [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]
 Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585 [Google Scholar]
 Ferrario, L., Pringle, J. E., Tout, C. A., & Wickramasinghe, D. T. 2009, MNRAS, 400, L71 [NASA ADS] [CrossRef] [Google Scholar]
 Flannery, B. P. 1976, ApJ, 205, 217 [NASA ADS] [CrossRef] [Google Scholar]
 Fossati, L., Castro, N., Schöller, M., et al. 2015, A&A, 582, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Frost, A. J., Sana, H., Mahy, L., et al. 2023, Science, submitted [Google Scholar]
 Fuller, G. M., Fowler, W. A., & Newman, M. J. 1985, ApJ, 293, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Gräfener, G., Vink, J. S., de Koter, A., & Langer, N. 2011, A&A, 535, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Grunhut, J. H., Wade, G. A., Neiner, C., et al. 2017, MNRAS, 465, 2432 [NASA ADS] [CrossRef] [Google Scholar]
 Hamann, W.R., Koesterke, L., & Wessolowski, U. 1995, A&A, 299, 151 [NASA ADS] [Google Scholar]
 Hazlehurst, J. 1985, A&A, 145, 25 [NASA ADS] [Google Scholar]
 Hazlehurst, J. 1993, A&A, 271, 209 [NASA ADS] [Google Scholar]
 Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [Google Scholar]
 Iglesias, C. A., & Rogers, F. J. 1993, ApJ, 412, 752 [Google Scholar]
 Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
 Itoh, N., Hayashi, H., Nishikawa, A., & Kohyama, Y. 1996, ApJS, 102, 411 [NASA ADS] [CrossRef] [Google Scholar]
 Jermyn, A. S., Schwab, J., Bauer, E., Timmes, F. X., & Potekhin, A. Y. 2021, ApJ, 913, 72 [NASA ADS] [CrossRef] [Google Scholar]
 Jermyn, A. S., Bauer, E. B., Schwab, J., et al. 2023, ApJS, 265, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Kähler, H. 1989, A&A, 209, 67 [NASA ADS] [Google Scholar]
 Kähler, H. 2004, A&A, 414, 317 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kippenhahn, R., Ruschenplatt, G., & Thomas, H.C. 1980, A&A, 91, 175 [NASA ADS] [Google Scholar]
 Köhler, K., Langer, N., de Koter, A., et al. 2015, A&A, 573, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kuiper, G. P. 1941, ApJ, 93, 133 [NASA ADS] [CrossRef] [Google Scholar]
 Langanke, K., & MartínezPinedo, G. 2000, Nucl. Phys. A, 673, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Langer, N. 2022, ArXiv eprints [arXiv:2209.04165] [Google Scholar]
 Langer, N., Fricke, K. J., & Sugimoto, D. 1983, A&A, 126, 207 [NASA ADS] [Google Scholar]
 Lorenzo, J., Negueruela, I., Baker, A. K. F. V., et al. 2014, A&A, 572, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lorenzo, J., SimónDíaz, S., Negueruela, I., et al. 2017, A&A, 606, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lubow, S. H., & Shu, F. H. 1977, ApJ, 216, 517 [NASA ADS] [CrossRef] [Google Scholar]
 Lubow, S. H., & Shu, F. H. 1979, ApJ, 229, 657 [NASA ADS] [CrossRef] [Google Scholar]
 Lucy, L. B. 1968, ApJ, 151, 1123 [Google Scholar]
 Lucy, L. B. 1976, ApJ, 205, 208 [NASA ADS] [CrossRef] [Google Scholar]
 Mahy, L., Almeida, L. A., Sana, H., et al. 2020, A&A, 634, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mandel, I., & de Mink, S. E. 2016, MNRAS, 458, 2634 [NASA ADS] [CrossRef] [Google Scholar]
 Marchant, P. 2016, PhD Thesis, Rheinischen FriedrichWilhelmsUniversität Bonn [Google Scholar]
 Marchant, P., Langer, N., Podsiadlowski, P., Tauris, T. M., & Moriya, T. J. 2016, A&A, 588, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marchant, P., Pappas, K. M. W., GallegosGarcia, M., et al. 2021, A&A, 650, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Menon, A., Langer, N., de Mink, S. E., et al. 2021, MNRAS, 507, 5013 [NASA ADS] [CrossRef] [Google Scholar]
 Nieuwenhuijzen, H., & de Jager, C. 1995, A&A, 302, 811 [NASA ADS] [Google Scholar]
 Oda, T., Hino, M., Muto, K., Takahara, M., & Sato, K. 1994, At. Data Nucl. Data Tables, 56, 231 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
 Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [Google Scholar]
 Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
 Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10 [Google Scholar]
 Potekhin, A. Y., & Chabrier, G. 2010, Contrib. Plasma Phys., 50, 82 [NASA ADS] [CrossRef] [Google Scholar]
 Poutanen, J. 2017, ApJ, 835, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Rogers, F. J., & Nayfonov, A. 2002, ApJ, 576, 1064 [Google Scholar]
 Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444 [Google Scholar]
 Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713 [NASA ADS] [CrossRef] [Google Scholar]
 Schneider, F. R. N., Podsiadlowski, P., Langer, N., Castro, N., & Fossati, L. 2016, MNRAS, 457, 2355 [Google Scholar]
 Schneider, F. R. N., Ohlmann, S. T., Podsiadlowski, P., et al. 2019, Nature, 574, 211 [Google Scholar]
 Schootemeijer, A., Langer, N., Grin, N. J., & Wang, C. 2019, A&A, 625, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sen, K., Langer, N., Marchant, P., et al. 2022, A&A, 659, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Shu, F. H., Lubow, S. H., & Anderson, L. 1976, ApJ, 209, 536 [NASA ADS] [CrossRef] [Google Scholar]
 Shu, F. H., Lubow, S. H., & Anderson, L. 1979, ApJ, 229, 223 [NASA ADS] [CrossRef] [Google Scholar]
 Szymanski, M., Kubiak, M., & Udalski, A. 2001, Acta Astron., 51, 259 [NASA ADS] [Google Scholar]
 Tassoul, J.L. 2000, Stellar Rotation (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
 Timmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501 [NASA ADS] [CrossRef] [Google Scholar]
 Vanbeveren, D., De Donder, E., Van Bever, J., Van Rensbergen, W., & De Loore, C. 1998, New Astron., 3, 443 [CrossRef] [Google Scholar]
 Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 von Zeipel, H. 1924, MNRAS, 84, 665 [NASA ADS] [CrossRef] [Google Scholar]
 Wickramasinghe, D. T., Tout, C. A., & Ferrario, L. 2014, MNRAS, 437, 675 [NASA ADS] [CrossRef] [Google Scholar]
 Yang, Y., Yuan, H., & Dai, H. 2019, AJ, 157, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Yoon, S.C., & Langer, N. 2005, A&A, 443, 643 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yoon, S.C., Langer, N., & Norman, C. 2006, A&A, 460, 199 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Appendix A: Energy transfer location
Because our ET formalism is implemented explicitly and the radius of a stellar model can change during the solver iterations, which is also coupled to mass added or subtracted due to MT, a nonnegligible mismatch between the location of energy input and the RL can occur. In particular, when mass is removed due to MT, the ET location is then evaluated on the remaining cells, without those cells being adjusted since the previous step. Therefore, in phases with a high MT rate, it is possible that all mass above the RL is removed in one step (although after the solver performs the integration, the star will expand due to thermal relaxation). This means that the location of ET in terms of radius is rather poorly defined. To mitigate this, instead of choosing directly to put energy in cells with radii:
we opt to do following. At the beginning of the evolutionary step, we evaluate:
as the depth fraction of the current RL in mass coordinate. Then during the next integration, the location of energy input is in the cells with radii:
This treatment mitigates mismatches in ET locations introduced by MT during evolutionary steps to the degree that the depth fraction changes from step to step.
Appendix B: Convergence study
This appendix details the convergence study we used to validate our stellar model with and without energy transfer. We have implemented the ET method in both MESA version r15140, as well as version r22.11.1.
B.1. Explicit evaluation
In MESA r15140, both the tidal deformation corrections of Paper I and the ET of Sect. 3.4 are evaluated explicitly in the differential equation solver, meaning that only information on the previous evolutionary step was used.
We performed a time resolution study where we took a base template and increase the time resolution by successive factors of two. Figure B.1 shows the evolution of the MT rate, ET rate, the mass ratio, and the donor radius during the fast case A MT phase. While the time of the onset of MT, which coincides with the donor reaching the RL, remains roughly constant for the different resolutions, the MT rate evolves significantly differently. For higher time resolutions, the MT rate is higher, resulting in slightly shorter contact phase and a higher final mass ratio. This effect is independent of our implementation of ET since the changes occur before the contact phase (signified by the bump in the donor radius). Furthermore, these models do not yet converge to a certain value for the quantities of interest, which are the contact phase duration and the mass ratio after the MT phase. This means we would need to increase the time resolution further to get to a convergent solution. This is prohibitive however keeping in mind the aim is to run a full population of models in later work. We therefore elected to implement the tidal deformation corrections in MESA version 22.11.1 to allow for implicit evaluation during integrations.
Fig. B.1. Evolution of the MT rate and ET rate (top panel) and the fractional RL radius and mass ratio (bottom panel) during fast case A MT for a 20 + 14 M_{⊙} binary, computed using the explicit evaluation of the tidal deformation corrections in MESA r15140. The color indicates what time resolution scale has been used. 
B.2. Implicit evaluation
In MESA r22.11.1, infrastructure is provided to supply numerically approximated derivatives necessary for the solver to implicitly evaluate the tidal deformation corrections during an integration step. Generally, implicit evaluation is beneficial for numerical stability, and allows for longer time steps. We note that the amount of energy to be transferred and its location in the stellar structure is still evaluated explicitly, as detailed in Sect. 3.4.
We carried out a similar resolution study as the one described in Sect. B.1 of a 20 + 14 M_{⊙} binary, were we progressively increase the time resolution between simulations. Compared to version r15140, we observe a much better convergence of the MT rate evolution, as pictured in Fig. B.2.
Fig. B.2. Evolution of the MT rate and ET rate (top panel) and the fractional RL radius and mass ratio (bottom panel) during fast case A MT for a 20 + 14 M_{⊙} binary, but, as opposed to Fig. B.1, now the models were computed using implicit evaluation of the tidal deformation corrections in MESA r22.11.1. The color indicates what time resolution scale has been used. 
Finally, we performed a spacial resolution convergence study for the same initial conditions as before. Figure B.3 shows that our models are well converged in the spacial coordinate already at the 1x level. Even during the contact phase where ET occurs, and a sharp increase in luminosity is introduced, the models depend very weakly on the spacial resolution.
Fig. B.3. Evolution of the MT rate and ET rate (top panel) and the fractional RL radius and mass ratio (bottom panel) during fast case A MT for a 20 + 14 M_{⊙} binary, using implicit evaluation of the tidal deformation corrections in MESA r22.11.1. Here the color indicates the increase in spacial resolution. 
All Figures
Fig. 1. HRD of a 25 + 20 M_{⊙} binary with an initial period of 1.4 d. Red lines indicate a simulation with the inclusion of ET, while blue lines ignore ET. Moving left to right, the columns highlight four subsequent phases in the evolution (see Sect. 4.1 for details). Contact phases are highlighted with thick outlines. 

In the text 
Fig. 2. MT rate (top panel) and mass ratio (bottom panel) evolution of the 25 + 20 M_{⊙} binary studied in Sect. 4.1. Contact phases are highlighted with thick outlines. 

In the text 
Fig. 3. Contact duration distributions as function of observed mass ratio . Top panel: differential duration (in bins of ). Bottom panel: cumulative duration in configurations more extreme than . 

In the text 
Fig. 4. Profiles of outer layers of the binary components at the onset of inverse MT. It shows (top to bottom) the temperature, density, and luminosity profiles for both components in the ET (red, left column) and no ET (blue, right column) cases, all as function of the scaled radius coordinate . The gray vertical lines show the location of the RL. 

In the text 
Fig. 5. Thickness of the ET layer, d, with respect to the binary separation a as function of the radius of the primary star during the long lived contact phase. The gray line gives the physical size of the overflowing layers, and corresponds to the maximal width the ET layer can assume. 

In the text 
Fig. 6. Luminosity ratio versus mass ratio during the nuclear timescale MT of the 19 + 14 M_{⊙} systems of Sect. 4.1. Measurements from observed massive (near)contact systems from AbdulMasih et al. (2021), Yang et al. (2019), and Lorenzo et al. (2014) are overplotted. The systems from Mahy et al. (2020) were classified as ‘uncertain configurations’. 

In the text 
Fig. B.1. Evolution of the MT rate and ET rate (top panel) and the fractional RL radius and mass ratio (bottom panel) during fast case A MT for a 20 + 14 M_{⊙} binary, computed using the explicit evaluation of the tidal deformation corrections in MESA r15140. The color indicates what time resolution scale has been used. 

In the text 
Fig. B.2. Evolution of the MT rate and ET rate (top panel) and the fractional RL radius and mass ratio (bottom panel) during fast case A MT for a 20 + 14 M_{⊙} binary, but, as opposed to Fig. B.1, now the models were computed using implicit evaluation of the tidal deformation corrections in MESA r22.11.1. The color indicates what time resolution scale has been used. 

In the text 
Fig. B.3. Evolution of the MT rate and ET rate (top panel) and the fractional RL radius and mass ratio (bottom panel) during fast case A MT for a 20 + 14 M_{⊙} binary, using implicit evaluation of the tidal deformation corrections in MESA r22.11.1. Here the color indicates the increase in spacial resolution. 

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.