EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 672 (April 2023) A&A, 672 (2023) A175 Full HTML
Open Access
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/0004-6361/202346277
Published online 19 April 2023

© The Authors 2023

Licence Creative CommonsOpen 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 (Minitial ≳ 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 so-called neck and form a peanut-like 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 long-duration gamma-ray bursts and magnetar-driven superluminous supernovae (Yoon & Langer 2005; Yoon et al. 2006; Aguilera-Dena 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, low-mass 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 mass-radius 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 mass-radius 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 (L1). 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, long-lived 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 M2/M1 = 0.8 − 1 and the effect was thought to be minimal.

Sen et al. (2022) used the models of Marchant (2016) to study the semi-detached 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 M2/M1 = 0.6 with the study of massive, long-lived 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 mass-radius 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:

(1)

where we defined q ≡ M2/M1. However, the condition of contact in a binary following the Roche geometry constrains the surface of the stars to same equipotential, which leads to:

(2)

for stars not overflowing their RL to excess. Clearly, Eqs. (1) and (2) cannot be satisfied simultaneously unless M1 = M2. 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 mass-radius 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 mass-radius 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 (M2/M1)β, with β = 0.96. Using the approximation β ≃ 1, combined with von Zeipel’s theorem of gravity darkening, von Zeipel (1924), and the Stefan-Boltzmann law, , results in the simple expectation that in contact binaries, the luminosity ratio follows the mass ratio (Lucy 1968; Tassoul 2000):

(3)

Single main sequence stars, on the other hand, follow the well known mass-luminosity relation:

(4)

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 Ls, 1 and the luminosity L1 of a star of the same mass in a contact binary of

(5)

and for the companion

(6)

since L2 ≃ qL1 and we require ΔL1 + ΔL2 = 0 to conserve energy. This defines f as

(7)

Therefore, the two stars in a contact binary can fulfill the single star mass-luminosity relation of Eq. (4) in their cores and the contact binary mass-luminosity 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:

(8a)

(8b)

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

(9)

Comparing Eqs. (8a) and (8b) against Eqs. (5) and (6), we find for the fraction f:

(10)

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 1H, 3He, 4He, 12C, 14N, 16O, 20Ne, and 24Mg, 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ínez-Pinedo (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 CO-enhanced 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 bi-stability 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 Bondi-Hoyle 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öhm-Vitense (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 one-dimensional (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 self-consistently 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 L1, 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:

(11)

Given a certain splitting strategy, and the requirement of the stellar surface being on a common equipotential, a relation of the form:

(12)

with q = M2/M1, 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 L1, Marchant et al. (2016) found the following fit:

(13)

where we have denoted rRL = 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, Lout, which is L2 for the less massive and L3 for the more massive component, we constructed a new fit from the integration results of Paper I, which is:

(14a)

(14b)

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 L2, 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, horizontal1 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:

(15)

where is the smoothed transferred luminosity at an evolutionary step, n, Ltrans 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:

(16a)

(16b)

where Δm1, 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:

(17)

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

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 qini = 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 Abdul-Masih 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 pre-interaction 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 L1. 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 semi-detached 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 L2 overflow occurs.

thumbnail 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 vice-versa for the primary). The point at which the primary and secondary tracks touch (at roughly log L/L⊙ = 4.8 and log Teff/K = 4.545) corresponds to the situation where M1 = M2. This is consistent with the theory since we expect that for M1 = M2 we have similar luminosities per Eq. (3) and also similar effective temperatures, following the law of Stefan-Boltzmann. 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 non-energy-transferring models during the final, inverted MT phase, which is long-lived 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 no-ET 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.

thumbnail 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 long-lived 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.

thumbnail 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 M1 = 18.3 M⊙, while the secondary accreted mass to M2 = 26.4 M⊙. We note that these masses are close to the component masses Abdul-Masih 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).

thumbnail 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:

(18)

with a as the binary separation, and ρ, h, and cs 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.

thumbnail 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 L1 on the primary star (location i) toward L1 (location f), we have

(19)

where we have already canceled the potential terms Ψi, Ψf since we move along an equipotential surface, and the initial velocity, vi, is assumed to be negligible. Making the estimation:

(20)

and assuming that the transferred energy through the binary neck of width, b, by a mass flow of Ṁ = ρivfb2 is:

(21)

we compute for the minimal thickness of the ET layer:

(22)

Finally the thickness of the layer at the neck is related to the thickness far away from L1 via:

(23)

since the Roche potential varies quadratically near L1 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, Ltrans. Conversely, it can be interpreted as Ltrans 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 − RRL 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 mass-luminosity 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 q2.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 mass-luminosity ratio changes drastically from the q2.2 line in near contact to the q1 relation in full contact.

thumbnail 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 Abdul-Masih 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 Abdul-Masih et al. (2021), Yang et al. (2019) and Lorenzo et al. (2014) (see also Fig. 2 of Langer 2022). Curiously, the luminosity ratios from Abdul-Masih et al. (2021), while they do follow an L ∝ M trend, are offset to lower L2/L1 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 mass-luminosity 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.

1

Horizontal in the context of this work means perpendicular to the local effective gravity g = ∇Ψ. Conversely, vertical means parallel to ∇Ψ.

2

We use here since in an observed system, one has no a priori information on what the most massive component is.

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

  1. Abdul-Masih, M., Sana, H., Hawcroft, C., et al. 2021, A&A, 651, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Aguilera-Dena, D. R., Langer, N., Moriya, T. J., & Schootemeijer, A. 2018, ApJ, 858, 115 [NASA ADS] [CrossRef] [Google Scholar]
  3. Aguilera-Dena, D. R., Langer, N., Antoniadis, J., & Müller, B. 2020, ApJ, 901, 114 [NASA ADS] [CrossRef] [Google Scholar]
  4. Alcock, C., Allsman, R. A., Alves, D., et al. 1997, AJ, 114, 326 [NASA ADS] [CrossRef] [Google Scholar]
  5. Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
  6. Biermann, P., & Thomas, H. C. 1972, A&A, 16, 60 [NASA ADS] [Google Scholar]
  7. Binnendijk, L. 1970, Vistas Astron., 12, 217 [NASA ADS] [CrossRef] [Google Scholar]
  8. Blouin, S., Shaffer, N. R., Saumon, D., & Starrett, C. E. 2020, ApJ, 899, 46 [NASA ADS] [CrossRef] [Google Scholar]
  9. Böhm-Vitense, E. 1958, ZAp, 46, 108 [NASA ADS] [Google Scholar]
  10. Bondi, H., & Hoyle, F. 1944, MNRAS, 104, 273 [Google Scholar]
  11. Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  12. Cassisi, S., Potekhin, A. Y., Pietrinferni, A., Catelan, M., & Salaris, M. 2007, ApJ, 661, 1094 [NASA ADS] [CrossRef] [Google Scholar]
  13. Chugunov, A. I., Dewitt, H. E., & Yakovlev, D. G. 2007, Phys. Rev. D, 76, 025028 [NASA ADS] [CrossRef] [Google Scholar]
  14. Cox, J. P., & Giuli, R. T. 1968, Principles of Stellar Structure (New York: Gordon and Breach) [Google Scholar]
  15. Cyburt, R. H., Amthor, A. M., Ferguson, R., et al. 2010, ApJS, 189, 240 [NASA ADS] [CrossRef] [Google Scholar]
  16. de Mink, S. E., Sana, H., Langer, N., Izzard, R. G., & Schneider, F. R. N. 2014, ApJ, 782, 7 [Google Scholar]
  17. Eggen, O. J. 1967, Mem. R. Astron. Soc., 70, 111 [NASA ADS] [Google Scholar]
  18. Evans, C. J., Taylor, W. D., Hénault-Brunet, V., et al. 2011, A&A, 530, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  19. Fabry, M., Marchant, P., & Sana, H. 2022, A&A, 661, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  20. Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585 [Google Scholar]
  21. Ferrario, L., Pringle, J. E., Tout, C. A., & Wickramasinghe, D. T. 2009, MNRAS, 400, L71 [NASA ADS] [CrossRef] [Google Scholar]
  22. Flannery, B. P. 1976, ApJ, 205, 217 [NASA ADS] [CrossRef] [Google Scholar]
  23. Fossati, L., Castro, N., Schöller, M., et al. 2015, A&A, 582, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Frost, A. J., Sana, H., Mahy, L., et al. 2023, Science, submitted [Google Scholar]
  25. Fuller, G. M., Fowler, W. A., & Newman, M. J. 1985, ApJ, 293, 1 [NASA ADS] [CrossRef] [Google Scholar]
  26. Gräfener, G., Vink, J. S., de Koter, A., & Langer, N. 2011, A&A, 535, A56 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Grunhut, J. H., Wade, G. A., Neiner, C., et al. 2017, MNRAS, 465, 2432 [NASA ADS] [CrossRef] [Google Scholar]
  28. Hamann, W.-R., Koesterke, L., & Wessolowski, U. 1995, A&A, 299, 151 [NASA ADS] [Google Scholar]
  29. Hazlehurst, J. 1985, A&A, 145, 25 [NASA ADS] [Google Scholar]
  30. Hazlehurst, J. 1993, A&A, 271, 209 [NASA ADS] [Google Scholar]
  31. Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [Google Scholar]
  32. Iglesias, C. A., & Rogers, F. J. 1993, ApJ, 412, 752 [Google Scholar]
  33. Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
  34. Itoh, N., Hayashi, H., Nishikawa, A., & Kohyama, Y. 1996, ApJS, 102, 411 [NASA ADS] [CrossRef] [Google Scholar]
  35. Jermyn, A. S., Schwab, J., Bauer, E., Timmes, F. X., & Potekhin, A. Y. 2021, ApJ, 913, 72 [NASA ADS] [CrossRef] [Google Scholar]
  36. Jermyn, A. S., Bauer, E. B., Schwab, J., et al. 2023, ApJS, 265, 15 [NASA ADS] [CrossRef] [Google Scholar]
  37. Kähler, H. 1989, A&A, 209, 67 [NASA ADS] [Google Scholar]
  38. Kähler, H. 2004, A&A, 414, 317 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Kippenhahn, R., Ruschenplatt, G., & Thomas, H.-C. 1980, A&A, 91, 175 [NASA ADS] [Google Scholar]
  40. Köhler, K., Langer, N., de Koter, A., et al. 2015, A&A, 573, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Kuiper, G. P. 1941, ApJ, 93, 133 [NASA ADS] [CrossRef] [Google Scholar]
  42. Langanke, K., & Martínez-Pinedo, G. 2000, Nucl. Phys. A, 673, 481 [NASA ADS] [CrossRef] [Google Scholar]
  43. Langer, N. 2022, ArXiv e-prints [arXiv:2209.04165] [Google Scholar]
  44. Langer, N., Fricke, K. J., & Sugimoto, D. 1983, A&A, 126, 207 [NASA ADS] [Google Scholar]
  45. Lorenzo, J., Negueruela, I., Baker, A. K. F. V., et al. 2014, A&A, 572, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Lorenzo, J., Simón-Díaz, S., Negueruela, I., et al. 2017, A&A, 606, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Lubow, S. H., & Shu, F. H. 1977, ApJ, 216, 517 [NASA ADS] [CrossRef] [Google Scholar]
  48. Lubow, S. H., & Shu, F. H. 1979, ApJ, 229, 657 [NASA ADS] [CrossRef] [Google Scholar]
  49. Lucy, L. B. 1968, ApJ, 151, 1123 [Google Scholar]
  50. Lucy, L. B. 1976, ApJ, 205, 208 [NASA ADS] [CrossRef] [Google Scholar]
  51. Mahy, L., Almeida, L. A., Sana, H., et al. 2020, A&A, 634, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Mandel, I., & de Mink, S. E. 2016, MNRAS, 458, 2634 [NASA ADS] [CrossRef] [Google Scholar]
  53. Marchant, P. 2016, PhD Thesis, Rheinischen Friedrich-Wilhelms-Universität Bonn [Google Scholar]
  54. Marchant, P., Langer, N., Podsiadlowski, P., Tauris, T. M., & Moriya, T. J. 2016, A&A, 588, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Marchant, P., Pappas, K. M. W., Gallegos-Garcia, M., et al. 2021, A&A, 650, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Menon, A., Langer, N., de Mink, S. E., et al. 2021, MNRAS, 507, 5013 [NASA ADS] [CrossRef] [Google Scholar]
  57. Nieuwenhuijzen, H., & de Jager, C. 1995, A&A, 302, 811 [NASA ADS] [Google Scholar]
  58. Oda, T., Hino, M., Muto, K., Takahara, M., & Sato, K. 1994, At. Data Nucl. Data Tables, 56, 231 [NASA ADS] [CrossRef] [Google Scholar]
  59. Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
  60. Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [Google Scholar]
  61. Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
  62. Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34 [NASA ADS] [CrossRef] [Google Scholar]
  63. Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10 [Google Scholar]
  64. Potekhin, A. Y., & Chabrier, G. 2010, Contrib. Plasma Phys., 50, 82 [NASA ADS] [CrossRef] [Google Scholar]
  65. Poutanen, J. 2017, ApJ, 835, 119 [NASA ADS] [CrossRef] [Google Scholar]
  66. Rogers, F. J., & Nayfonov, A. 2002, ApJ, 576, 1064 [Google Scholar]
  67. Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444 [Google Scholar]
  68. Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713 [NASA ADS] [CrossRef] [Google Scholar]
  69. Schneider, F. R. N., Podsiadlowski, P., Langer, N., Castro, N., & Fossati, L. 2016, MNRAS, 457, 2355 [Google Scholar]
  70. Schneider, F. R. N., Ohlmann, S. T., Podsiadlowski, P., et al. 2019, Nature, 574, 211 [Google Scholar]
  71. Schootemeijer, A., Langer, N., Grin, N. J., & Wang, C. 2019, A&A, 625, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Sen, K., Langer, N., Marchant, P., et al. 2022, A&A, 659, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  73. Shu, F. H., Lubow, S. H., & Anderson, L. 1976, ApJ, 209, 536 [NASA ADS] [CrossRef] [Google Scholar]
  74. Shu, F. H., Lubow, S. H., & Anderson, L. 1979, ApJ, 229, 223 [NASA ADS] [CrossRef] [Google Scholar]
  75. Szymanski, M., Kubiak, M., & Udalski, A. 2001, Acta Astron., 51, 259 [NASA ADS] [Google Scholar]
  76. Tassoul, J.-L. 2000, Stellar Rotation (Cambridge: Cambridge University Press) [CrossRef] [Google Scholar]
  77. Timmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501 [NASA ADS] [CrossRef] [Google Scholar]
  78. Vanbeveren, D., De Donder, E., Van Bever, J., Van Rensbergen, W., & De Loore, C. 1998, New Astron., 3, 443 [CrossRef] [Google Scholar]
  79. Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  80. von Zeipel, H. 1924, MNRAS, 84, 665 [NASA ADS] [CrossRef] [Google Scholar]
  81. Wickramasinghe, D. T., Tout, C. A., & Ferrario, L. 2014, MNRAS, 437, 675 [NASA ADS] [CrossRef] [Google Scholar]
  82. Yang, Y., Yuan, H., & Dai, H. 2019, AJ, 157, 111 [NASA ADS] [CrossRef] [Google Scholar]
  83. Yoon, S.-C., & Langer, N. 2005, A&A, 443, 643 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. 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 non-negligible 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:

(A.1)

we opt to do following. At the beginning of the evolutionary step, we evaluate:

(A.2)

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:

(A.3)

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.

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

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

thumbnail 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

thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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
thumbnail 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 Abdul-Masih 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
thumbnail 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
thumbnail 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
thumbnail 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 (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.