EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 695 (March 2025) A&A, 695 (2025) A223 Full HTML
Open Access
Issue
A&A
Volume 695, March 2025
Article Number A223
Number of page(s) 12
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/202453345
Published online 20 March 2025

© The Authors 2025

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

Observations indicate that the majority of massive stars (M  ≳  8 M⊙) are in binary systems (Sana et al. 2014; Villaseñor et al. 2021; Banyard et al. 2022). Sana et al. (2012) show that half of all massive stars interact before the end of the main sequence; this is called Case A evolution (Kippenhahn et al. 1967; Pols 1994; Vanbeveren et al. 1998; Nelson & Eggleton 2001).

Case A systems have initial orbital periods of around ten days or less. When such a massive binary evolves, an extreme form of interaction can occur when both stars overflow their Roche lobes simultaneously. This is the contact phase, and the system is then classified as a contact binary1. Previous studies show that a contact phase is expected for at least 40% of all massive binaries with a mass-transfer event (Henneco et al. 2024).

Massive contact binaries are likely to merge (Menon et al. 2021; Henneco et al. 2024; Fabry et al. 2025), which has been identified as a possible process for creating magnetic massive stars (Schneider et al. 2016, 2019). Recently, this scenario has been observationally confirmed in HD 148937 (Frost et al. 2024). Mergers are also thought to create rapidly rotating neutron stars, long gamma-ray bursts (Yoon & Langer 2005; Yoon et al. 2006), and super-luminous supernovae (Aguilera-Dena et al. 2018, 2020). At low metallicity, massive contact binaries are expected to make a dominating contribution to the binary black hole population that is observed through gravitational waves (Mandel & de Mink 2016; Marchant et al. 2016; Mandel & Farmer 2022).

Despite many efforts to find massive contact binaries, only around two dozen are known. There are several systems in the Galactic environment (e.g. Hilditch & Evans 1985; Lorenzo et al. 2017; Abdul-Masih et al. 2021; Li et al. 2022). Additionally, large surveys including the Massive Compact Halo Objects (MACHO) Project (Alcock et al. 1997) and the Tarantula Massive Binary Monitoring (Almeida et al. 2015, 2017; Mahy et al. 2020a,b) detected such systems. These systems are well studied. For example, recent studies investigated the orbital variations of a sample of O- and B-type massive contact binaries and studied the evolutionary timescales of these systems (Abdul-Masih et al. 2022; Vrancken et al. 2024). A more extensive overview of observed massive contact binaries can be found in Appendix B.

The theoretical description of contact systems was initiated by Kuiper (1941), who identified that contact binaries cannot be stable if the components obey a single-star internal structure, except at a mass ratio of unity. Contact binaries are preferentially observed at mass ratios away from unity, however (especially in the low-mass W UMa regime), which led to Kuiper’s argument being called Kuiper’s paradox. The solution to the paradox is to consider stellar models that are affected by binary interaction, for example through mass and energy transfer, as well as tidal distortion.

Lucy (1968) was the first to consider energy transfer within the contact binary as a possible solution. Other models followed these works (e.g. Biermann & Thomas 1972; Flannery 1976; Shu et al. 1976; Webbink 1977). However, some were found to contain inconsistencies (Papaloizou & Pringle 1979; Hazlehurst & Refsdal 1980), or focused only on those low-mass contact binaries with convective envelopes.

In addition to analytical models, computational studies of contact binaries can be performed. This was done with 1-dimensional binary evolution codes for example by Pols (1994), Wellstein et al. (2001), Nelson & Eggleton (2001), de Mink & Pols (2007), and Marchant et al. (2016) for massive stars. More recently, the population synthesis calculation by Menon et al. (2021) showed an overestimation of equal mass ratio using current mass transfer models. They found that unequal mass ratio systems equalize on the nuclear timescale, but the models still show that the time spent as equal mass contact binaries dominates the population. However, Fabry et al. (2023) shows that, under certain conditions, the inclusion of energy transfer can delay the equalization of the mass ratio. The models also account for the tidal deformation of the components (Fabry et al. 2022). Despite these improvements, Fabry et al. (2025) shows that this only has a limited effect on the overall properties of the contact binary population. Therefore, other processes in stellar and binary physics should be explored to tackle this discrepancy.

Mass transfer significantly influences the evolution of the stars in a binary system. As a main sequence accretor grows during mass transfer, its core grows in parallel and rejuvenation of the accretor can take place. Due to the strong mixing induced by convection in massive stars, additional hydrogen gets mixed in and the hydrogen mass fraction (X) increases. The star will look younger than predicted by single stellar evolution (Dray & Tout 2007). Schneider et al. (2016) found that rejuvenation can make a merger product appear younger by a considerable fraction of its nuclear timescale. The strength of the rejuvenation increases for more evolved binaries, for lower mass binaries, and for binaries with a larger mass ratio. It is a key factor in the evolution of the binary system, and limiting the rejuvenation is expected to have large effects. One way to inhibit the rejuvenation is by limiting the mixing of additional fuel into the core through convective overshooting. The efficiency of rejuvenation is expected to have a large impact on the evolution of a close binary. The size and composition of the core is related to the radius of a star, which determines the mass transfer events in the system when overflowing the Roche lobe.

In evolutionary models of massive stars, convective boundary mixing (CBM) is still a large source of uncertainty (Johnston et al. 2024). One type of CBM is overshooting. However, both the extent and efficiency of the mixing in this region are still debated. In the last decade, using high-cadence photometry from space missions, asteroseismic analyses have provided accurate constraints on the convective core masses and overshooting parameters (Moravveji et al. 2015, Angelou et al. 2020, Viani & Basu 2020, Michielsen et al. 2021, Pedersen et al. 2021, Burssens et al. 2023).

Convective regions in a stellar model are determined by the Schwarzschild (Schwarzschild & Härm 1958) or Ledoux (Ledoux 1947) criterion. The Schwarzschild criterion compares the adiabatic gradient with the radiative gradient to determine stability against convection. The Ledoux criterion takes the effect of the molecular weight gradient into account, together with the adiabatic gradient. A molecular weight gradient is built up in the near core region as a massive star ages and its core recedes. During mass transfer events, the core needs to adjust, and it is possible to affect the size of the core post-mass transfer with an overshooting prescription dependent on the molecular weight gradient. Förster (2023) reported promising results, showing that when modifying the CBM prescriptions in massive contact binaries, it is possible to alter their mass ratio evolution. Their work was limited to a set of models of a single initial primary mass, initial mass ratio and initial period. In this work, we aim to study the effect of rejuvenation on a wide parameter range of contact binary progenitors.

This paper is organized as follows. Section 2 describes the convective core overshoot and implemented limit based on the molecular gradient. Section 3 describes the MESA implementation and the models. Section 4 contains the results. We describe one model in detail and discuss the results over the whole grid. Section 5 discusses the results and investigates the numerical behaviour of the grid. We summarize our conclusions in Sect. 6.

2. Convective mixing

Convection is a mode of heat transport that occurs when material is unstable under adiabatic displacement. It occurs in the stellar interior when radiative heat transport becomes inefficient. It is a very efficient mixer that leads to chemically homogeneous regions inside stars, for example in the cores of massive stars on the main sequence. The size of the convective region determines the total amount of hydrogen that is available inside the core. This has an immediate effect on the lifetime and evolution of massive stars.

When using the Ledoux criterion (Ledoux 1947), the stability condition takes the chemical gradient ∇ μ ≡ d ln μ d ln P $ {\boldsymbol{\nabla}}_\mu \equiv \frac{\mathrm{d}\ln \mu}{\mathrm{d}\ln P} $ into account. A region is stable to convection if

∇ rad ≤ ∇ ad + ϕ δ ∇ μ , $$ \begin{aligned} {\boldsymbol{\nabla }}_{\mathrm{rad}} \le {\boldsymbol{\nabla }}_{\mathrm{ad}} + \frac{\phi }{\delta }{\boldsymbol{\nabla }}_{\mu }, \end{aligned} $$(1)

where ϕ ≡ ( ∂ ln ρ ∂ ln μ ) P , T $ \phi \equiv \left( \frac{\partial\ln \rho}{\partial\ln \mu}\right)_{P, T} $, δ ≡ − ( ∂ ln ρ ∂ ln μ ) T , μ $ \delta \equiv - \left( \frac{\partial\ln \rho}{\partial\ln \mu} \right)_{T, \mu} $, ∇ ad ≡ ( d ln T d ln P ) ad $ {\boldsymbol{\nabla}}_{\mathrm{ad}} \equiv \left( \frac{\mathrm{d}\ln T}{\mathrm{d}\ln P}\right)_{\mathrm{ad}} $ is the temperature gradient under adiabatic conditions, and ∇ rad ≡ ( d ln T d ln P ) rad $ {\boldsymbol{\nabla}}_{\mathrm{rad}} \equiv \left( \frac{\mathrm{d}\ln T}{\mathrm{d}\ln P} \right)_{\mathrm{rad}} $ is the radiative temperature gradient.

A common prescription of convection for stellar models is local mixing length theory (Vitense 1953; Böhm-Vitense 1958; Bohm-Vitense & Nelson 1976). It calculates the efficiency of convection and the diffusion coefficient. The efficiency of convection is a measure of the ratio of the thermal adjustment timescale to the lifetime of the convective bubble. The diffusion coefficient quantifies the rate at which the convection makes the density, composition, and temperature gradient homogeneous. However, as this is a local theory, non-local effects such as convective core overshoot must be added ad hoc.

As fluid parcels rise in the convective zone of a star, they gain momentum. This non-zero momentum at the boundary keeps them from immediately dissolving, and they overshoot into the radiative zone. This leads to the additional mixing of envelope material into the core of the star. As such, the strength of the overshoot can have a large influence on the composition and evolution of the core.

Two common ways to implement the convective overshoot are the step- and exponential overshoot. The former keeps the convective diffusion coefficient at the inside of the convective boundary constant into the radiative zone for some prescribed length, while the latter has the diffusion coefficient exponentially decay beyond the boundary. Brott et al. (2011) calibrated the step overshoot length (lP) as a fraction (αov) of the pressure scale height (HP):

l P = α ov H P . $$ \begin{aligned} l_{P} = \alpha _{\mathrm{ov}} H_{\mathrm{P}}. \end{aligned} $$(2)

Using data from the Survey of Massive Stars (Evans et al. 2005), they found the best fitting result to be αov = 0.335.

The disadvantage of using the above formulation is that the overshooting region is entirely defined by the conditions at the convective boundary. It leaves no freedom for potentially changing conditions in the overshooting region. However, since the Ledoux criterion of convection takes into account molecular weight gradients, which can inhibit convective motion, we consider the possibility for the overshoot to be similarly limited by molecular weight gradients.

We define an upper limit (B) of the molecular weight gradient, above which we do not allow the overshoot to penetrate. The length of the overshoot region is thus given by solving

∇ μ ( l B ) = B , $$ \begin{aligned} {\boldsymbol{\nabla }}_\mu (l_{\mathrm{B}}) = B , \end{aligned} $$(3)

where lB is the gradient dependent overshoot length and B is a free parameter. To make sure we do not extend the overshooting region, we take the overshoot length as the minimum of lP and lB:

l ov = min { l P , l B } . $$ \begin{aligned} l_{\mathrm{ov}} = \min \{l_{P}, l_{\mathrm{B}} \}. \end{aligned} $$(4)

In the case where Eq. (3) has no solution, that is when ∇μ < B everywhere beyond the convective boundary, we revert to using lP for the overshooting length. We expect the evolution of our stellar models to be sensitive to the value of B chosen. A low value of B will lead to a small lB and as such will inhibit overshooting. The limit on the overshoot will inhibit the growth of the core during mass accretion, which will stall the rejuvenation of the core. All else being equal in a mass transfer event, the core of the accretor will be more evolved and the stellar radius will be larger.

3. Methods

We computed a grid of binary evolution models using Modules for Experiments in Stellar Astrophysics (MESA Paxton et al. 2011, 2013, 2015, 2018, 2019; Jermyn et al. 2023) version r23.05.1 with MesaSDK version 23.7.3 (Townsend 2023). Section 3.1 summarizes the microphysics, wind prescription, and binary interaction physics. Section 3.2 describes the mixing parameters used, while Sect. 3.3 gives the grid parameters and initial and termination conditions. Finally, Sect. 3.4 details the population synthesis calculations.

3.1. Microphysics, winds, rotation, and binary interaction

We used a nuclear net containing eight isotopes: 1H, 3He, 4He, 12C, 14N, 16O, 20Ne, and 24Mg. The appropriate nuclear reaction rates were from Cyburt et al. (2010) and Angulo et al. (1999). The tabulated weak reaction rates were from Fuller et al. (1985), Oda et al. (1994), and Langanke & Martínez-Pinedo (2000). Plasma screening was included via the prescription of Chugunov et al. (2007). Thermal neutrino loss rates were from Itoh et al. (1996). The MESA equation of state that we used was a blend of Rogers & Nayfonov (2002), Saumon et al. (1995), Timmes & Swesty (2000), Potekhin & Chabrier (2010), and Jermyn et al. (2021). Iglesias & Rogers (1993, 1996) provided the radiative opacities, with low-temperature data from Ferguson et al. (2005). The Compton-scattering dominated regime with high temperatures was from Poutanen (2017). Cassisi et al. (2007) and Blouin et al. (2020) gave the electron conduction opacities.

Mass loss through winds was dependent on the surface hydrogen fraction (X), following Brott et al. (2011). The prescription of Hamann et al. (1995) was applied when X < 0.4. For X > 0.7 and temperatures above the iron bi-stability jump (as calibrated by Vink et al. 2001), the mass loss through winds was calculated following Vink et al. (2001). For temperatures below the iron bi-stability jump and high hydrogen surface abundance, we took the maximum of the rates described by Vink et al. (2001) and Nieuwenhuijzen & de Jager (1995). For surface abundances between these regimes, we linearly interpolated the above cases. The wind was allowed to be accreted by the companion following the description of Hurley et al. (2002). All models assumed solar metallicity Z⊙ = 0.0142 and metal fractions from Asplund et al. (2009). The rotation rate of both stars in the models was uniformly synchronized on the timescale of an orbital period, as we expected short-period binaries to be rapidly synchronized (Zahn 1975).

We implemented conservative mass transfer. For a semi-detached configuration, the donor was kept just below the Roche lobe and the mass was transferred onto the accretor, while in a contact phase the surface of both stars was kept on the same equipotential surface. We took into account the effect of tidal deformation on the internal structure of the star following the corrections of Fabry et al. (2022), and implemented efficient energy transfer in contact phases as in Fabry et al. (2023).

3.2. Mixing parameters

The convective regions were determined using the Ledoux criterion. We used the mixing length theory of Böhm-Vitense (1958), as described by Cox & Giuli (1968), and set the mixing length parameter as α = 2. We allowed for convective boundary mixing through step overshooting, where the diffusion coefficient at 0.01 pressure scale heights within the convective zone is kept constant out to lov out of the convective zone. We limited the overshoot zone to lov = min{lP, lB} as described in Sect. 2. We followed the model of Langer et al. (1983) to implement semi-convection, with a high efficiency of αsc = 100 following Schootemeijer et al. (2019). Thermohaline mixing followed the scheme of Kippenhahn et al. (1980), with an efficiency parameter of αth = 1. We included rotational mixing by including dynamical shear instability, secular shear instability, the Goldreich-Schubert-Fricke instability, and Eddington-Sweet circulation (Heger et al. 2000). The molecular gradient was calculated following Paxton et al. (2013). To numerically stabilize derivatives of the composition, which may exhibit large steps over cell boundaries, a smoothing by two adjacent cells in both directions was used.

3.3. Initialization and termination of models

The grid consists of 4896 models. We considered massive binaries with an initial primary mass M1, init = 8 M⊙ − 40 M⊙ with a step ΔM = 2 M⊙. The mass ratio is the ratio of the secondary to the primary mass and takes initial values q init = M 2 , init M 1 , init = 0.40 − 0.95 $ q_{\mathrm{init}} =\frac{M_{\mathrm{2, init}}}{M_{\mathrm{1, init}}} = 0.40 - 0.95 $ with a step Δq = 0.05. We considered initial periods between 0.4 d and 4.159 d with a step Δlog10P = 0.04. The overshoot limit defined above in Eq. (3) was taken to be B = 10−4.

The models were spun up to be synchronous with the orbital period, afterwards they relaxed onto the zero-age main sequence (ZAMS), which is defined as the moment when the luminosity is within 1% of nuclear luminosity. The eccentricity was kept at zero, as is expected for close binaries (Zahn 1975).

There were multiple terminations possible for the simulations:

  1. L2 overflow: The envelope of the contact binary overflows the L2 Lagrange point, in which case the system is expected to quickly evolve to a merger.

  2. Convergence issues: The numerical solver cannot find an acceptable solution for a finite time step.

  3. High mass transfer rate: If the mass transfer rate exceeds 10−2 M⊙yr−1, the simulation tends toward dynamical timescale mass transfer, which we do not attempt to resolve.

  4. Survival: When one of the companions leaves the main sequence, we terminate the simulation.

If the model overflowed L2 before the ZAMS condition was met or Roche-lobe overflow at ZAMS occured, we discarded it from the population. When any of the termination conditions was met, the simulation was quitted. We made no distinction between the terminations and performed our further analysis on all models in the grid, regardless of termination.

3.4. Population synthesis

We compare the results of the performed simulation to observations of candidate massive contact binaries (see Appendix B for an overview). This was done by means of a cumulative distribution function (CDF).

We applied the observational mass ratio q obs = min ( M 2 M 1 , M 1 M 2 ) $ q_{\mathrm{obs}} = \min \left( \frac{M_2}{M_1}, \frac{M_1}{M_2}\right) $ to compare to the observational values. We divided the observational mass ratio and the period into bins and considered the time each model spent in each bin. We obtained a duration of the form t(qobs, pobs, Mi, Pi, qi), where qobs denotes the observational mass ratio bin and pobs denotes the period bin. Symbols Mi, Pi, and qi denote the initial primary mass, initial period, and initial mass ratio, respectively. To construct a synthetic population from our grid of models, we needed to make assumptions about the distributions of initial masses, periods, and mass ratios.

To weigh the different initial periods (Pi), we applied Öpik’s law (Öpik 1924). This law states that the distribution is flat in logarithmic space, which is supported by multiple surveys of massive binaries (Dunstall et al. 2015; Almeida et al. 2017; Villaseñor et al. 2021; Banyard et al. 2022). This means that the initial periods are weighted as

w P i ( P i ) d P i = d log 10 P i = 1 P i ln ( 10 ) d P i , $$ \begin{aligned} { w}_{P_i}(P_i) \mathrm{d} P_i = \mathrm{d} \log _{10} P_i = \frac{1}{P_i \ln (10)} \mathrm{d} P_i, \end{aligned} $$(5)

where wPi represents the weight for the initial periods, Pi.

To correctly weigh the different initial masses (Mi), we applied the Salpeter initial-mass function (IMF) (Salpeter 1955; Kroupa 2001; Bastian et al. 2010)

w M i ( M i ) = M i − 2.35 d M i , $$ \begin{aligned} { w}_{M_i} (M_i) = M_i^{-2.35} \mathrm{d} M_i, \end{aligned} $$(6)

where wMi represents the weight for the initial primary mass, Mi.

We assumed the initial mass ratio weight to be uniform:

w q i = q i 0 d q i = d q i . $$ \begin{aligned} { w}_{q_i} = q_i^0 \mathrm{d} q_i = \mathrm{d} q_i. \end{aligned} $$(7)

As we assumed the star-formation rate (SFR) to also be uniformly distributed, we can write

t ( q obs , p obs ) ∝ ∫ t ( q obs , p obs , M i , P i , q i ) w P i w M i w q i d P i d M i d q i = ∫ t ( q obs , p obs , M i , P i , q i ) M i − 2.35 P i ln ( 10 ) d P i d M i d q i $$ \begin{aligned} t(q_{\mathrm{obs}}, p_{\mathrm{obs}})&\propto \int t(q_{\mathrm{obs}}, p_{\mathrm{obs}}, M_i, P_i, q_i) { w}_{P_i} { w}_{M_i} { w}_{q_i} \mathrm{d} P_i \mathrm{d} M_i \mathrm{d} q_i \\&= \int t(q_{\mathrm{obs}}, p_{\mathrm{obs}}, M_i, P_i, q_i) \frac{M_i^{-2.35}}{P_i \ln (10)}\mathrm{d} P_i \mathrm{d} M_i \mathrm{d} q_i \end{aligned} $$

as we integrate over the initial primary masses, periods, and mass ratios.

The models take discrete values for the initial values as described above. We therefore needed to discretise the integral. We took each model to represent a finite area of the initial distribution. Let us consider the initial period distribution. It is uniformly distributed in logarithmic space. This gives

W P i = ∫ lower upper d log 10 P i = log 10 P i upper − log 10 P i lower . $$ \begin{aligned} W_{P_i} = \int _{\mathrm{lower}}^{\mathrm{upper}} \mathrm{d} \log _{10} P_i = \log _{10} P_i^{\mathrm{upper}} - \log _{10} P_i^{\mathrm{lower}}. \end{aligned} $$(8)

As we have taken our initial period distribution to be evenly spread in logarithmic space, this will always take the constant value of WPi = 0.04. For the initial mass distribution, the bins at the outer ends of the model grid are chosen as such to conserve the constant bin size.

By applying a normalization, we can consider t(qobs, pobs) to be probabilities 𝒫(qobs, pobs). This allows us to compute the CDF

CDF ( q , p ) = ∫ q 0 = 0.4 , p 0 = 0.4 d q , p P ( q obs , p obs ) d q obs d p obs . $$ \begin{aligned} \text{ CDF} (q, p) = \int _{q_0 = 0.4, p_0 = 0.4 \mathrm{d}}^{q, p} \mathcal{P} (q_{\mathrm{obs}}, p_{\mathrm{obs}}) \mathrm{d} q_{\mathrm{obs}} \mathrm{d} p_{\mathrm{obs}}. \end{aligned} $$(9)

The marginal distributions can be obtained by integrating over all observed periods or mass ratios.

4. Results

In this section we analyse and discuss our evolution models. We start by discussing a single evolutionary model in detail, with a focus on the mass transfer rates and the rejuvenation. Finally, we compare our synthetic population against observed massive contact binaries (listed in Table B.1).

4.1. Detailed model evolution

Here, we describe in detail an example of our case A interacting models. We selected the initial conditions of M1, initial = 24 M⊙, Pinitial = 1.377 d, and qinitial = 0.8. As for all the models, we set the overshoot limit at B = 10−4. To set a baseline model, we selected the MESA models of Fabry (2024) and Fabry et al. (2025). They are set up in a similar way, with the exception that they follow the overshooting lP = 0.335HP as defined in Eq. (2). The evolution of these models closely follows those of Menon et al. (2021). We refer to this model as the Full Overshoot (FO) model.

Figure 1 shows the mass ratio evolution of the model of this work together with the mass ratio evolution of the FO model. The FO model closely follows System 2 of Menon et al. (2021) as expected, with a brief contact phase (fast case A mass transfer) around 1 Myr followed by a short detached phase. Afterwards, a semi-detached phase follows as mass is transferred from the primary to the secondary (slow case A). A new contact phase is formed when the secondary reaches its Roche lobe and the system equalizes over a nuclear timescale. The stars merge around 3.9 Myr after the ZAMS.

thumbnail Fig. 1.

Example of the mass ratio evolution with initial conditions of M1, initial = 24 M⊙, Pinitial = 1.377 d, and qinitial = 0.8. Top panel: Mass ratio evolution of the system. The coloured lines show the new model, the grey line shows the FO mass ratio evolution. The dashed lines show the stellar radii normalized by the respective Roche lobe radii. Bottom panel: Mass transfer rate of the system of the new model (blue) and the FO model (grey). The dashed black line denotes the initial thermal mass transfer rate of the secondary.

The new overshoot limited model deviates considerably from the FO model. This is immediately noticeable on Fig. 1. We go over some important differences.

The new model has a longer lifetime by about 0.6 Myr. Additionally, the fast case A mass transfer starts about 40 kyr earlier than the reference model. As we limit the overshoot, we limit the mixing at the core boundary. As less envelope hydrogen is mixed into the core, the stars will age faster. The faster ageing is accompanied by expansion, and the donor reaches the Roche lobe faster.

The fast case A leads to a larger mass ratio before detachment. While the FO model has q = 1.19 after the fast case A, the new model is at mass ratio q = 1.25. The slow case A that follows pushes the mass ratio further from unity and takes about 1.1 Myr. This slow case A is notably not smooth. This seems to be a consequence of the system repeatedly going in and out of a contact phase. We suspect this behaviour to be of numerical origin. The numerical behaviour of our models is discussed in Sect. 5.

After the slow case A mass transfer event, multiple rapid mass ratio inversions take place over multiple contact phases and the mass ratio is pushed to more extreme values. A final contact phase brings the mass ratio down, but due to L2 overflow, the system is considered to be merged.

The new model has a similar early evolution to the FO models, but deviates strongly from about 1.5 Myr onwards. The slow convergence towards unity is replaced by multiple rapid mass transfer events. Each transfer seems to deviate the mass ratio further from unity. The evolution of the system in the HR diagram is shown in Appendix A.

We see that for the fast case A and following slow case A mass transfer events up to around 1.5 Myr, both models seem to show similar mass transfer rates. This is confirmed by the bottom panel of Fig. 1, which shows the mass transfer rate for the new model in blue and the mass transfer for the FO model in grey. Similar to above, the plot shows staggered behaviour as a consequence of the numerical implementation. We do see, however, that overall the mass transfer rate is higher in the new model especially at later times in the evolution, which we suspect to be a consequence of the larger radii that accompany the more evolved cores. The three mass ratio inversion events are accompanied by strong mass transfer rates, which approach the initial thermal mass transfer rate of the secondary log M ˙ KH , init , 2 = M init , 2 τ KH , 2 $ \log \dot{M}_{\mathrm{KH, init, 2}} = \frac{M_{\mathrm{init, 2}}}{\tau_{\mathrm{KH, 2}}} $ (black line), where τKH is the Kelvin-Helmholtz timescale.

4.2. Rejuvenation

The new method inhibits convective boundary mixing when the molecular gradient becomes larger than our boundary value. Less near-core mixing will lead to a smaller rejuvenation of the core. This can be seen in Fig. 2, which shows the evolution of the helium mass fraction of the convective cores (Ycc) of both the primary and the secondary (for the same model presented in Fig. 1).

thumbnail Fig. 2.

Evolution of the helium mass fraction in the core of the star. The upper plot shows the primary star, the lower plot shows the secondary star. The grey line denotes the FO model.

The effect of our model is clear before the first contact phase is initiated at 1 Myr, as the helium mass fraction rises more rapidly than in the FO case. This behaviour continues after the first fast mass transfer event. The mass transfer event seems to lower the molecular gradient near the core as a rejuvenation happens for both companions. Similar behaviour is present throughout the lifetime of the system, where the helium core mass fraction grows steadily until the convective core overshoot can break through the molecular gradient. This then leads to a quick rejuvenation of the core. The staggered behaviour is an effect of the current implementation. We make use of a step overshoot with a hard limit on the molecular gradient (see Sect. 5 for a discussion of the numerical behaviour).

In general, both helium mass fractions are larger than in the FO models, indicating that our overshooting scheme reduces the efficiency of rejuvenation, which has a direct effect on the evolution. As discussed in the previous section, the slow case A mass transfer is extended and multiple inversions of the mass ratio take place.

4.3. Full grid results

We evolved 4896 models (for grid description, see Sect. 3.3) and analysed the contact phases. Of the 4896 models, about 27.7% ran into numerical trouble. Most often, the time step became too small because MESA could not find a convergent model. This means that 72.3% of the models evolved as expected. Of all the models, 14.4% models contained at least one star that left the main sequence before merging. These are mostly models with large initial periods. The models with L2 overflow before the ZAMS account for 21.7%. These all have initial periods below one day. Only 5.4% of the models attained the high mass transfer termination as described above. In total, 1507 models reached L2 overflow on the main sequence (30.8%).

The weighted distribution 𝒫(qobs, pobs) is visualized in Fig. 3. This plot shows the probability distribution of the time spent in each observed mass ratio and period bin. The known contact binaries are denoted by gold stars. The uncertainties on the periods are negligible. The probability distribution for the mass ratio is shown in the upper plot. This is obtained by integrating over all period bins.

thumbnail Fig. 3.

Two-dimensional distribution of t(qobs, pobs). The gold stars denote the observed massive contact binaries. The upper plot shows the PDF of 𝒫(qobs).

The 2D histogram in Fig. 3 has the same shape as that seen in Menon et al. (2021) and Fabry (2024), as apparent from the strong correlation between q and p. However, the mass ratio distribution in this work reaches more extreme values. The probability distribution of qobs is less skewed to unity compared to previous studies. We have a clear bimodal distribution with a secondary peak around qobs = 0.5. The bimodality is also clearly visible in Fig. 4, which shows the cumulative distribution of the observational mass ratios. The CDF of the observations is shown in black. Figure 5 shows the cumulative distribution function of the observed period. The general trend of the periods is explained by our model, although we overestimate the short-period systems. We suspect this is a mass effect, as the observations contain more O+O systems, while the population is dominated by B+B systems on the grounds of the initial mass function.

thumbnail Fig. 4.

Cumulative distribution function of the observational mass ratios. The black line denotes the CDF of the observed contact binaries. The grey lines denote CDFs of the observations with the mass ratios shifted by their errors.

thumbnail Fig. 5.

Cumulative distribution function of the observational periods. The black line denotes the CDF of the observed contact binaries.

Figure 3 shows that quite a few observations fall in regions of low probability. This indicates that the model still does not fully explain the observations. A Kolmogorov-Smirnov test to compare the observed sample to the predicted distribution can be performed on both the marginal period and the mass ratio distribution. For the period distribution, we obtain a p-value of 0.024. For the mass ratio distribution, we obtain a p-value of 0.058. Considering a confidence level of 95%, we can reject the period distribution, but not the mass ratio distribution. However, this is a large increase from a similar test performed by Fabry (2024) and Fabry et al. (2025) who report a p-value of 6.2 × 10−7.

To quantify the different rejuvenation efficiency between the FO models and this work, we consider the ratio of the helium mass fraction of the core at the moment the initial secondary becomes the mass donor of the system2, and before the first interaction:

R = Y core, initial secondary becomes donor Y core, before interaction . $$ \begin{aligned} R = \frac{Y_{\text{ core,} \text{ initial} \text{ secondary} \text{ becomes} \text{ donor}}}{ Y_{\text{ core,} \text{ before} \text{ interaction}}}. \end{aligned} $$(10)

We calculate the ratio for a slice of the grid above with a single initial mass value M1, init = 24 M⊙ and display the resulting values in Fig. 6. For the models with an initial mass ratio close to unity and a short initial period, we expect the ratio R to be sensitive to the core boundary processes. We see that in this regime, the new models have noticeably higher ratios, indicating a less efficient rejuvenation. As the inhibited overshoot leads to smaller convective cores, we expect there to be lower mixing of the envelope hydrogen into the core, which is especially relevant for the accretor during mass transfer of the system and leads to weaker rejuvenation. We suspect this to be the dominant effect for short initial period systems on the ratio R. For longer initial periods, a higher helium core mass fraction and smaller core radii are expected when the mass transfer occurs, as the stars are further along their evolution. This leads to a higher denominator in the ratio, which leads to lower values.

thumbnail Fig. 6.

Ratio of the helium mass fraction of the core at the moment the initial secondary becomes the donor, and before the first interaction for the grid models with M1, init = 24 M⊙. The left plot gives the values for the Full Overshoot (FO) models of Fabry et al. (2025) and the right plot gives the parameter values for the models of this work.

4.4. Surviving systems

Of the full grid, 14.4% of the models were terminated due to one companion leaving the main sequence. These surviving systems are of interest as they contain gravitational-wave progenitors, and the remaining helium core masses are of special interest. We compare the helium core masses of the surviving systems to the helium core masses of the FO models, using the models of the grid presented in Fabry et al. (2025) within the initial condition ranges of this study. The helium core masses of both models are shown in Fig. 7, with the new models in blue and the FO models in grey compared to the initial mass. We show the core masses for the companion leaving the main sequence. For the FO models, we took the helium core to be the convective core mass at termination, while for the new models, we defined the helium core masses to be the mass coordinate where the hydrogen mass fraction drops below 10%.

thumbnail Fig. 7.

Helium core masses of the surviving systems compared to the initial mass of the stars. The new models are shown in blue and the FO in grey. A square (top panel) is used when the primary leaves the main sequence, a star (bottom panel) when the secondary leaves the main sequence.

For both implementations, the primary stars are most often the first to leave the main sequence, with the overshoot limited implementation generally leading to similar helium core masses, especially at the higher mass. At the lower initial mass end, we see a difference in the helium core masses, which is most likely due to the different definition used for the helium core mass and, as such, we see no significant difference in the helium core masses of the initial primaries that survive the main sequence. The bottom panel suggests a systematic decrease in the helium cores of the secondaries that leave the main sequence, which is expected due to inhibited core rejuvenation as the accretor during case A mass transfer.

5. Discussion and numerical behaviour

The above results on the predicted mass ratio distribution of massive contact binaries are different to what is currently presented in the literature. The additional dependence of the overshoot on the chemical gradient limits the rejuvenation of the accretor and affects the size of the receding convective cores. The presented method overpredicts systems between the observed mass ratios qobs = 0.4 − 0.6, with a clear peak in the probability distribution function at qobs = 0.5. The value of the overshoot limit (B) chosen is expected to have a large effect on the behaviour of the system; the exact impact of this value should be studied further.

One large caveat is the precise implementation of the dependence on the molecular gradient. As this is obtained through numerical differentiation, there could be a large dependence on the refinement of the grid. This is clear from the following parameter study. To study the effects of varying the smoothing parameters, we selected the model with initial values M1, init = 24 M⊙, Pinitial = 1.377 d, and qinitial = 0.8, the same initial values as the model shown in Sect. 4.1.

The parameter with the largest effect on the numerical behaviour of the MESA models is the smoothing used to calculate composition gradients (called num_cells_for_smooth_gradL_composition_term3), which takes a value of two in the grid presented above. This parameter defines the number of grid cells taken on each side over which a weighted average is taken to compute ∇μ.

Figure 8 shows the molecular gradient of the primary star at 1.75 Myr after ZAMS, which puts the model in the slow case A mass transfer event. The figure shows that there are large differences depending on the smoothing parameter. As expected, a low number of cells leads to peaked molecular gradients with smaller cores, while a larger number of cells leads to a smoother molecular gradient. However, it is unclear what the optimal value for this parameter should be. The currently implemented default of MESA is 3, which is likely to increase the numerical stability of the evolution models.

thumbnail Fig. 8.

Molecular gradients of the primary star for different values of the smoothing parameter. The model has initial values M1, initial = 24 M⊙, Pinitial = 1.377 d, and qinitial = 0.8. The profiles are selected close to 1.75 Myr. The dashed grey line denotes the overshoot limit B = 10−4.

To investigate the numerical behaviour of the smoothing further, single star models of masses 10 M⊙, 25 M⊙, 30 M⊙, and 35 M⊙ were computed with a fixed rotation rate of 0.75ωcrit and with the microphysics given in Sect. 3.1. In these models, we applied no overshoot limiting to isolate as best as possible the effect of the smoothing parameter. This was done to approach the situation of a binary model without considering any mass transfer. Figure 9 gives the core mass evolution of the different single star models. The solid lines give num_cells_for_smooth_gradL_composition_term = 2 and the dashed lines give num_cells_for_smooth _gradL_composition_term = 5.

thumbnail Fig. 9.

Core mass evolution of a single star model with fixed rotational velocity for different initial masses. The solid lines denote num_cells_for_smooth_gradL_composition_term = 2, and the dashed lines have num_cells_for_smooth_gradL _composition_term = 5.

The core mass shrinks for the models with an initial mass Minitial ≤ 30 M⊙ and a value of 2 for the smoothing parameter. The model with an initial mass Minitial = 35 M⊙ has an increasing core mass for the same model parameters. The initial masses for which such an increase happens are dependent on the value of the smoothing parameter. This can be seen by the model with an initial mass of 30 M⊙ that decreases for num_cells_for_smooth_gradL_composition_term = 2, but increases when the value is set to 5.

The large increase in the core mass during the main sequence is a numerical effect. This can be seen in Fig. 10, which shows different mesh refinements for the single star model with an initial mass Minitial = 24 M⊙ and num_cells_for_smooth_gradL_composition_term = 10. For a larger refinement, core growth does not happen, and the evolution reverts to the core receding, as expected for massive stars.

thumbnail Fig. 10.

Core mass evolution of a single star model with different mesh refinements with an initial mass of M = 24 M⊙. All models apply num_cells_for_smooth_gradL_composition_term = 10.

We find that the high mass models are sensitive to the smoothing on the molecular gradient. The sensitivity increases with higher mass, where the core can artificially increase in mass. This can be helped with a larger mesh refinement. These findings are similar to those of Moore & Garaud (2016) for the lower mass stars with a convective core (1.2 M⊙ – 1.7 M⊙.) The authors find that the core artificially shrinks for a model with smoothing enabled, which is unexpected for stars in that mass range. In the mass regime studied here, the cores artificially grow. In both cases, the cores do not behave correctly when implementing a non-zero value for the smoothing parameter. We note that, similar to Moore & Garaud (2016), we considered the Ledoux criterion to determine the core size.

The grid implemented num_cells_for_smooth _gradL_composition_term = 2. From this discussion, we see that models with an initial mass over 30 M⊙ can behave unexpectedly. However, excluding these models has little impact on the probability distribution of the observed mass ratio and periods, since the IMF heavily favours lower mass stars. We therefore chose to include all models in the grid in the analysis.

6. Conclusions

In this work, we investigated the effects of limiting the convective core overshoot on a massive contact binary population. We limited the overshoot based on the molecular weight gradient. This was accomplished by using a step overshoot with an additional requirement that the molecular weight gradient not exceed the threshold value B = 10−4. The goal of this implementation was to limit the rejuvenation efficiency during mass transfer.

We calculated a grid of 4896 models with different initial masses, mass ratios, and periods. This allowed us to investigate the resulting population and compare with observations. Our main conclusions are as follows.

  1. By limiting the overshoot, the mixing of envelope material into the core is reduced. This inhibits the rejuvenation present in contact binaries. We notice a lower rejuvenation that is staggered due to the hard limit presented by the current implementation. Throughout the grid, the rejuvenation is less efficient.

  2. Because of the reduced rejuvenation, we considerably change the evolution of a massive contact binary, with multiple contact phases and more extreme mass ratios. This is in contrast to the previously presented models in the literature. The population synthesis shows an additional peak in the PDF of the mass ratios, and we overpredict the number of unequal-mass contact binaries.

  3. High mass stellar models are sensitive to smoothing on the Ledoux gradient. This sensitivity increases with stellar mass and can make the core artificially grow. When using smoothing to aid numerical stability, the mesh density of the stellar models needs to be high enough to assure the correct behaviour of the convective boundary.

We showed that the implementation of stellar physics can have large impacts on binary population synthesis. Follow-up work could include extending the grid to different values of the overshoot limit, B. Additionally, a more gradual dependence on the molecular weight gradient could be explored. Other processes to inhibit the core rejuvenation of the accretor, such as magnetism and radiative transfer, could also be considered.

Nevertheless, to make progress in the understanding of rejuvenation, and thus of the evolution of massive contact binaries, we lack observationally motivated prescriptions for its efficiency. Asteroseismology has mainly focused on single (rotating) stars to obtain reliable convective core masses, but as of yet no asteroseismic properties of accretors in binary systems have been published. Such efforts will likely bear fruit, however, as Henneco et al. (2024) have identified that merged stars carry asteroseismic imprints that are distinguishable from true single stars.

Finally, we note that when computing massive star models, the mesh refinement and Ledoux smoothing should be carefully considered. This is of high importance for the asteroseismic estimation of the core overshoot and the core mass of massive stars (e.g. Aerts et al. 2003; Ausseloos et al. 2004; Aerts et al. 2006; Mazumdar et al. 2006; Burssens et al. 2023; Fritzewski et al. 2024). In particular, β Cephei pulsators are in the regime where this numerical behaviour is present. As shown, the choice of smoothing parameter value can strongly affect the modelled core size. Therefore, when comparing models to asteroseismic observables, systematic errors could be introduced by this artificial smoothing.

1

We make no distinction between the term ‘contact’ and ‘overcontact’.

2

We ignore small numerical jumps in the evolution.

Acknowledgments

M.F. thanks the Flemish research foundation (FWO, Fonds voor Wetenschappelijk Onderzoek) PhD fellowship No. 11H2421N for its support. The data was processed and visualised using Matplotlib (Hunter 2007), SciPy (Virtanen et al. 2020) and NumPy (Harris et al. 2020).

References

  1. Abdul-Masih, M., Sana, H., Hawcroft, C., et al. 2021, A&A, 651, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Abdul-Masih, M., Escorza, A., Menon, A., Mahy, L., & Marchant, P. 2022, A&A, 666, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Aerts, C., Thoul, A., Daszyńska, J., et al. 2003, Science, 300, 1926 [Google Scholar]
  4. Aerts, C., Marchenko, S. V., Matthews, J. M., et al. 2006, ApJ, 642, 470 [NASA ADS] [CrossRef] [Google Scholar]
  5. Aguilera-Dena, D. R., Langer, N., Moriya, T. J., & Schootemeijer, A. 2018, ApJ, 858, 115 [NASA ADS] [CrossRef] [Google Scholar]
  6. Aguilera-Dena, D. R., Langer, N., Antoniadis, J., & Müller, B. 2020, ApJ, 901, 114 [NASA ADS] [CrossRef] [Google Scholar]
  7. Alcock, C., Allsman, R. A., Alves, D., et al. 1997, AJ, 114, 326 [NASA ADS] [CrossRef] [Google Scholar]
  8. Almeida, L. A., Sana, H., de Mink, S. E., et al. 2015, ApJ, 812, 102 [Google Scholar]
  9. Almeida, L. A., Sana, H., Taylor, W., et al. 2017, A&A, 598, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Angelou, G. C., Bellinger, E. P., Hekker, S., et al. 2020, MNRAS, 493, 4987 [NASA ADS] [CrossRef] [Google Scholar]
  11. Angulo, C., Arnould, M., Rayet, M., et al. 1999, Nucl. Phys. A, 656, 3 [Google Scholar]
  12. Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
  13. Ausseloos, M., Scuflaire, R., Thoul, A., & Aerts, C. 2004, MNRAS, 355, 352 [NASA ADS] [CrossRef] [Google Scholar]
  14. Banyard, G., Sana, H., Mahy, L., et al. 2022, A&A, 658, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Bastian, N., Covey, K. R., & Meyer, M. R. 2010, ARA&A, 48, 339 [Google Scholar]
  16. Biermann, P., & Thomas, H. C. 1972, A&A, 16, 60 [NASA ADS] [Google Scholar]
  17. Blouin, S., Shaffer, N. R., Saumon, D., & Starrett, C. E. 2020, ApJ, 899, 46 [NASA ADS] [CrossRef] [Google Scholar]
  18. Böhm-Vitense, E. 1958, Z. Astrophys., 46, 108 [Google Scholar]
  19. Bohm-Vitense, E., & Nelson, G. D. 1976, ApJ, 210, 741 [Google Scholar]
  20. Brott, I., de Mink, S. E., Cantiello, M., et al. 2011, A&A, 530, A115 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  21. Burssens, S., Bowman, D. M., Michielsen, M., et al. 2023, Nat. Astron., 7, 913 [NASA ADS] [CrossRef] [Google Scholar]
  22. Çakırlı, Ö., Ibanoglu, C., & Sipahi, E. 2014, MNRAS, 442, 1560 [CrossRef] [Google Scholar]
  23. Cassisi, S., Potekhin, A. Y., Pietrinferni, A., Catelan, M., & Salaris, M. 2007, ApJ, 661, 1094 [NASA ADS] [CrossRef] [Google Scholar]
  24. Chugunov, A. I., Dewitt, H. E., & Yakovlev, D. G. 2007, Phys. Rev. D, 76, 025028 [NASA ADS] [CrossRef] [Google Scholar]
  25. Cox, J. P., & Giuli, R. T. 1968, Principles of stellar structure (New York: Gordon and Breach) [Google Scholar]
  26. Cyburt, R. H., Amthor, A. M., Ferguson, R., et al. 2010, ApJS, 189, 240 [NASA ADS] [CrossRef] [Google Scholar]
  27. de Mink, S. E., & Pols, O. R. 2007, in Binary Stars as Critical Tools& Tests in Contemporary Astrophysics, eds. W. I. Hartkopf, P. Harmanec, & E. F. Guinan, 240, 384 [Google Scholar]
  28. Dray, L. M., & Tout, C. A. 2007, MNRAS, 376, 61 [Google Scholar]
  29. Dunstall, P. R., Dufton, P. L., Sana, H., et al. 2015, A&A, 580, A93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Evans, C. J., Smartt, S. J., Lee, J. K., et al. 2005, A&A, 437, 467 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Fabry, M. 2024, Ph.D. Thesis, KU Leuven, Belgium. Available at, https://fys.kuleuven.be/ster/pub/phd-thesis-matthias-fabry/phd-thesis-final.pdf [Google Scholar]
  32. Fabry, M., Marchant, P., & Sana, H. 2022, A&A, 661, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Fabry, M., Marchant, P., Langer, N., & Sana, H. 2023, A&A, 672, A175 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Fabry, M., Marchant, P., Langer, N., & Sana, H. 2025, A&A, 695, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585 [Google Scholar]
  36. Flannery, B. P. 1976, ApJ, 205, 217 [NASA ADS] [CrossRef] [Google Scholar]
  37. Förster, K. 2023, Bachelor’s thesis, Universität Bonn, Germany. Available at https://astro.uni-bonn.de/~nlanger/siu_web/pre.html [Google Scholar]
  38. Fritzewski, D. J., Vanrespaille, M., Aerts, C., et al. 2024, A&A, submitted [arXiv:2408.06097] [Google Scholar]
  39. Frost, A. J., Sana, H., Mahy, L., et al. 2024, Science, 384, 214 [NASA ADS] [CrossRef] [Google Scholar]
  40. Fuller, G. M., Fowler, W. A., & Newman, M. J. 1985, ApJ, 293, 1 [NASA ADS] [CrossRef] [Google Scholar]
  41. Hamann, W. R., Koesterke, L., & Wessolowski, U. 1995, A&A, 299, 151 [Google Scholar]
  42. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
  43. Hazlehurst, J., & Refsdal, S. 1980, A&A, 84, 200 [NASA ADS] [Google Scholar]
  44. Heger, A., Langer, N., & Woosley, S. E. 2000, ApJ, 528, 368 [NASA ADS] [CrossRef] [Google Scholar]
  45. Henneco, J., Schneider, F. R. N., & Laplace, E. 2024, A&A, 682, A169 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Hilditch, R. W., & Evans, T. L. 1985, MNRAS, 213, 75 [CrossRef] [Google Scholar]
  47. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
  48. Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [Google Scholar]
  49. Iglesias, C. A., & Rogers, F. J. 1993, ApJ, 412, 752 [Google Scholar]
  50. Iglesias, C. A., & Rogers, F. J. 1996, ApJ, 464, 943 [NASA ADS] [CrossRef] [Google Scholar]
  51. Itoh, N., Hayashi, H., Nishikawa, A., & Kohyama, Y. 1996, ApJS, 102, 411 [NASA ADS] [CrossRef] [Google Scholar]
  52. Jermyn, A. S., Schwab, J., Bauer, E., Timmes, F. X., & Potekhin, A. Y. 2021, ApJ, 913, 72 [NASA ADS] [CrossRef] [Google Scholar]
  53. Jermyn, A. S., Bauer, E. B., Schwab, J., et al. 2023, ApJS, 265, 15 [NASA ADS] [CrossRef] [Google Scholar]
  54. Johnston, C., Michielsen, M., Anders, E. H., et al. 2024, ApJ, 964, 170 [NASA ADS] [Google Scholar]
  55. Kippenhahn, R., Kohl, K., & Weigert, A. 1967, Z. Astrophys., 66, 58 [NASA ADS] [Google Scholar]
  56. Kippenhahn, R., Ruschenplatt, G., & Thomas, H. C. 1980, A&A, 91, 175 [Google Scholar]
  57. Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
  58. Kuiper, G. P. 1941, ApJ, 93, 133 [NASA ADS] [CrossRef] [Google Scholar]
  59. Langanke, K., & Martínez-Pinedo, G. 2000, Nucl. Phys. A, 673, 481 [NASA ADS] [CrossRef] [Google Scholar]
  60. Langer, N., Fricke, K. J., & Sugimoto, D. 1983, A&A, 126, 207 [NASA ADS] [Google Scholar]
  61. Ledoux, P. 1947, ApJ, 105, 305 [NASA ADS] [CrossRef] [Google Scholar]
  62. Leung, K. C., Sistero, R. F., Zhai, D. S., Grieco, A., & Candellero, B. 1984, AJ, 89, 872 [NASA ADS] [CrossRef] [Google Scholar]
  63. Li, F. X., Qian, S. B., Jiao, C. L., & Ma, W. W. 2022, ApJ, 932, 14 [NASA ADS] [CrossRef] [Google Scholar]
  64. Linnell, A. P., & Scheick, X. 1991, ApJ, 379, 721 [NASA ADS] [CrossRef] [Google Scholar]
  65. Lorenz, R., Mayer, P., & Drechsel, H. 1999, A&A, 345, 531 [NASA ADS] [Google Scholar]
  66. Lorenzo, J., Simón-Díaz, S., Negueruela, I., et al. 2017, A&A, 606, A54 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Lucy, L. B. 1968, ApJ, 151, 1123 [Google Scholar]
  68. Mahy, L., Almeida, L. A., Sana, H., et al. 2020a, A&A, 634, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Mahy, L., Sana, H., Abdul-Masih, M., et al. 2020b, A&A, 634, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  70. Mandel, I., & de Mink, S. E. 2016, MNRAS, 458, 2634 [NASA ADS] [CrossRef] [Google Scholar]
  71. Mandel, I., & Farmer, A. 2022, Phys. Rep., 955, 1 [NASA ADS] [CrossRef] [Google Scholar]
  72. Marchant, P., Langer, N., Podsiadlowski, P., Tauris, T. M., & Moriya, T. J. 2016, A&A, 588, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  73. Martins, F., Mahy, L., & Hervé, A. 2017, A&A, 607, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Mayer, P., Drechsel, H., Harmanec, P., Yang, S., & Šlechta, M. 2013, A&A, 559, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  75. Mazumdar, A., Briquet, M., Desmet, M., & Aerts, C. 2006, A&A, 459, 589 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Menon, A., Langer, N., de Mink, S. E., et al. 2021, MNRAS, 507, 5013 [NASA ADS] [CrossRef] [Google Scholar]
  77. Michielsen, M., Aerts, C., & Bowman, D. M. 2021, A&A, 650, A175 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Moore, K., & Garaud, P. 2016, ApJ, 817, 54 [Google Scholar]
  79. Moravveji, E., Aerts, C., Pápics, P. I., Triana, S. A., & Vandoren, B. 2015, A&A, 580, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  80. Nelson, C. A., & Eggleton, P. P. 2001, ApJ, 552, 664 [NASA ADS] [CrossRef] [Google Scholar]
  81. Nieuwenhuijzen, H., & de Jager, C. 1995, A&A, 302, 811 [NASA ADS] [Google Scholar]
  82. Oda, T., Hino, M., Muto, K., Takahara, M., & Sato, K. 1994, Atomic Data Nucl. Data Tables, 56, 231 [NASA ADS] [CrossRef] [Google Scholar]
  83. Öpik, E. 1924, Publications of the Tartu Astrofizica Observatory, 25, 1 [Google Scholar]
  84. Ostrov, P. G. 2001, MNRAS, 321, L25 [NASA ADS] [CrossRef] [Google Scholar]
  85. Papaloizou, J., & Pringle, J. E. 1979, MNRAS, 189, 5P [NASA ADS] [CrossRef] [Google Scholar]
  86. Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
  87. Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [Google Scholar]
  88. Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
  89. Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34 [NASA ADS] [CrossRef] [Google Scholar]
  90. Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10 [Google Scholar]
  91. Pedersen, M. G., Aerts, C., Pápics, P. I., et al. 2021, Nat. Astron., 5, 715 [NASA ADS] [CrossRef] [Google Scholar]
  92. Penny, L. R., Ouzts, C., & Gies, D. R. 2008, ApJ, 681, 554 [Google Scholar]
  93. Pols, O. R. 1994, A&A, 290, 119 [Google Scholar]
  94. Potekhin, A. Y., & Chabrier, G. 2010, Contrib. Plasma Phys., 50, 82 [NASA ADS] [CrossRef] [Google Scholar]
  95. Poutanen, J. 2017, ApJ, 835, 119 [NASA ADS] [CrossRef] [Google Scholar]
  96. Raucq, F., Gosset, E., Rauw, G., et al. 2017, A&A, 601, A133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  97. Rogers, F. J., & Nayfonov, A. 2002, ApJ, 576, 1064 [Google Scholar]
  98. Salpeter, E. E. 1955, ApJ, 121, 161 [Google Scholar]
  99. Sana, H., de Mink, S. E., de Koter, A., et al. 2012, Science, 337, 444 [Google Scholar]
  100. Sana, H., Le Bouquin, J. B., Lacour, S., et al. 2014, ApJS, 215, 15 [Google Scholar]
  101. Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713 [NASA ADS] [CrossRef] [Google Scholar]
  102. Schneider, F. R. N., Podsiadlowski, P., Langer, N., Castro, N., & Fossati, L. 2016, MNRAS, 457, 2355 [Google Scholar]
  103. Schneider, F. R. N., Ohlmann, S. T., Podsiadlowski, P., et al. 2019, Nature, 574, 211 [Google Scholar]
  104. Schootemeijer, A., Langer, N., Grin, N. J., & Wang, C. 2019, A&A, 625, A132 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  105. Schwarzschild, M., & Härm, R. 1958, ApJ, 128, 348 [NASA ADS] [CrossRef] [Google Scholar]
  106. Shu, F. H., Lubow, S. H., & Anderson, L. 1976, ApJ, 209, 536 [NASA ADS] [CrossRef] [Google Scholar]
  107. Timmes, F. X., & Swesty, F. D. 2000, ApJS, 126, 501 [NASA ADS] [CrossRef] [Google Scholar]
  108. Townsend, R. 2023, https://doi.org/10.5281/zenodo.10624843 [Google Scholar]
  109. Vanbeveren, D., De Donder, E., Van Bever, J., Van Rensbergen, W., & De Loore, C. 1998, New Astron., 3, 443 [CrossRef] [Google Scholar]
  110. Viani, L. S., & Basu, S. 2020, ApJ, 904, 22 [NASA ADS] [CrossRef] [Google Scholar]
  111. Villaseñor, J. I., Taylor, W. D., Evans, C. J., et al. 2021, MNRAS, 507, 5348 [CrossRef] [Google Scholar]
  112. Vink, J. S., de Koter, A., & Lamers, H. J. G. L. M. 2001, A&A, 369, 574 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  113. Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
  114. Vitense, E. 1953, Z. Astrophys., 32, 135 [NASA ADS] [Google Scholar]
  115. Vrancken, J., Abdul-Masih, M., Escorza, A., et al. 2024, A&A, 691, A150 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  116. Webbink, R. F. 1977, ApJ, 215, 851 [Google Scholar]
  117. Wellstein, S., Langer, N., & Braun, H. 2001, A&A, 369, 939 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  118. Yang, Y., Yuan, H., & Dai, H. 2019, AJ, 157, 111 [NASA ADS] [CrossRef] [Google Scholar]
  119. Yaşarsoy, B., & Yakut, K. 2014, New Astron., 31, 32 [CrossRef] [Google Scholar]
  120. Yoon, S. C., & Langer, N. 2005, A&A, 443, 643 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  121. Yoon, S. C., Langer, N., & Norman, C. 2006, A&A, 460, 199 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  122. Zahn, J. P. 1975, A&A, 41, 329 [NASA ADS] [Google Scholar]

Appendix A: HR diagram of the example evolution

The evolution of a contact binary system is heavily influenced by the adjusted rejuvenation efficiency presented above. Figure A.1 shows the HR diagram of the model presented in Sect. 4.1, where we see the evolution of the primary (orange) and the secondary (red). The tracks in the HR diagram are strongly deviated from models previously presented in the literature due to the multiple strong mass transfer events. We recognize the initial discontiuity of the case-A mass transfer, where the first contact phase takes place. The inversions of the mass ratios are clearly visible in the global evolution of the system, with the primary and secondary track switching between the upper left and lower right of the HR diagram. Around 3 Myr, the system undergoes a ‘bump’ in the mass ratio, where there is a quick evolution to equal mass ratio, but it does not invert and decreases back to a lower mass ratio. This characteristic is also visible in the HR diagram, where there are two distict loops halfway through the evolution, which cross around log Teff ∼ 4.515 and log L/L⊙ ∼ 4.8.

thumbnail Fig. A.1.

HR diagram of the model presented in Sect. 4.1 for both the primary (orange) and the secondary (red). The contact phases are shown with a lower opacity. The ZAMS is highlighted by a star. The black lines denote the three mass transfer events approximating thermal mass transfer rate.

Appendix B: Observations

Table B.1.

Parameters of observed, massive (near-)contact systems in the Milky Way, M31, and the Magellanic Clouds.

All Tables

Table B.1.

Parameters of observed, massive (near-)contact systems in the Milky Way, M31, and the Magellanic Clouds.

In the text

All Figures

thumbnail Fig. 1.

Example of the mass ratio evolution with initial conditions of M1, initial = 24 M⊙, Pinitial = 1.377 d, and qinitial = 0.8. Top panel: Mass ratio evolution of the system. The coloured lines show the new model, the grey line shows the FO mass ratio evolution. The dashed lines show the stellar radii normalized by the respective Roche lobe radii. Bottom panel: Mass transfer rate of the system of the new model (blue) and the FO model (grey). The dashed black line denotes the initial thermal mass transfer rate of the secondary.
In the text
thumbnail Fig. 2.

Evolution of the helium mass fraction in the core of the star. The upper plot shows the primary star, the lower plot shows the secondary star. The grey line denotes the FO model.
In the text
thumbnail Fig. 3.

Two-dimensional distribution of t(qobs, pobs). The gold stars denote the observed massive contact binaries. The upper plot shows the PDF of 𝒫(qobs).
In the text
thumbnail Fig. 4.

Cumulative distribution function of the observational mass ratios. The black line denotes the CDF of the observed contact binaries. The grey lines denote CDFs of the observations with the mass ratios shifted by their errors.
In the text
thumbnail Fig. 5.

Cumulative distribution function of the observational periods. The black line denotes the CDF of the observed contact binaries.
In the text
thumbnail Fig. 6.

Ratio of the helium mass fraction of the core at the moment the initial secondary becomes the donor, and before the first interaction for the grid models with M1, init = 24 M⊙. The left plot gives the values for the Full Overshoot (FO) models of Fabry et al. (2025) and the right plot gives the parameter values for the models of this work.
In the text
thumbnail Fig. 7.

Helium core masses of the surviving systems compared to the initial mass of the stars. The new models are shown in blue and the FO in grey. A square (top panel) is used when the primary leaves the main sequence, a star (bottom panel) when the secondary leaves the main sequence.
In the text
thumbnail Fig. 8.

Molecular gradients of the primary star for different values of the smoothing parameter. The model has initial values M1, initial = 24 M⊙, Pinitial = 1.377 d, and qinitial = 0.8. The profiles are selected close to 1.75 Myr. The dashed grey line denotes the overshoot limit B = 10−4.
In the text
thumbnail Fig. 9.

Core mass evolution of a single star model with fixed rotational velocity for different initial masses. The solid lines denote num_cells_for_smooth_gradL_composition_term = 2, and the dashed lines have num_cells_for_smooth_gradL _composition_term = 5.
In the text
thumbnail Fig. 10.

Core mass evolution of a single star model with different mesh refinements with an initial mass of M = 24 M⊙. All models apply num_cells_for_smooth_gradL_composition_term = 10.
In the text
thumbnail Fig. A.1.

HR diagram of the model presented in Sect. 4.1 for both the primary (orange) and the secondary (red). The contact phases are shown with a lower opacity. The ZAMS is highlighted by a star. The black lines denote the three mass transfer events approximating thermal mass transfer rate.
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.