The majority of close massive binary stars with initial periods of a few days experience a contact phase in which both stars overflow their Roche lobes simultaneously. We perform the first dedicated study of the evolution of massive contact binaries and provide a comprehensive prediction of their observed properties. We compute 2790 detailed binary models for the Large and Small Magellanic Clouds each assuming mass transfer to be conservative. The initial parameter space for both grids span total masses from 20 to 80 mathrmModot  orbital periods of 0.62 d and mass ratios of 0.61.0. We find that models that remain in contact over nuclear timescales evolve towards equal masses echoing the mass ratios of their observed counterparts. Ultimately the fate of our nucleartimescale models is to merge on the main sequence. Our predicted periodmass ratio distributions of Otype contact binaries are similar for both galaxies and we expect 10 such systems together in both Magellanic Clouds. While we can largely reproduce the observed distribution we overestimate the population of equalmass contact binaries. This situation is somewhat remedied if we also account for binaries that are nearly in contact. Our theoretical distributions work particularly well for contact binaries with periods lt2 d and total masses lessapprox 45 mathrmModot . We expect stellar winds nonconservative mass transfer and envelope inflation to have played a role in the formation of the more massive and longerperiod contact binaries.
['binaries: close', 'binaries: spectroscopic', 'stars: early-type', 'stars: evolution', 'stars: fundamental parameters', 'stars: fundamental', 'Astrophysics - Solar and Stellar Astrophysics', 'Astrophysics - Astrophysics of Galaxies']
Detailed evolutionary models of massive contact binaries  I. Model grids and synthetic populations for the Magellanic Clouds
{'Detailed evolutionary models of massive contact binaries I. Model grids and synthetic populations for the Magellanic Clouds': 'Athira Menon 1 GLYPH<150> 2 ★ , Norbert Langer 1 GLYPH<150> 3 , Selma E. de Mink 4 GLYPH<150> 2 GLYPH<150> 5 , Stephen Justham 2 GLYPH<150> 6 GLYPH<150> 7 GLYPH<150> 4 , Koushik Sen 1 GLYPH<150> 3 , Dorottya Szécsi 8 GLYPH<150> 9 , Alex de Koter 2 GLYPH<150> 12 , Michael Abdul-Masih 10 , Hugues Sana 12 , Laurent Mahy 11 GLYPH<150> 12 , Pablo Marchant 12  \n- 1 Argelander-Institut für Astronomie, Universität Bonn, Auf dem Hügel 71, 53121 Bonn, Germany\n- 2 Anton Pannekoek Institute of Astronomy and GRAPPA, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, Netherlands\n- 3 Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany\n- 4 Max Planck Institute for Astrophysics, Karl-Schwarzschild-Strasse 1, 85748 Garching, Germany\n- 5 Center for Astrophysics, Harvard-Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA\n- 6 School of Astronomy and Space Science, University of the Chinese Academy of Sciences, Beĳing 100012, China\n- 7 National Astronomical Observatories, Chinese Academy of Sciences, Beĳing 100012, China\n- 8 Institute of Astronomy, Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University, Grudziądzka 5, 87-100 Torun, Poland\n- 9 Physikalisches Institut, Universität zu Köln, Züpicher-Strasse 77, D-50937 Cologne, Germany  \n10  \nEuropean Southern Observatory, Alonso de Cordova 3107, Vitacura, Casilla 19001, Santiago de Chile, Chile  \n- 11 Royal Observatory of Belgium, Avenue circulaire 3, B-1180 Brussel, Belgium\n- 12 Institute of Astrophysics, KU Leuven, Celestĳnenlaan 200D, 3001 Leuven, Belgium.  \nAccepted XXX. Received YYY; in original form ZZZ', 'ABSTRACT': 'The majority of close massive binary stars with initial periods of a few days experience a contact phase, in which both stars overflow their Roche lobes simultaneously. We perform the first dedicated study of the evolution of massive contact binaries and provide a comprehensive prediction of their observed properties. We compute 2790 detailed binary models for the Large and Small Magellanic Clouds each, assuming mass transfer to be conservative. The initial parameter space for both grids span total masses from 20 to 80 M GLYPH<12> , orbital periods of 0.6 to 2 days and mass ratios of 0.6 to 1.0. We find that models that remain in contact over nuclear timescales evolve towards equal masses, echoing the mass ratios of their observed counterparts. Ultimately, the fate of our nuclear-timescale models is to merge on the main sequence. Our predicted period-mass ratio distributions of O-type contact binaries are similar for both galaxies, and we expect 10 such systems together in both Magellanic Clouds. While we can largely reproduce the observed distribution, we over-estimate the population of equal-mass contact binaries. This situation is somewhat remedied if we also account for binaries that are nearly in contact. Our theoretical distributions work particularly well for contact binaries with periods GLYPH<159> 2 days and total masses / 45 M GLYPH<12> . We expect stellar winds, non-conservative mass transfer and envelope inflation to have played a role in the formation of the more massive and longer-period contact binaries.  \nKey words: stars: massive - stars: early-type - binaries:close - binaries: general - stars: evolution - stars: fundamental parameters', '1 INTRODUCTION': 'A significant fraction of all stars in the Universe live in binary systems with companions that they are gravitationally bound to. Surveys of the Milky Way find that the fraction of low-mass stars found in such binary systems is 30-40% (e.g., Duquennoy & Mayor 1991; Eggleton & Tokovinin 2008; Raghavan et al. 2010; Duchêne  \n&Kraus 2013), while for massive stars this fraction goes up to 50100% (e.g., Vanbeveren et al. 1998; Chini et al. 2012; Sana et al. 2012, 2014; Dunstall et al. 2015; Kobulnicky et al. 2014; Aldoretta et al. 2015; Moe & Di Stefano 2017).  \nThelargest homogeneous sample of O-type stars (M GLYPH<21> 15 M GLYPH<12> ) wehave so far is from the young star-forming region, 30 Doradus (or Tarantula nebula), nestled in the Large Magellanic Cloud. The VLTFLAMESTarantula survey (VFTS, Evans et al. 2011) obtained data for nearly 800 massive stars in 30 Dor, 360 of which are O-type stars  \n(Sana et al. 2013) and 116 of which are binaries. A follow up study called the Tarantula Massive Binary Monitoring (TMBM) project calculated orbital solutions for nearly 82 of these systems (Almeida et al. 2017; Mahy et al. 2020a,b), thus providing us with a rich, unbiased database to study the lives of O stars in binary systems.  \nA large fraction of early-type massive binaries in both, the Galactic and the VFTS samples, have orbital periods of a few days. These systems are expected to evolve into a contact configuration (Wilson 2001), wherein both stars overflow their Roche lobes simultaneously and become bound by a common equipotential surface that extends between the inner and outer Lagrangian points of either star, lending them their unique peanut-shaped structure.  \nInvestigating the evolution of massive contact binaries is essential to estimate the fraction of the binary population that will merge on the main sequence or evolve further until they explode or collapse to black holes. Those that do merge before either star forms a compact object may give rise to fast-rotating single stars which are considered to be progenitors of long duration 𝛾 GLYPH<0> ray bursts (Yoon &Langer 2005; Woosley & Bloom 2006; Yoon et al. 2006; Meynet & Maeder 2007; Dessart et al. 2008; van Marle, A. J. et al. 2008; Szécsi 2017; Aguilera-Dena et al. 2018), B[e] and sgB[e] stars (e.g., Podsiadlowski et al. 2006; Vanbeveren et al. 2013; Justham et al. 2014; Wu et al. 2020), blue supergiant progenitors of Type II supernovae such as SN 1987A (Podsiadlowski et al. 1992; Menon & Heger 2017; Urushibata et al. 2018; Menon et al. 2019), progenitors of superluminous supernovae (Justham et al. 2014; Sukhbold et al. 2016; Aguilera-Dena et al. 2020), pulsational pair instability supernovae (Langer 1991; Heger et al. 2003; Langer et al. 2007; Woosley et al. 2007; Chatzopoulos & Wheeler 2012; Vigna-Gómez et al. 2019) and magnetic stars (Schneider et al. 2016, 2019).  \nLow-mass contact binaries, called W UMa binaries, have been observed quite frequently- data is available for hundreds of these systems, and they have orbital periods between 0.3 and 1 day, mass ratios between 0.1 and 0.8, and total masses of approximately 1 M GLYPH<12> (Szymanski et al. 2001; Selam 2004; Rucinski et al. 2007). In comparison, fewer massive contact binaries are known from observations. Surveys of O& B-type stars have so far reported nearly 40 such systems in the Large Magellanic Cloud (LMC), 2 such systems in the Small Magellanic Cloud (SMC) and 17 in the Milky Way. They have total masses from 17 to 85 M GLYPH<12> , orbital periods from 0.45 to 6.6 days and mass ratios between 0.3 and 1. The Magellanic Cloud contact binaries are much less wider than those in the Milky Way, and are also more confined in the period-mass ratio space than the more scattered Milky Way systems.  \nObservationally, contact binaries are identified by their distinct light curve shapes. Unfortunately in most systems, the signature of contact is not always clear-cut since their degree of contact is uncertain, and could imply that they may only be nearing contact and not be in contact at present. In some cases the same system has a different status assigned by different works that studied them. We will look at these uncertainties in more detail later in the paper. Prominent among the handful of systems which are definitively in deep (over) contact, is VFTS 352 in 30 Dor, one of the most massive contact binaries known. The timescale over which the orbital period of this over-contact system varies is / 𝑃 GLYPH<157> 𝑃 GLYPH<25> 0 GLYPH<149> 6 Myr, indicating an ongoing nuclear-timescale mass transfer (Almeida et al. 2015). Similar time derivatives have also been inferred for several other massive contact binaries as well (e.g., Qian et al. 2006, 2007; Martins et al. 2017). The majority of these deep contact binaries have equal-mass components, including VFTS 352.  \nDespite the data available for massive contact binaries, there are at present no known detailed models that thoroughly investigate  \ntheir evolution or explain their observed stellar and orbital parameter distribution. In contrast, there have been several evolutionary models for low-mass contact binaries which explain their observed orbital period and mass ratio distribution and, their spread along the main sequence (e.g., Webbink 1976; Yakut & Eggleton 2005; Stepien 2006; Gazeas & Ste¸pień 2008; Stępień & Gazeas 2012; Zhang et al. 2020). An important aspect of these models is the heat transfer in their common convective envelope, which leads to large-scale circulations of mass and thermal energy in the envelope and helps explain the temperature and luminosity differences between the component stars of the binary (Lucy 1968; Flannery 1976; Webbink 1977; Ste¸pień 2009).  \nContrary to the low-mass contact binaries, massive contact binaries share a common envelope that is largely radiative and the majority of them are located close to ZAMS. The physics of their common radiative envelope is not well understood, especially the heat transfer that may occur between their components, and is hence not included while calculating the evolution of massive contact binaries.  \nThere have been some exploratory calculations for massive contact binaries, the earliest among them published in the series of papers by Sybesma (Sybesma 1985, 1986a,b). The next set of works such as those of Nelson & Eggleton (2001), Wellstein et al. (2001) and de Mink et al. (2007) explored the parameter space in which binaries can enter contact during Case A mass transfer and the impact of semi-convection and mass-transfer efficiency on the evolution of their models. de Mink et al. (2007) computed a large grid of binary models at the SMC metallicity and identified which main-sequence binaries will undergo contact and which will avoid contact. They also evolved some models that were in shallow contact. Vanbeveren et al. (1998); Mennekens & Vanbeveren (2017) also investigated the evolution of contact binaries, but only of thermal-timescale contact binaries. They stopped the evolution of those models that were expected to enter contact over a nuclear-timescale, as they expected these systems to merge eventually. Of particular relevance to our work are the works of Marchant et al. (2016, 2017), who computed grids of binary models for a range of metallicities (Z GLYPH<12> /4...Z GLYPH<12> /50) and initial binary parameters. We will discuss their work shortly.  \nAn important open question concerning the evolution of close binaries, is the efficiency of internal mixing due to rotational instabilities. Rotationally-induced mixing is predominantly due to the circulation of meridional currents in the outer radiative layers of a star (Eddington 1926, 1929) and is expected to become stronger with increasing rotational velocity, increasing initial mass, and decreasing metallicity (e.g.,Heger et al. 2000; Heger & Langer 2000; Maeder 2009; Brott et al. 2011a,b.) Mixing in the radiative outer layers causes the transport of nuclear-burning products such as helium and nitrogen from the convective core to the surface of the star. In its most extreme form, rotational mixing can cause the star to undergo chemically homogeneous evolution (CHE), wherein mixing is so efficient that the stars lose their core-envelope stratification soon after ZAMS, and remain hot and compact throughout their evolution. This channel of evolution is considered to be particularly important in the formation of rapidly-rotating metal-poor single stars (Maeder 1987; Yoon & Langer 2005; Yoon et al. 2006; Brott et al. 2011a; Szécsi et al. 2015).  \nIn the case of binaries, CHE has been suggested to occur in very short-period systems on the main sequence, in which the stars become deformed by tides (de Mink et al. 2009; Song et al. 2013, 2016; Hastings et al. 2020). In such systems, this type of evolution can prevent the stars from merging and they evolve as massive (and compact) helium stars which end their lives as heavy black holes  \n(de Mink et al. 2009; Song et al. 2016; de Mink & Mandel 2016; Mandel & de Mink 2016; Marchant et al. 2016; du Buisson et al. 2020). CHE has also been considered for near-contact binaries (de Mink et al. 2009; Song et al. 2016; Riley et al. 2020) and overcontact binaries (Marchant et al. 2016). In particular, the models computed by Marchant et al. (2016, 2017) were found to undergo CHE during long contact phases, as their metallicity decreased (Z GLYPH<159> Z GLYPH<12> /2), primary masses increased (M 1 GLYPH<21> 40 M GLYPH<12> ) and for initial mass ratios closer to unity.  \nThis paper is intended as the first in a series that will have the ultimate goal of (1) improving our understanding of the physics of massive contact binaries, (2) gaining insight in their evolutionary pathways and observable properties and (3) providing reliable predictions for their final fate as stellar mergers or possibly massive double helium stars and gravitational-wave progenitors.  \nIn this first paper we focus on very close binaries and follow their evolution through contact on the main sequence. We perform stellar evolution calculations for a large grid of short-period binaries varying their initial masses, orbital periods and mass ratios, with an initial composition that is typical for young stars in the LMC and SMC. We study the evolution of massive binaries from ZAMS and through the contact phase and, in particular, the changes in their orbital properties prior and during the contact phase. We make predictions for the distributions of their orbital parameters and provide a first comparison with observed massive contact systems. We show where our predictions agree well with observations, but also discuss tensions that exist that can provide us with further insight.  \nOur paper is organized as follows: Section 2 describes our model assumptions, the initial parameter space of the grids and our method for computing orbital parameter distributions. There are two results sections: Section 3 discusses how massive contact systems typically form and evolve in our grid and how the choice of initial parameters affect their evolution. Section 4 compares the distributions calculated from our models with the observed distributions, including the number of massive contact binaries we expect in the Magellanic Clouds. In Section 5 we present a discussion of our results and point out avenues where we can improve them, and in Section 6, we summarize our findings and indicate directions for future work.', '2 METHODS': 'To compute a large grid of binary-star models suitable for comparison with observed systems in the LMC and SMC, we use version 10398 of the 1D stellar evolution code MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019) 1 . In the remainder of this section we discuss our physical assumptions (Section 2.1), the initialization and termination of our models (Section 2.2.2) and our assumptions for modelling a population of contact binaries (Section 2.3).', '2.1.1 Mixing, mass loss, microphysics': 'Various mixing processes are implemented as time-dependent diffusive processes in MESA , where the user has the freedom to vary  \ntheir efficiencies. Convection is modelled using the standard mixinglength theory (Böhm-Vitense 1958; Cox & Giuli 1968) prescription with a mixing length parameter of 𝛼 MLT = 1 GLYPH<149> 5, following Pols et al. (1998). We adopted the Ledoux criterion for convective stability (Ledoux 1958). We allow for mixing beyond the edge of convectively unstable regions by assuming a step overshooting parameter 𝛼 ov = 0 GLYPH<149> 335 as calibrated by Brott et al. (2011a). Semiconvection is included with an efficiency parameter of 𝛼 sc = 1 GLYPH<149> 0 as in Langer et al. (1983). This value of this parameter can significantly affect the evolution of the accreting binary components (Braun & Langer 1995). Schootemeĳer et al. (2019) find a good agreement between stellar models and the distribution of massive stars in the HertzsprungRussel diagram of the SMC for 𝛼 sc = 1 GLYPH<149> 0. Thermohaline mixing is included by adopting a value of 1.0 for the dimensionless parameter 𝛼 th . The sources of rotationally-induced mixing and angular momentum transport considered in MESA include the EddingtonSweet circulation, secular and dynamic shear instabilities, and the Goldreich-Schubert-Fricke instability (Heger et al. 2000, 2005). The efficiency factor for the total rotationally-induced mixing coefficient is set to fc = 1 GLYPH<157> 30 as determined theoretically by Chaboyer & Zahn (1992). The factor controlling the inhibition of rotational mixing against gradients in mean molecular weight is set to f 𝜇 = 0 GLYPH<149> 1 as in Yoon et al. (2006). Angular momentum transport due to magnetic fields is implemented as in Heger et al. (2005) and Petrovic et al. (2005), which is motivated by observations (e.g.,Suĳs et al. 2008).  \nStellar winds are implemented depending on the surface helium mass fraction ( 𝑌 s) as in Yoon et al. (2006); if 𝑌 s GLYPH<159> 0 GLYPH<149> 4, i.e., for hydrogen-rich stars, mass-loss rates are computed as in Vink et al. (2001), while for hydrogen-poor Wolf-Rayet stars (with 𝑌 s ¡ 0 GLYPH<149> 7) the Hamann et al. (1995) prescription is used but reduced by a factor of 10. For both rates, we use a metallicity-dependent wind whose strength is proportional to 𝑍 0 GLYPH<149> 86 . For stars with 0 GLYPH<149> 4 GLYPH<159> 𝑌 s GLYPH<159> 0 GLYPH<149> 7, the mass-loss rates are interpolated between the above two prescriptions.  \nNucleosynthesis reactions are calculated using the MESA approx21.net network containing 21 isotopes from 1 H until 56 Ni and encompasses the necessary reactions for the hydrogen CNO burning cycle, thereby allowing us to follow the evolutionary phase of interest. Radiative opacities are calculated using CO-enhanced opacity tables from the OPAL project (Iglesias & Rogers 1996).', '2.1.2 Binary star physics assumptions': 'We use the binary module of MESA as described in Paxton et al. (2015) and Marchant et al. (2016) to compute the evolution of interacting binary models in our grid. Tidal synchronization is implemented according to the Hut (1981) and Hurley et al. (2002) prescriptions. Mass transfer due to Roche lobe overflow is implicitly calculated using the Ritter scheme (Ritter 1988) and is switched to the contact scheme when both stars overfill their Roche lobes, as will be discussed in Section 2.1.3. As we are studying close binaries in this work, we assume that both stars in the system are synchronized to their initial orbital period at the beginning of the evolution but we allow for differential rotation.  \nMass transfer is treated conservatively except due to mass loss through winds. We apply the standard mass transfer physics of MESA (Paxton et al. 2015) which computes the mass transfer rate through an implicit scheme, including angular momentum accretion and tidal spin-orbit coupling. The kinetic energy of the accretion stream is neglected (Ulrich & Burger 1976) and the entropy of the accreted material is assumed to be the same as that of the surface of the mass gainer.', '2.1.3 Treatment of the contact phase': "Mass transfer during contact, i.e., when both stars in a binary system simultaneously overflow their Roche-lobe volumes, is implemented in MESA as described in Marchant et al. (2016) and Marchant (2017). We summarize the main features of the contact scheme here.  \nDuring contact we assume that the surfaces of both stars lie on a common equipotential surface. The amount of mass transferred from one star to the other is adjusted such that the following relationship holds:  \n𝑅 2 ' Φ ' GLYPH<0> 𝑅 RL,2 𝑅 RL,2 = 𝐹 ' 𝑞GLYPH<150> 𝑅 1 ' Φ ' GLYPH<0> 𝑅 RL,1 𝑅 RL,1 ' GLYPH<149> (1)  \nwhere Φ is the equipotential surface shared by both stars, R j and R RL,j are the volume-equivalent stellar and Roche lobe radius respectively of either star (j=1,2) and 𝑞 = 𝑀 2 GLYPH<157> 𝑀 1 where 𝑀 2 is the current less massive star in the binary. The function 𝐹 ' 𝑞GLYPH<150> 𝑥 ' , where 𝑥 = 𝑅 1 ' Φ 'GLYPH<0> 𝑅 RL,1 𝑅 RL,1 , is solved by numerically integrating the equipotential volume of each star through the inner Lagrangian point L1 for different mass ratios. 𝐹 ' 𝑞GLYPH<150> 𝑥 ' is thereby approximated as 𝐹 ' 𝑞GLYPH<150> 𝑥 ' = 𝑞 0 GLYPH<149> 52 𝑥 for 𝑥 ¡ 0 such that 𝐹 ' 𝑞GLYPH<150> 𝑥 ' = 0 if 𝑥 = 0 and is equal to 𝑥 if 𝑞 = 1, i.e., both stars have equal radii when they attain equal masses.  \nThe volume equivalent radius corresponding to the outer Lagrangian point L2 (R L2,2 ) depends on the Roche lobe radius of the less massive star (R RL,2 ) and the mass ratio 𝑞 at the time considered, and is calculated as:  \n𝑅 L2,2 GLYPH<0> 𝑅 RL,2 𝑅 RL,2 = 0 GLYPH<149> 299 tan GLYPH<0> 1 GLYPH<16> 1 GLYPH<149> 84 𝑞 0 GLYPH<149> 397 GLYPH<17> GLYPH<149> (2)  \nEnergy transport in the common envelope is not accounted for during the contact phase. Note that in this subsection (and this subsection only) we follow the notation used by Marchant et al. (2016) where subscript 1 refers to the star that is the more massive one at a given time during the evolution of the binary model. This maybetheprimary (the initially more massive star) or the secondary (the initially less massive star) if the mass ratio has been reversed as a result of mass exchange. Throughout the rest of this paper we will use the common convention where 1 (2) refers to the initially more (less) massive star or the primary (secondary) star.", '2.2.1 Initial parameters': "We compute a total of 2790 binary models for each galaxy grid and begin the evolution of each model with both stars on the ZAMS (we describe this initialization in Section 2.2.2). The initial parameter range of our grids and the grid spacing ( Δ values) are:  \n- GLYPH<15> Initial mass ratio : q i GLYPH<17> M 2,i GLYPH<157> M 1,i = 0.6, 0.7, ...,1.0 with Δ q i = 0 GLYPH<149> 1 and we add an extra grid for q i = 0 GLYPH<149> 95.\n- GLYPH<15> Initial total mass : M T,i GLYPH<17> M 1,i , M 2,i = 20, 22, ...,80 M GLYPH<12> with Δ M T,i = 2 M GLYPH<12>\n- GLYPH<15> Initial period : P i = 0 GLYPH<149> 6 GLYPH<150> 0 GLYPH<149> 7 GLYPH<150> GLYPH<149> GLYPH<149> GLYPH<149> GLYPH<150> 2 GLYPH<149> 0 days with Δ P i = 0 GLYPH<149> 1 days  \nWe restrict the lower limit of the initial mass ratio to q i = 0 GLYPH<149> 6 as convergence errors become more common for lower initial mass ratios and the upper limit of the total mass to M T,i = 80 M GLYPH<12> as the observed contact binaries have total masses less than this value (except for one). Our choice to space our grid in equal steps of the initial total mass, M T GLYPH<150> i , instead of the more usual choice of the  \ninitial primary mass M 1 GLYPH<150> i is a natural one when considering contact binaries that experience near conservative mass transfer. The total mass stays approximately constant during the early evolution of these models while the masses of the individual components can change substantially as the system experiences one or more phases of mass transfer (cf., Fig 1). We note that our choice of sampling the grid evenly in total mass and mass ratio leads to effects at the edges of our grid towards the lowest and highest initial primary masses (M 1,i GLYPH<159> 12 GLYPH<149> 5 M GLYPH<12> and M 1,i ¡ 40 M GLYPH<12> , respectively) which are only possible for a limited combination of choices of the mass ratio. However, these 'edge effects' are found to be very small and do not affect our main conclusions.  \nThe range in initial orbital period is motivated to cover the systems of interest for the formation of long-lived contact systems. Systems with shorter periods are already in contact at zero age and merge immediately. Although systems with wider periods may contribute, especially for the most massive systems that start with almost equal masses, we expect that their contributions do not significantly affect our main conclusions at the metallicities of our grids. Sen et al. (in prep) also report a similar finding for binary models with initial orbital periods above 2 days.  \nOur choice of initial parameter space is different from that of Marchant et al. (2016), who were mainly interested in the parameter space in which binaries can evolve chemically homogeneously and form double black hole systems. They hence evolved very lowmetallicity massive binary models with initial primary masses between 25-502 M GLYPH<12> and initial mass ratios equal to 1 (except for their most metal-poor grid). In comparison, we explore much smaller masses with initial mass ratios from 0.6-1 in our grids, as we intend to cover the parameter range of the spectroscopic contact binaries identified in the Magellanic Clouds. In addition, our method of initializing our binary models is also different from Marchant et al. (2016), as we explain in the next section.", '2.2.2 Initialization of models': "Before initiating the binary evolution of models in our grid, we first build a separate library of starting models. For each binary model we compute single star ZAMS models of masses M 1,i and M 2,i . We define the ZAMS model as the point of evolution at which a stellar model has contracted from the pre-main sequence (pre-MS) branch to its smallest radius before expanding on the MS.  \nThis initial setup differs from that of Marchant et al. (2016) where the binary evolution is initiated by allowing MESA to find the appropriate pre-MS models for the required masses and initial composition. In their setup, the two stellar models are puffed up at the beginning of the evolution before shrinking and expanding again during their MS evolution. As a consequence, the stars may interact prior to reaching the ZAMS point by transferring mass and thereby changing the orbital period before their MS expansion actually begins. Some of the models in this setup of Marchant et al. (2016) even merge during the transition from the pre-MS to MS branch.  \nThe 'initial' masses and orbital period of the binary models presented in this paper refer to the masses of the primary and secondary stars and the orbital period of the binary system at the beginning of the ZAMS. By using carefully constructed starting models for the primary and secondary stars, our initial setup allows us to provide more direct predictions of the evolutionary trends of close binaries. We also find that more binary models avoid merging initially on the ZAMS and their overall lifetime on the MS is longer  \ncompared to the models from the setup of Marchant et al. (2016) since fewer of them come in contact at zero age.  \nFor the initial composition we follow Brott et al. (2011a) and choose a chemical mixture that is representative of stars in the SMC and LMC, with Z SMC = 0 GLYPH<149> 0021 and Z LMC = 0 GLYPH<149> 0047 which are GLYPH<25> Z GLYPH<12> GLYPH<157> 2 and GLYPH<25> Z GLYPH<12> GLYPH<157> 5 respectively, where Z GLYPH<12> = 0 GLYPH<149> 014 (Asplund et al. 2005). The individual abundances of 7 elements: H, He, C, N, O, Mg, Si and Fe are included as reported in Table 1 of Brott et al. (2011a) and the abundances of all other elements are scaled-down to the solar abundances of Asplund et al. (2005).  \nFor each binary model in a grid corresponding to a given metallicity, we need three main inputs to begin their evolution: their initial period and the ZAMS models corresponding to the individual masses, M 1,i and M 2,i . At the start of the simulation we assume the binary is tidally synchronized, hence within a short timescale (less than 0.01 Myr), the rotational periods of the stars become equal to the orbital period of the system.", '2.2.3 Outcomes': "Since we are only modelling the MS binary phase of each model, the evolution is terminated when either star leaves the MS or when the system undergoes L2 overflow. The latter is expected to result in mass loss through the outer Lagrangian point and a corresponding loss of angular momentum from the binary, leading to a common envelope phase and the eventual merger of the binary system. Our models experience one of these three final outcomes:  \n- (i) Overflows L2 initially : The stars are initially so close that the binary system experiences mass loss through the L2 Lagrangian point while the stars have barely evolved (the secondary has depleted less than 3% of its central hydrogen abundance).\n- (ii) Overflows L2 on MS : If the binary system experiences overflowthrough the L2 point post ZAMS but while still on the MS, they are denoted as 'Late L2 overflow' systems. These systems merge while still burning hydrogen in their cores.  \n(iii) Survive the MS : If one or both stars evolves past the MS, we call these 'MS survive' binaries.  \nOf the 2790 models in the LMC grid, nearly 13% faced convergence errors i.e., the binary evolution was unable to initiate correctly or terminated sooner than their proper end (one of the three outcomes above). The contribution from these non-convergent models to our synthetic populations is negligible since they fall in the parameter space where they will either experience L2 overflow at the start of the simulation or will undergo contact only on a thermal timescale ( GLYPH<159> 0 GLYPH<149> 01 Myr). Of the models that do converge, nearly 16% of them experience 'initial L2 overflow' and 63% of them 'overflow L2 on MS'; thus nearly 80% of our LMC models merge on the MS. These fractions are similar for the SMC grid of models as well.", '2.3 Modelling the population of contact binaries distributions': "Each evolutionary model is characterized by a unique combination of initial orbital period P i , initial total mass M T,i and initial mass ratio q i (and hence a choice of M 1,i ) sampled from the parameter space as described in Section 2.2.1. The birth distributions are given by:  \ndN dM 1,i GLYPH<24> M 𝛼 1,i GLYPH<150> dN dq i GLYPH<24> q i 𝜅 GLYPH<150> dN dP i GLYPH<24> P i 𝛾 GLYPH<150> (3)  \nFor the purpose of a first order calculation, we assume the initial mass function (IMF) has an exponent of 𝛼 = GLYPH<0> 2 GLYPH<149> 35 according to  \nthe Salpeter IMF (Salpeter 1955; Kroupa 2002). We also assume an Öpik law in period, i.e., flat in log P (which is 𝛾 = GLYPH<0> 1 for P), and a flat distribution in q, i.e., 𝜅 = 0. The birth weight of a given initial choice of M 1,i , P i and q i is then calculated as:  \n𝑤 M1 GLYPH<150> i = C M1 h M GLYPH<0> 1 GLYPH<149> 35 1,i i M1,i , dM1,l M1,i GLYPH<0> dM1,u (4a)  \n𝑤 Pi = C P GLYPH<2> log 10 P i GLYPH<3> Pi , dP Pi GLYPH<0> dP (4b)  \n𝑤 q i = Cq GLYPH<2> q i GLYPH<3> q i , dq q i GLYPH<0> dq GLYPH<150> (4c)  \nwhere C M1 GLYPH<150> C P , and Cq are normalization factors multiplied with their respective weights to ensure that the sums of 𝑤 M1 GLYPH<150> i , 𝑤 Pi and 𝑤 q i are each equal to 1.  \nThe mass, period, and mass ratio intervals in Eqs. 4 (dM 1 , dP, dq) define the grid cell dimensions for each considered binary evolution model in the initial parameter space.  \nFor the initial period, dP = 0 GLYPH<149> 05 days, which is half the grid spacing. For the mass ratio, we choose dq = 0 GLYPH<149> 05 for all choices of q i except for the highest mass ratios where our models are spaced more densely. We adapt the cell boundaries such that they lie exactly in between our choices of q i . In particular, we choose dq for our q i = 1 GLYPH<149> 0 grid of models such that they represent the systems between q i = 0 GLYPH<149> 975 and q i = 1. While our grid is evenly spaced in total mass (see Section 2.2.1), the spacing in initial primary mass is not even as M 1,i = M T,i GLYPH<157> ' 1 , q j ' for any combination of i and j. Therefore, we define the upper and lower limits for M 1,i as:  \ndM 1,u = ' M 1,i+1 GLYPH<0> M 1,i ' 2 M GLYPH<12> GLYPH<150> dM 1,l = ' M 1 GLYPH<0> M 1,i-1 ' 2 M GLYPH<12> GLYPH<149> (5)  \nwhere the cell size does depend on the mass ratio q j . However, with this definition, the parameter space is covered without gaps or overlaps. For the primary masses at the edge of the grids, we use: at M T,i = 20 M GLYPH<12> , dM 1,l = 0 and at M T,i = 80 M GLYPH<12> , dM 1,u = 0 while for all other primary masses the limits are calculated as in Eq. 5.  \nFinally, the statistical weight 𝑤 s of a model 's' with its particular combination of the above initial parameters, which is the convolution of the birth weight and the dimension of the cell it represents, is:  \n𝑤 s = 𝑤 M1 GLYPH<150> i GLYPH<2> 𝑤 q,i GLYPH<2> 𝑤 P,i GLYPH<149> (6)", '2.3.1 Probability distributions of contact binaries according to their initial parameters': 'One of our goals is to study the distribution of the fraction of time systems spend in contact compared to their overall MS lifetime, weighted by their birth weights. We define this fractional contact time per system in two ways:  \n𝑭 1,s = 𝑤 s 𝜏 contact,s 𝜏 MS-binary,s GLYPH<150> 𝑭 2,s = 𝑤 s 𝜏 contact,s 𝜏 MS-single,s GLYPH<150> (7)  \nwhere 𝜏 contact,s is the net time spent in contact by the system and 𝑤 s is the birth weight of the system as calculated in Eq 4. 𝜏 MS-binary,s is the lifetime of the binary model before it merges on the MS or before either of its components leaves the MS. In practice since many of these close binaries merge before they exhaust their central hydrogen, we also define 𝜏 MS-single,s , which is the MS lifetime of a star had it not been part of a binary system that is close enough to interact. We take 𝜏 MS-single,s for each binary system to be the lifetime of a single star of 1/2 the total binary mass, M T,i . The values of 𝜏 MS-single,s for each M T,i /2 in our grids, were supplied', 'Menon et al.': 'from the grid of single star models computed by Schootemeĳer, A. via private communication. The SMC set of these models are published in Schootemeĳer et al. (2019).  \nWethencompute the overall fraction of contact systems among the MS binaries in our population, 𝒇 contact/MS . Using this number we can predict the number of massive contact binaries in a given population of MS binaries. It is calculated as follows:  \n𝒇 contact/MS = " n s=1 𝑤 s 𝜏 contact,s " n s=1 𝑤 s 𝜏 MS-binary,s GLYPH<150> (8)  \nwhere s is an index that loops through all 2790 binary models in each grid.  \nA detailed explanation for how we compute the distributions shown in this paper can be found in Sections A1 and A2.', '3 RESULTS': 'In the following section, we will describe the properties of our binary evolution models, which will be mainly presented in the form of distribution functions. We start by discussing two of our models in detail, which are representative of the evolutionary paths followed by the majority of contact binary models in our study.', '3.1 Example models': 'The binary models in our grids go through at least two of three mass transfer episodes which occur in the following order: i) contact on a thermal timescale ii) semi-detached phase or second contact on a thermal timescale iii) contact on a nuclear timescale. The two example models we use to illustrate the major evolutionary channels in our grids only differ in their initial periods, but otherwise have the same initial mass ratio of q i = 0 GLYPH<149> 8 and total mass of M T,i = 26 M GLYPH<12> . \'System 1\' has an initial period of 0.8 days and \'System 2\' has an initial period of 1.0 day. As defined in Section 2.1.3, the primary is the initially more massive star and the secondary is the initially less massive star. Please refer to Figs. 1 and 2 (left panels for system 1 and right for system 2) in tandem with this section.  \nThe first episode of mass transfer in system 1 occurs soon after initializing the binary evolution. A fast Case A mass transfer occurs from the primary to the secondary for the first GLYPH<25> 0 GLYPH<149> 01 Myrs. The orbital period dips slightly during this mass transfer phase and then rises again when the mass ratio reverses, as is expected during conservative Case A mass transfer (point A, left panel set of Fig. 1). Stars with radiative envelopes typically shrink in radii as they lose mass and experience a drop in their surface luminosities and effective temperatures, while the reverse happens when they gain mass. These effects can be seen in the trajectory of the primary and secondary stars on the HR diagram in Fig. 2. Thereafter, the more massive secondary transfers mass on a nuclear timescale and the system enters contact. The orbital period decreases as the system converges towards equal masses. Both stars expand on the MS and remain in contact as a slow mass transfer continues from the secondary until it overflows the L2 Lagrangian point at about 6 Myr, after which the evolution is terminated (point B). Models that evolve like system 1 are in contact for most of their binary lifetimes before merging on the MS. Models that evolve like system 1 have the longest contact durations in both our metallicity grids and their contribution dominate the synthetic binary populations we will see in later sections.  \nSince system 2 begins initially wider than system 1, its first  \nmass transfer episode occurs at about 3 GLYPH<149> 6 Myr. The system enters contact and transfers mass on a thermal timescale (points A to B, right panel set of Fig. 1). The primary loses mass and shrinks in radius while the secondary gains mass and expands (right panel, Fig. 2). The orbital period drops and rises again in this fast masstransfer episode, and both stars detach and break contact. Both stars continue to expand and at nearly 4.3 Myr, a second mass transfer phase ensues from the now less-massive primary on a nuclear timescale, and the system enters a semi-detached configuration for nearly 3 Myr (points C to D). As the donor is the less massive star, the orbital separation widens during this mass transfer phase. The mass-transfer scheme employed by MESA during the semi-detached phase fixes the radius of the donor to its Roche lobe radius such that the donor radius remains constant during this phase. In the meantime, the secondary gains mass and continues to expand but does not overflow its Roche lobe until nearly 6.5 Myr. At this point, the third mass transfer episode ensues with both stars overflowing their Roche lobes and the system thus enters contact.  \nThe mass transfer during this contact phase occurs from the now more massive secondary to the less massive primary. Initially the reverse mass transfer rate is high, reaching about 10 GLYPH<0> 5 M GLYPH<12> yr GLYPH<0> 1 but it gradually decreases to a nuclear timescale and the mass ratio steadily approaches unity as the orbital period decreases. This type of evolution is reminiscent of the typical rapid and slow phases of Case A systems described in the literature, albeit with significant differences. The initial rapid phase lasts longer than the thermal timescale that is typically seen in standard Case A systems and the mass ratio does not (re)reverse in this system after it achieves q = 1. After nearly 3 GLYPH<149> 4 Myr in this slow contact phase, the system experiences L2 overflow (point E) with a final mass ratio very close to unity. Some models in our grids that evolve like system 2 experience L2 overflow before attaining a mass ratio of 1. In such cases, the final mass ratio of the system prior to merger is close to, but slightly different from 1.  \nIt is not only important to consider the contact phase where both stars overflow their Roche lobe volumes, but also phases where both stars are nearly in contact. As mentioned in Section 1, observationally it can be difficult to distinguish over-contact systems from those approaching contact. We assume that this \' near-contact phase \' in binaries, occurs when both stars simultaneously have R/R L GLYPH<21> 0 GLYPH<149> 9, while not yet achieving contact (R/R L GLYPH<21> 1 for both stars). Models that follow the evolution of system 2 contain an extended phase of being in near contact during the semi-detached phase (between points B and C) of their evolution. While the actual contact phase of system 2 lasts 3.4 Myr, it spends an additional 3.2 Myr in near contact. We hence introduce the definition of the \'relaxed contact phase\' , which is the sum of the near and actual contact phases of a system, during which it may be classified as a \'true" contact system. In the case of system 2, the duration of this relaxed contact phase is 6.6 Myrs.  \nThe effect of including these relaxed-contact binaries in our synthetic stellar populations significantly impacts our predictions for the observed contact binary population, as we shall see in Section 4.', '3.2 Impact of the initial mass ratio on systems 1 and 2': 'In general, the initial mass ratio (q i ) of the models in our grids determine the duration of their 𝜏 MS,binary and 𝜏 contact; the closer their q i is to 1, the longer they are in contact before merging.  \nFig. 3 shows the evolution of systems 1 and 2 but now with varying initial mass ratios. Independent of q i , models that evolve  \n˙  \nFigure 1. The evolution of the two example systems described in Section 3.1, with initial masses of 14.4 and 11.6 M GLYPH<12> for the primary (solid blue lines) and secondary (dashed blue lines) respectively. System 1 (left panels) begins with Pi = 0 GLYPH<149> 8 days and System 2 (right panels) begins with Pi = 1 GLYPH<149> 0 days. In the topmost panels, the left y-axis represents the period evolution and the right y-axis the mass ratio evolution. The black dotted horizontal line in the topmost panels is a reference for q = 1 GLYPH<149> 0 while in the third panel, the dotted line indicates when both stars have R/RL ¡ = 0 GLYPH<149> 9. The letters in each panel correspond to the respective phases of each system, as described in Section 3.1. The true contact phases (where both stars have R/RL ¡ = 1) are shaded in grey (nuclear timescale) and red (thermal timescale).  \n<!-- image -->  \nlike system 1 enter contact almost immediately after ZAMS and attain a mass ratio of unity within about 1 Myr, staying in contact until the end of their MS evolution (see also Marchant et al. 2016).  \nIn contrast, the final mass ratio of models that evolve like system 2 does depend on their initial mass ratio. The further their initial mass ratios are further away from 1, the shorter the duration of their contact phase and the time they spend as equal-mass contact binaries. While all the models asymptotically approach q = 1, the majority of models with q i = 0 GLYPH<149> 6 and 0.7 experience L2 overflow prior to attaining q = 1. With decreasing q i the duration of the semi-  \ndetached phase, which is responsible for the second increase in mass ratio of the models, becomes smaller and altogether disappears in the q i = 0 GLYPH<149> 6 model. This is because the system widens considerably after the first thermal-timescale mass transfer episode and remains detached for about 2 Myrs. The next mass transfer phase occurs from the massive secondary on a thermal timescale thereby reversing the mass ratio from q = 1 GLYPH<149> 65 to 1.2 and the system enters its second contact phase. Thereafter, the mass transfer rate slows down to the nuclear timescale of the expanding secondary and undergoes L2 overflow before attaining a mass ratio of unity.', '8 Menon et al.': '<!-- image -->  \nFigure 2. HR diagram of the two systems from Fig. 1, including the phases of evolution marked by their corresponding letters (see text in Section 3.1). Solid line represents the primary star and dashed line represents the secondary star. The grey region denotes the nuclear timescale contact phase and the red region denotes the thermal timescale.  \n<!-- image -->  \n<!-- image -->  \nFigure 3. The impact of the initial mass ratio on the evolution of binary systems that come in contact. We show the mass ratio q as a function of time for systems 1 and 2 from Fig. 1, for a fixed total mass MT,i = 26 M GLYPH<12> and respective initial periods (0.8 days and 1.0 day) but with different initial mass ratios qi = 0 GLYPH<149> 6 GLYPH<150> 0 GLYPH<149> 7 GLYPH<150> 0 GLYPH<149> 8 GLYPH<150> 0 GLYPH<149> 9. The colored shaded regions indicates the contact phase corresponding to each qi. The dashed black line indicates q = 1.  \n<!-- image -->  \nFigure 4. Reverse cumulative distribution of three timescales of all models in the SMC grid (left plot) and LMC grid (right plot), of the contact phase 𝜏 contact (hatched purple histogram), 𝜏 MS,single (orange filled histogram) and 𝜏 MS,binary (green histogram). See Section 2.3.1 for their definitions. Only those models which have 𝜏 MS,binary GLYPH<21> 0 GLYPH<149> 1 Myr are included in these histograms.  \n<!-- image -->  \nIn both types of evolution, the net contact duration decreases as the initial mass ratio decreases. This is because the expansion of the primary star drives the evolution of the system; the more massive the primary is initially, the sooner the system attains contact  \nFigure 5. Weighted histogram of fractional contact duration for all the models in the SMC (grey striped histogram) and LMC (purple filled histogram), as a fraction of the equivalent single-star MS lifetime (left plot) and the binary model MS lifetime (right plot). Each model is weighted by its birth parameters. As in Fig. 4, we do not include models with 𝜏 MS,binary GLYPH<159> 0 GLYPH<149> 1 Myr.  \n<!-- image -->  \nand more rapid its overall evolution. In the case of models that evolve like system 2, the dip in orbital period corresponding to the (re)reversal of the mass ratio (during the second rapid mass transfer phase) becomes more pronounced as the initial mass ratio decreases, causing the system to merge sooner as well.  \nWe summarize our findings thus far: a mass ratio of one is the equilibrium configuration that all our binary models proceed towards when they come in contact over a nuclear timescale. As their initial mass ratio increases, the more time the binary model spends in contact and the likelier it is to attain a final mass ratio of one before merging.', '3.3 A census of the contact binary models in our grids: their life-span and fates.': "3.3.1 The duration of the contact phase in our models  \nWe briefly examine the duration of the contact and MS phases of the binary models in our grid. Figure 4 shows the contact duration for the models as an unweighted cumulative histogram, along with their MS life time before merging ( 𝜏 MS, binary ) and their equivalent lifetime had they been single stars ( 𝜏 MS GLYPH<150> single ), which is the MS lifetime of a star with half the total mass of the binary system (cf., Section 2.3.1). For the sake of clarity, we exclude systems that merge within 0.1 Myr ( 𝜏 MS,binary GLYPH<21> 0 GLYPH<149> 1 Myr).  \nWe see that about 65% of our models experience contact for less than 1 Myr. These systems are born in contact or very close to contact and merge soon after. The net contact duration of a system also decreases when its total mass increases, owing to their shorter overall nuclear timescale, or with initial mass ratios further from 1 (cf., Section 3.2). The longest contact durations are 10 GLYPH<149> 1 Myr and 10 GLYPH<149> 7 Myr in the LMC and SMC grid respectively, which belong to the models with the lowest initial total mass and an initial mass ratio of one.  \nFigure 5 shows the distribution of the fractional contact durations, normalized to their MS durations and weighted by their birth parameters. On average, the weighted contact lifetimes of our models is about 45% of their MS binary lifetimes, for both the LMC and SMC (right plot in Fig. 5). We also see that the weighted distribution of 𝜏 contact GLYPH<157> 𝜏 MS,binary has two peaks, at 0-10% and at 90-100%. While the first peak corresponds to the similar peak in the left plot of Fig. 5 (between 0 and 0.1), the second peak is produced by systems that remain in contact for most of their MS binary lifetimes similar to system 1 or merge shortly after being born in contact (but with 𝜏 MS,binary ¡ 0 GLYPH<149> 1 Myr). To discern the contributions of the latter type of models that merge soon after ZAMS, the normalisation to equivalent single star lifetimes, 𝜏 MS GLYPH<150> single , gives a better idea about the expected number of contact binaries in a population of mixed single stars and binaries (see Section 6 below). From Fig. 5 we find that our models are in contact for at most 50% of their single-star MS lifetime.  \nWe further note that the distributions in Figs. 4 and 5 for our SMCand LMC models are remarkably similar, despite a difference of a factor of 2.5 in their metallicity. This is an indication that metallicity is not a primary factor of importance, at least in the parameter range studied here.  \nFigure 6 shows two slices of our three dimensional LMC grid with the contact duration as a function of the initial period (y-axis) and initial total mass (x-axis) for the models with the most extreme mass ratios: q i = 0 GLYPH<149> 6 and q i = 0 GLYPH<149> 95. We do not include the q i = 1 models in this analysis as they are special and will be discussed later in this section.  \nSystems that experience L2 overflow at initialization (yellow squares in Fig. 6) are expected to merge quickly and not contribute to the observable population of contact binaries. The fraction of models with convergence issues (grey shaded squares) increases for total masses more than 40 M GLYPH<12> and periods shorter than 1.3 days,  \nwhere we typically also find L2 overflow at the start of the evolution. This fraction increases for our more extreme initial mass ratios, where we experience problems at wider periods for the higher total masses.  \nModels that evolve similar to our example system 1 (blue hatched squares) become less common as the total mass increases (for a fixed period) and as the period increases (for a fixed total mass). For the lowest total mass systems in our grid, models follow a system 1 like evolution for initial periods up to 0.8 days for the least massive models in our grid, and initial periods between 1.1 and 1.5 days for the most massive systems in our grid. These models spend at least 70% of their MS binary lifetimes as contact systems and have the longest contact durations in our grid. They also attain equal masses before merging during core hydrogen burning.  \nModels that evolve in a similar way as our example system 2 are very common, and typically come from models with initial periods larger than 0.8 days and become more common as q i decreases. These models evolve through a semi-detached phase and may merge while they still have unequal masses.  \nModels in which at least one star evolves past the MS (the 'MS survive systems', marked as green stars in Fig. 6) have orbital periods larger than 1.5 days and become increasingly common as the total masses and initial separation increases. For the more massive systems, the initial period at which the binary will survive the MS without merging lies beyond the edge of our grid. These models may still merge when the star that leaves the main sequence swells up during H-shell burning. Since however, the post-MS phase of a star is at least an order shorter in duration than its MS lifetime, we do not expect models undergoing contact after their MS phase to contribute significantly to the observed population of contact binaries. We also report that none of our models evolve completely chemically homogeneously.  \nOur models with q i = 1 are special in that, both stars have exactly the same mass. They hence evolve simultaneously and fill their Roche lobes at the same time. They do not exchange mass but keep expanding on their nuclear timescale until L2 overflow occurs. Hence the duration of the contact phase from these models is longer than their corresponding models with unequal initial mass ratios (cf., Fig. A1). The contribution of binaries with stellar masses that are exactly equal at birth is not well known, hence we reduce the importance of the q i = 1 models by our choice of the integration boundaries (cf., Section 2.3)  \nGenerally we see that systems with the longest contact duration (those with the darkest shade of purple in Figs. 6 and A1) are concentrated at initial orbital periods at 0.8 days and total masses of 20 M GLYPH<12> (the lowest in our grid) across all initial mass ratios. The model with the longest contact duration in our LMC grid has P i = 0 GLYPH<149> 8 days, M T,i = 20 M GLYPH<12> and q i = 1 and spends about 10.1 Myr in contact. The same model in the SMC grid spends about 10.7 Myr in contact.", '3.3.2 The impact of initial binary parameters on the overall grid population': "In Fig. 7 we fold the contact durations for all systems in each grid (Fig. A1) along with the birth weight of each system, to give us the likelihood of finding a contact binary in our synthetic population with a given initial period and mass ratio (cf., Eq. A1 in Section 2.3.1). The more purple the background color for each pixel is, the more likely that combination of P i and M T,i produces a contact binary. Between our LMC and SMC populations the main difference is the curves of 'Initial L2 overflow' below which all models merge  \nFigure 7. 2D probability distribution showing the contributions of systems to a synthetic contact binary population as function of their initial period (Pinitial) and initial total mass (MT,initial), across all initial mass ratios. The background color of each pixel is the normalized probability of a binary with a particular initial orbital period and total mass contributing to the contact binary population (cf., Section A1). The yellow rectangular box emphasizes the pixel which provides the largest contribution. Turquoise star symbols indicate pixels which contain binary models which avoid merging during the MS evolution. The area below the dotted black line indicates those systems that experience initial L2OF for an initial mass ratio of qi = 1 GLYPH<149> 0.  \n<!-- image -->  \nglyph[circledot]  \n<!-- image -->  \nglyph[circledot]  \nFigure 6. Contact duration ( 𝜏 contact) of our LMC models as function of their initial orbital period (Pinitial) and initial total mass (MT,initial) for two initial mass ratios, q i = 0 GLYPH<149> 6 (left) and q i = 0 GLYPH<149> 95 (right). The background colors represent the contact duration of each system in the grid; the more purple a pixel is, the longer the contact phase of the corresponding system is. The parameter space colored in yellow undergo Initial L2 overflow while the one colored in grey experienced convergence errors. Grey pixels in the lower half of the diagram are also expected to experience Initial L2 overflow. Pixels with blue hatching represent models that evolve as system 1, while those are that are neither hatched nor shaded in grey or yellow, represent models that evolve like system 2. Systems with white backgrounds spend less than 0.01 Myr in contact. Green star symbols represent 'MS survive' models, i.e., the system avoids merging during core hydrogen burning. All other systems experience L2 overflow and merge during their MS evolution.  \n<!-- image -->  \nglyph[circledot]  \nglyph[circledot]  \nat the beginning of their binary evolution. We show these curves in Fig. 7 for the q i = 1 models. Owing to their higher metallicity, stars of the same initial mass are bigger in radii in the LMC compared to the SMC, due to which more binaries merge initially at larger separations in the LMC than the SMC. Apart from the shift in this boundary of 'Initial L2 overflow' we find remarkably similar results between the SMC and LMC grids (see also our comparison between these grids in Fig. 4 and 5).  \nThe general trends already observed in the panels in Fig. 6 are still visible. Contact binaries originate predominantly from the longer lived lower mass systems with short orbital periods. Models that have total masses of 20 GLYPH<0> 24 M GLYPH<12> dominate the population of contact binaries across all initial periods and mass ratios considered in our study. This is because along with their actual contact lifetimes, these systems are also weighted higher by the IMF, although it should be noted that we have used a log scale for the color bar in Fig. 7 and this may enhance features that could only be seen in a population that contains hundreds of contact binaries.  \nSystems which avoid merging during the MS evolution occur for a larger range of initial orbital periods and total masses in the SMC grid compared to the LMC, as their overall radii are smaller than their LMC counterparts. The large majority of these systems have initial orbital periods larger than 1.3 days for the SMC grid and 1.5 days for the LMC grid. We also see from Fig. 6 that the binaries which avoid merging on the MS only come in contact over a thermal timescale or do not come in contact at all. The implication for the observed contact binaries which are in contact over nuclear timescales thus is that, they will merge before either of the two stars in the binary finishes core hydrogen burning.", '4.1 The observed sample': "We have collected the data of 26 spectroscopic massive contact binaries from the Magellanic Clouds (Table 1) and the Galaxy (Table 2), for which the orbital and stellar parameters are available. These include 2 systems from the SMC, 7 from the LMC and 17 from the MW. The majority of these are O+O contact binaries, with a few B+B and O+B systems. Most of these systems are located close to ZAMS and a few are located much further along the MS (Fig. 8). A few of the Galactic systems have tertiary companions (or systems) reported as well. However, we do not model the effect of tertiary systems on the evolution of contact binaries as this is beyond the scope of our work.  \nOf the 7 O-type spectroscopic contact binaries in the LMC, 4 are from the VFTS-TMBM sample. This is because the VFTS sample is largely biased (and complete) in O stars (Sana et al. 2013; Almeida et al. 2017), which in turn may be because 30 Dor itself has a larger concentration of very massive stars than the average LMC population (Schneider et al. 2018a,b). In addition to the VFTS sample, 29 'extremely-blue' contact binaries have been identified from the MAssive Compact Halo Objects (MACHO) survey data (Rucinski 1999) , all located close to the Zero Age Main Sequence (ZAMS), with periods of 0.45-1.3 days and a few outliers with periods up to 3.5 days. From the Optical Gravitational Lensing Experiment (OGLE)-III catalogue, 28 early-type contact binaries with periods between 0.5 and 1.6 days have been reported (Pawlak 2016). There is an overlap between the two data sets and roughly, we estimate there to be 30 unique O & B-type contact binaries together from the MACHO and OGLE database for the LMC. Unfortunately,  \nsince these are photometric binaries, the stellar parameters of individual systems are not available.  \nMost contact systems are eclipsing binaries that are distinguished by their smooth, nearly-sinusoidal light curves which do not show plateaus between eclipses (Lorenzo et al. 2017), indicating that both stars have near-equal brightness and show signs of tidal distortion (Hilditch et al. 1998). By modelling the light curve and with the radial velocity measurements, one can discern the geometry of the binary and the degree by which the stars overflow their respective Roche lobe volumes- their 'fill-out factors'.  \nIt is in the uncertainty of these fill-out factors that the exact nature of a system becomes unclear; whether it is truly a contact binary, or if it is simply approaching contact. The main contributor to this uncertainty is the inclination of the system- the further away from 90 GLYPH<14> the inclination is, the larger the error bars are on the measurements of the stellar parameters.  \nThefill-out factors themselves are reported differently, depending on how they are calculated. One method to calculate the fill-out factor ( 𝑓 Ω ) is by comparing the difference between the surface potential ( Ω ) of the stars (assumed to be identical for either star) with the potentials at the inner and outer Lagrangian points of the system (Mochnacki & Doughty 1972). A value of 𝑓 Ω = 0 indicates a system exactly in contact and 𝑓 Ω ¡ 0 indicates an over-contact system (Lorenz et al. 1999). By this definition, V606 Cen which has 𝑓 Ω = 0 GLYPH<149> 01 GLYPH<0> 0 GLYPH<149> 04 is considered to be marginally in contact (Lorenz et al. 1999), V729 Cyg with 𝑓 Ω = 0 GLYPH<149> 17 GLYPH<0> 0 GLYPH<149> 22 as a binary in moderate contact (Yas,arsoy & Yakut 2014), and V701 Sco, GU Mon and CT Tau with 𝑓 Ω = 0 GLYPH<149> 55 GLYPH<150> 0 GLYPH<149> 70 and 0.99 respectively, as systems in deep contact (Yang et al. 2019).  \nA second approach, is to find the (over)filling factor f L , by comparing the average volume of each star to its respective Roche lobe volume, i.e., 𝑓 L = ' 𝑅 mean GLYPH<157> 𝑅 RL ' 3 . A value of 𝑓 L = 1 thus indicates a system exactly in contact. One of the hottest and most massive contact binaries known, VFTS 352, has 𝑓 L = 1 GLYPH<149> 29 for both stars (Almeida et al. 2015), thus assigning its status as a binary in deep contact. Other systems can also be classified as over-contact binaries based on this definition, such as V382 Cyg ( 𝑓 L = 1 GLYPH<149> 1 for both stars), OGLE SMC-SC10 108086 ( 𝑓 L = 1 GLYPH<149> 7 for both stars) and TU Mus ( 𝑓 L = 1 GLYPH<149> 2 GLYPH<150> 1 GLYPH<149> 3) (these 𝑓 L values are as reported by Almeida et al. 2015). Despite 𝑓 L being higher than 1 for V382 Cyg, Martins et al. (2017) who derived the parameters and geometry of V382 Cyg, report it as a system in which both stars are barely filling their Roche lobes. Penny et al. (2008) who studied TU Mus, indicate it as a system that is either approaching contact or is in marginal contact. The status of another Galactic contact binary, LY Aur, is also uncertain, with Mayer et al. (2013) indicating that it may even be a semi detached system. SX Aur is reported as a nearcontact binary by Öztürk et al. (2014) who calculated a value of 𝑓 L = 1 GLYPH<149> 02 for both its stars. Hence although the system LSS 3074 is reported as an over-contact binary by Raucq et al. (2017), given its 𝑓 L = 1 GLYPH<149> 008 GLYPH<6> = 0 GLYPH<149> 01, we suspect this may also be a near-contact binary. Other reported near-contact binaries are: OGLE SMC-SC9 175323 (Harries et al. 1997), LZ Cen (Vaz et al. 1995), XZ Cep (Martins et al. 2017), V348 Car (Hilditch & Bell 1987) and BAT99 126 (Janssens et al. 2021).  \nFor some systems, the fill-out factors are not mentioned in the works that study them; these are MACHO*05:34:41.3 GLYPH<6> 69:31:39 (Ostrov 2001), VFTS 066 (Mahy et al. 2020a,b), HD 64315 B (Lorenzo et al. 2014) and MY Cam (Lorenzo et al. 2014). As their degree of contact is not reported, we assign these systems a 'C/NC' status (contact/near-contact) in Tables 1 and 2, along with V382 Cyg, TU Mus and LY Aur in the Milky Way and, VFTS 217 and  \n<!-- image -->  \nFigure 8. The Hertzsprung-Russel diagram of contact binaries from the Magellanic Clouds (left) and Milky Way (right). Where not directly reported for the system, we calculate the log L/L GLYPH<12> with the Teff and Rmean values in Tables 1 and 2. For reference, we also plot the evolutionary tracks of rotating (solid lines) and non-rotating (dotted lines) single-star models, with the LMC and solar metallicities, and with initial masses from 10...40 M GLYPH<12> . We compute these models until the end of their MS phase and with the same physics assumptions as in Section 2.1. For the rotating models, we choose vini = 330 km s GLYPH<0> 1 , appropriate for the typical rotation rate expected for tidally locked stars in contact binaries.  \n<!-- image -->  \nTable 1. Observed parameters for O&B type massive contact systems in the LMC and SMC, including the orbital period (Psystem), mass ratio qsystem, dynamical mass (M ★ ), mean radius (Rmean) and effective temperature Teff of individual components. We also provide the status known for the system: 'C' for a confirmed contact or over-contact system, 'NC' for a confirmed near-contact system and 'C/NC' implies the status is not explicitly mentioned in the reference or it could be in either configuration due to the error bars on its Roche-lobe filling factors (cf., Section 4.1 for more details). 'P' and 'S' stand for the primary (more massive) and secondary star of the system. The superscript ' T ' denotes the confirmed presence of a tertiary system. The system MACHO CB stands for MACHO*05:34:41.3 GLYPH<6> 69:31:39. Error bars are absent for some primary Teff values as they were fixed to calculate the stellar parameters of the system.  \nVFTS563 in the LMC, which are reported as systems with 'uncertain configurations' (Mahy et al. 2020b), but which have large enough errors on their Rmean/R L values to allow them to be contact systems as well. The same holds for VFTS 066 which is reported as a contact binary from its light curve, but again has a large range in Rmean/R L to also allow it to be a near-contact binary (Mahy et al. 2020a,b). The fill-out factors of contact binaries from the OGLE+MACHO database are also unknown.  \nFrom our literature survey thus, the list of confirmed massive over-contact binaries include: V701 Sco, CT Tau, GU Mon, V729 Cyg, VFTS 352 and OGLE SMC-SC 108086. Systems reported  \nwith equal masses (q GLYPH<25> 1) are V701 Sco, CT Tau, GU Mon, HD 64315, V348 Car and VFTS 352, while systems with near-equal masses within error bars (q ' 0 GLYPH<149> 9) are OGLE SMC SC10-108086, BH Cen and LZ Cen.", '4.2 Prediction for orbital parameters in the observed population of massive contact binaries': 'In this section we present the results of our 2D probability distributions for the currently observed parameters of contact binaries for the LMCandSMC,assumingthatstarformationoccursat a constant  \nglyph[circledot]  \nFigure 9. Probability distribution of observed period and total mass (Pcontact-MT,contact) and of observed period and mass ratio (Pcontact-qT,contact), for our LMCpopulation of contact binaries (top plots) and relaxed-contact binaries. The background colour represents the probability of finding a system with a given combination of orbital parameters (cf., Section A2). The 1D shaded histograms on top and to the right are the corresponding projections of the parameters plotted in the 2D histograms. Over-plotted in solid dark lines in the 1D histograms are the corresponding distributions from the SMC grid. Observed contact systems are also overplotted; red markers for LMC, purple marker for SMC (Table 1) and grey markers for MW systems (Table 2). The observed period range of contact binaries from the OGLE, MACHO, and VFTS samples are marked in the Pcontact histogram (blue, top panel). The grey arrows in the plots indicate the position of the two MW binaries with Pcontact ¡ 5 days.  \n<!-- image -->', '14 Menon et al.': "Table 2. Observed parameters of early-type massive contact binaries in the MW. Table headings are the same as Table 1 except for the last column where we also include other papers which have studied a given system. Additionally in the last column, the superscript ' T ' denotes the confirmed/possible presence of a tertiary companion or system and ' A ' denotes the presence of an accretion disk. Where possible we report the dynamical mass of the components else we report the evolutionary or spectroscopic masses. Systems marked with a y symbol are much further evolved from ZAMS. Finally, error bars on all parameters are only mentioned up to the first significant digit.  \nrate. Figure 9 shows the probability of finding a contact system in the LMC with a particular combination of observed period, Pcontact and total mass, M T,contact , or with a particular period and mass ratio qcontact (cf., Section A2 for the construction of these distributions). The upper two plots consider only those models which are strictly in contact (where both stars have R/R L GLYPH<21> 1), while the lower two plots include all models in relaxed contact, where both stars have R/R L GLYPH<21> 0 GLYPH<149> 9 (cf., Section 3.1). We also compare our predicted distributions with the available data of contact binaries in the Magellanic Clouds and the Galaxy (Tables 1 and 2). The distributions of the above parameters for the SMC are very similar to those of the LMC, as we can anticipate from the 2D histograms in Fig. 7. Therefore the discussion of our results for the LMC also applies to the SMC. We also include the Galactic contact binaries in our discussion for the same reason, although of course, a dedicated solar-metallicity grid would be desirable for a proper analysis of these binaries.  \nAs the contact duration of a binary is strongly dependent on its initial mass, we need to adapt our synthetic population to suit the observed populations which are dominated by very massive systems. Hence, we also compute Pcontact-qcontact distributions for  \nthe LMC and SMC grids by applying two different initial mass cuts on the synthetic contact binary population, i.e., with M T,i GLYPH<21> 30 M GLYPH<12> and M T,i GLYPH<21> 50 M GLYPH<12> . These are shown in Fig. A2.  \nIn both the top and bottom panels of Fig. 9, we see a tight correlation in the predicted distribution between M T,contact and Pcontact; the more massive a contact binary is, the larger its current orbital period will be. Another (anti) correlation can also be seen between the Pcontact and qcontact; wider systems are expected to have smaller mass ratios.  \nThe average orbital period expected for (over-)contact binaries is less than 1 day in the LMC and SMC, with a maximum of GLYPH<25> 1 GLYPH<149> 5 days for the lowest mass contact binaries. The mass ratio distribution from this population predicts that nearly 90% of overcontact systems must have mass ratios GLYPH<21> 0 GLYPH<149> 9. Not surprisingly, the most likely total mass of contact systems is close to the lower total mass boundary of our grids, due to the IMF and lifetime effects and also because the contact duration increases for smaller initial mass and initial period. The most likely orbital periods are found strongly skewed towards the smallest period values with Pcontact in the range  \n0.6 GLYPH<149> GLYPH<149> GLYPH<149> 1 day for both metallicities, and a peak slightly lower for the SMC (at 0.7 days) than the LMC (at 0.8 days).  \nWe see that the over-contact SMC system OGLE SMC SC10108086, falls well within the probable parameter space of overcontact binaries. Among the LMC contact binaries, the lowest mass system VFTS 066, has an orbital period which is about 30% larger than the most likely period at this total mass. While this could be significant, a comparison with the mass ratio distribution (right panel) shows that for its mass ratio of about 0.52, a larger period than the most likely period for this total mass is actually expected. The other five systems, lie in fact all very close to the most likely orbital period for a given system mass, which is especially clarified by Fig. A2. This even seems to hold for VFTS 217, whose mass exceeds that of our most massive models; however, this system is most likely to be a near-contact binary. We conclude that while statistics is sparse, our models reproduce the observed periods of the Magellanic Cloud contact binaries well. We also note a strong agreement between the Galactic over-contact binaries and our predictions.  \nLooking at the mass ratios of the observed contact binaries in Fig. 9, we see that indeed the confirmed over-contact binaries (VFTS352 and SMC SC10-108086) in the Magellanic Clouds have qcontact = 0 GLYPH<149> 9 GLYPH<0> 1 within their error bars, as do the three confirmed over-contact systems in the Galaxy, along with two others whose status is unclear. However, given how skewed the mass ratio distribution is towards equal masses, we would expect the majority of observed systems to cluster in the qcontact GLYPH<21> 0 GLYPH<149> 9 GLYPH<0> 1 space. This mismatch is curious, since, considering the mass dependence of the period distribution, most of the observed contact binaries follow the predicted tight correlation between orbital period and mass ratio rather closely. This is especially evident in Fig. 9 systems with M T,contact GLYPH<20> 40 M GLYPH<12> in the systems from the Magellanic Clouds as well as from the Galaxy, and in Fig. A2 for the more massive ones.  \nA more appropriate comparison for the observed systems would be with the predictions from our mixed population of nearcontact and over-contact binaries, shown in the bottom panels of Fig. 9. By including the near-contact binaries (those in which both stars have 0 GLYPH<149> 9 GLYPH<20> R GLYPH<157> R RL GLYPH<159> 1), the range of qcontact spreads down to 0.4 and that of Pcontact goes up to 2 days. The majority of the observed contact binaries are well covered in this Pcontact-qcontact space. The near-contact binaries OGLE SMC SC9-175323 and VFTS217, however, still remain outliers in this space, along with some of the Galactic systems. Hence while the likelihood for finding unequal-mass contact/near-contact binaries increases considerably, we see that, while the cut at 90% of the Roche-lobe radius is somewhat arbitrary, this effect may alone be insufficient to explain the observed q-distribution.  \nFrom the M T,contact -Pcontact diagram in Fig. 9, we see that many of the observed systems are heavier than what we expect from our synthetic population predictions. We attribute this to the fact that TMBM focuses on O-star binaries (Sana et al. 2013; Almeida et al. 2017), while in our grid, many models belong to the lowermass regime of early B stars. Hence another relevant factor to consider while comparing with observations, is the high mass bias in the observed sample. A comparison of Figs. 9 and A2 shows that the predicted preference for q contact GLYPH<21> 0 GLYPH<149> 9 decreases significantly for higher system masses-the preference for high mass ratios (qcontact GLYPH<21> 0 GLYPH<149> 9) is reduced from GLYPH<24> 90% for our complete grid to GLYPH<25> 60% for systems with initial masses above 50 M GLYPH<12> . We also note that whereas the distributions shown in Fig. 9 show very little metallicity dependence, this is different in Fig. A2 where the predicted LMCperiod distributions appear to be shifted to significantly larger  \nperiods, compared to that of the SMC for the highest mass cut as seen in Fig. A2.  \nWhile we did not model binaries with Galactic or solar metallicity, the similarity of our SMC and LMC results for the lower half of the considered mass range suggests that metallicity effects do not play a major role for system masses below GLYPH<24> 45 M GLYPH<12> . However, for more massive systems, differences become more prominent (cf., Fig. A2), and could even be larger when compared to the rather metal-rich MW.  \nWhile the LMC and SMC spectroscopic contact binaries have periods larger than 1.1 days our predicted period for contact binaries can be as low as 0.6 days. In this regard, we also consider the orbital period range of massive contact binaries from the MACHO (Rucinski 1999) and OGLE-III (Pawlak 2016) surveys, which also include near-contact binaries. Almost half of the massive contact binaries in these samples have periods GLYPH<159> = 1 day, which agrees with the predictions of our distributions. It is important to note that the smaller periods in these samples compared to the VFTS-TMBM sample may reflect the fact that the OGLE binaries are not restricted to O stars alone, and could therefore mirror the early B star systems in our model grid.  \nWesee that, in comparison with the Magellanic Cloud systems the MW contact binaries are much more scattered in the observed parameter space. In particular, we find 6 contact binaries with orbital periods in excess of 2 days (with a maximum value of 6.6 days), whereas the largest period in the Magellanic Cloud systems is nearly 1.89 days. The more extreme outliers in mass ratio and period in the MW sample, such as LY Aur, XZ Cep, V348 Car and V729 Cyg, are much further evolved from the ZAMS compared to those with periods close to or less than 2 days (Fig. 8). Notably, three of these outliers except for XZ Cep, have a total mass of more than 40 M GLYPH<12> . This may be because XZ Cep, although stated as a contact binary in Martins et al. (2017), is in fact a semi-detached system (Harries et al. 1997, private communication, Laurent Mahy).  \nWe speculate about possible reasons for this in Section 5. On the other hand, the observed MW systems with orbital periods below GLYPH<24> 1.5 days follow the predicted distributions quite closely, including a concentration of four of them with mass ratios of almost one.  \nFrom the distributions in Figs. 6 and 7, and in more detail from Fig. A1, we can obtain an idea about the completeness of our synthetic contact binary population. We see that in Fig. A1, towards the shortest orbital periods, our grid is certainly complete. For larger periods, we see that the darker purple shading remains restricted to initial orbital periods below GLYPH<24> 1 GLYPH<149> 5 days, while our upper period limit is 2 days. Only for q i = 1 we see contact times exceeding 1 Myr at initial periods of 2 days. However, at q i = 1, the nuclear timescale contact regime is limited by the maximum orbital period for Case A mass transfer, and below this limit the contact duration decreases stronger than linear with increasing period, reflecting the stellar radius evolution. Therefore, while we are missing a fraction of contact binaries starting with q i ' 1, this fraction is quite small. On the other hand, we may also be overestimating the number of contact binaries with q i ' 1, since our q i = 1 models have a statistical weight in our plots represented by the interval » 1 GLYPH<150> 0 GLYPH<149> 975 … , but perhaps the contact timescales depicted in the q i = 1 plot of Fig. A1 are representative of a smaller q i -interval, which may explain the apparent overabundance of symmetric contact binaries in our synthetic population.  \nConcerning contact binaries that may originate from smaller initial mass ratios than the lowest value in our grid, Fig. A1 shows that the dark purple island indicating the longest contact times becomes much paler towards the lowest considered initial mass ratios.", '16': 'We may miss a fraction of contact systems which originate from more extreme initial mass ratios than q i = 0 GLYPH<149> 6, but expect their integrated contact durations to be relatively small, compared to most other contact binaries in the population we model. We also see in Fig. A1 that with initial mass ratios further from one, the plots contain increasingly large areas in which the MESA models did not converge. Again, we expect the contact durations of each of those models to be relatively small.  \nThere may, however, be local regions of the predicted parameter space in our grid is incomplete, either due to numerical instability or the extent of the initial conditions. This concerns in particular the models that fail to converge with high total masses and initial mass ratios further from one, as shown in the panels of q i = 0 GLYPH<149> 6 and 0.7 of Fig. A1. The missing models in these regions may be the reason why the largest predicted orbital periods of contact binaries with M T GLYPH<150> i GLYPH<21> 50M GLYPH<12> , in the top-right panel of Fig. A2, is about 3.5 days. This is smaller than the largest orbital period predicted from our lower-mass contact binary models of about 5 days, as seen in the top-right panel of Fig. 9, which conflicts with the expectation that more massive contact binaries would generally have longer orbital periods for the same mass ratio. It is also worth noting systems with q i = 1 GLYPH<149> 0 could also contribute to even longer orbital periods, however, the probability of find binary systems that begin their evolution with with components of exactly the same mass is hard to estimate.  \nOverall, we argue that our coverage of contact binaries in initial period and mass ratio is sufficient to include the majority of the expected massive contact binary populations with SMC and LMC composition, for total system masses below 80 M GLYPH<12> .', '4.3 The number of O-type contact systems in the Magellanic Clouds': "We estimate the number of contact binaries in 30 Dor using our distributions, and thereby their numbers in the LMC and SMC. We only calculate the number of O+O type systems in 30 Dor as the VFTS-TMBM data with which we compare our results, is designed to have minimal biases for O-type stars, i.e., those with masses GLYPH<21> 15 M GLYPH<12> (Evans et al. 2011; Sana et al. 2013).  \nThe total number of O stars in the VFTS sample is N O-stars = 360 and the binaries among them span a range in period of 1.1 GLYPH<149> GLYPH<149> GLYPH<149> 3162 days, and in mass ratio of 0.1 GLYPH<149> GLYPH<149> GLYPH<149> 1.0 (Sana et al. 2013). The maximum mass of an early O-type star in the sample is 90 M GLYPH<12> (VFTS 061, Ramírez-Agudelo et al. 2017), thus lending us a mass range of 15 GLYPH<0> 90 M GLYPH<12> for O-type stars in the VFTS sample. The intrinsic binary fraction (f bin ) inferred for O-type stars in 30 Dor is 0.51 GLYPH<6> 0.04 (Sana et al. 2013), thus lending us a range in f bin of GLYPH<25> 0 GLYPH<149> 47 GLYPH<0> 0 GLYPH<149> 55. The orbital period range of the VFTS binaries is from 1.1 days to to 3162 days, while that of the mass ratio is from 0.1 to 1.0.  \nAs we are interested only in O+O-type binaries, we select only those models from our grid which have both M 1,i and M 2,i as 15 M GLYPH<12> . The fraction of O+O-type contact systems among MS binaries from our grid yields f contact/MS = 0 GLYPH<149> 31 (Eq. 8) for the strict definition of contact and f contact/MS = 0 GLYPH<149> 50 for the relaxed contact definition.  \nTo be consistent with our calculation of probability distributions (as in Section 2.3), we assume a Salpeter distribution for the IMF and a flat distribution in log period and also in mass ratio, to find the overlap between the initial parameter space coverage of our grid and that of the VFTS binaries. We calculate this by finding the fraction of the total birth weight of our grid, W grid , compared to that of the VFTS sample W VFTS . This fraction is W grid /W VFTS = 0 GLYPH<149> 03.  \nThe number of contact binaries in 30 Dor using the VFTS parameters 𝑁 contact, VFTS can be estimated as:  \n𝑁 contact, VFTS GLYPH<24> GLYPH<18> f contact/MS 0 GLYPH<149> 31 GLYPH<149> GLYPH<149> GLYPH<149> 0 GLYPH<149> 50 GLYPH<19> GLYPH<2> GLYPH<18> 𝑁 O-stars 360 GLYPH<19> GLYPH<2> GLYPH<18> 𝑓 bin 0 GLYPH<149> 47 GLYPH<149> GLYPH<149> GLYPH<149> 0 GLYPH<149> 55 GLYPH<19> GLYPH<2> GLYPH<18> W grid GLYPH<157> W VFTS 0 GLYPH<149> 03 GLYPH<19> GLYPH<149> (9)  \nWe expect a maximum of 2 O-type over-contact binaries or 3 relaxed-contact systems (comprising both the near and over-contact binaries) in 30 Dor. As there are currently 1 confirmed over-contact system and 3 contact/near-contact O-type binaries in the VFTSTMBM sample (see systems marked as 'C/NC' in Table 1), our numbers are in close agreement with the observations.  \nSince 30 Dor hosts about GLYPH<25> 25% of the total number of O stars in the LMC (Crowther, P., private comm.), we can expect up to 𝑁 contact, LMC GLYPH<25> 8 O-type over-contact and up to 12 relaxed-contact (near + over-contact) binaries in the LMC, assuming like 30 Dor that only GLYPH<25> 50% of O-stars are in binaries. If we instead assume that the binary fraction is 100% in the LMC as is the case for O-stars in the Milky Way (e.g., Sana et al. 2014), our predicted numbers go up to 13 over-contact or 21 relaxed O-type contact binaries in the LMC. The total number of observed massive contact binaries in the LMC is 38: 30 in the OGLE+MACHO database, 5 systems in the VFTS-TMBMsampleand3othermiscellaneous ones. Although we do not have the respective contributions of O and B-type stars to the OGLE+MACHOsample,weexpectthat the B-type contact binaries will dominate the sample owing to their higher lifetime. Assuming thus that less than 50% of the known LMC contact binaries are O+O-type systems, our estimated number of these systems is in agreement with the data.  \nThese calculations implicitly assume that there are no strong observational biases and that the star formation rate has been constant for a period of time longer than the lifetime of a 15 M GLYPH<12> star. However, Schneider et al. (2018a,b) suggest that the star formation rate of 30 Dor has been declining over the last 10 Myr. They also expect a shallower IMF for 30 Dor than the Salpeter one, which may be a reason for the abundance of high mass stars GLYPH<21> 30M GLYPH<12> in the VFTS sample. These numbers also depend on the power laws we assume for the birth parameters of our synthetic populations.  \nExtending our calcualtions to the SMC, we first estimate the number of O-type stars in the SMC by comparing its star formation rate of GLYPH<24> 0 GLYPH<149> 05 M GLYPH<12> /yr (Harris & Zaritsky 2009) to that of the LMC GLYPH<24> 0 GLYPH<149> 2 M GLYPH<12> /yr (Hagen et al. 2017). Given thus that the star formation rate of the SMC is 1/4th of the LMC and by assuming f bin = 0 GLYPH<149> 55, we expect a maximum of 2 over-contact or 3 relaxed-contact O-type binaries in the SMC.", '5 DISCUSSION': "From our binary evolution models, we find that nuclear-timescale massive contact binaries evolve towards equal masses. This result agrees with the confirmed over-contact systems, which indeed have equal or nearly equal masses (q GLYPH<25> 1), thus implying that a mass ratio of unity is the equilibrium configuration that contact binaries tend to achieve. We further infer from our synthetic contact-binary populations that these observed equal-mass systems must have had initial mass ratios GLYPH<21> 0 GLYPH<149> 8 and that their fate is to eventually merge on the main sequence. Two of the confirmed over-contact binaries though do not have exactly equal masses- OGLE SMC SC10-108086 which has q = 0 GLYPH<149> 85 GLYPH<6> , 0 GLYPH<149> 06 and an even more exceptional V729 Cyg which  \nhas the lowest mass ratio (q = 0 GLYPH<149> 29) in our examined sample. The latter though is known to harbour a tertiary companion (Kennedy et al. 2010; Rauw et al. 2019).  \nAlthough our predicted mass ratio distribution is heavily skewed towards a value of one, the majority of observed massive contact binaries have unequal masses (q GLYPH<159> 0 GLYPH<149> 9). Upon inspection wefindthat these unequal mass systems are in fact, either confirmed near-contact binaries or those whose degree of contact is unclear. To study their contribution, we constructed a synthetic population that also include binaries in near-contact, in which both stars have R/R L GLYPH<21> 0 GLYPH<149> 9, i.e., they overflow at least 72% of their Roche lobe volumes ((R GLYPH<157> R L ' 3 ). While our definition of near-contact binaries is generous, we find that the distributions constructed with these models can accommodate observed systems with mass ratios as low as q = 0 GLYPH<149> 4. We thus conclude that the observed unequal-mass contact systems are likely to be binaries just nearing contact.  \nThe overall distribution of the 7 LMC spectroscopic contact binaries is better explained by considering the contribution of our near-contact binary models and considering the fact that the VFTS data is skewed towards O-type stars with masses GLYPH<21> 15 M GLYPH<12> . We make very specific predictions about the distribution of their orbital periods, masses and mass ratios. For over-contact O and early B-type systems in both Magellanic Clouds, i.e., where both stars in the binary have R/R RL GLYPH<21> 1, we expect the total mass of the majority of observed contact binaries to be / 45 M GLYPH<12> and to have orbital periods less than 1.5 days, with the most likely period to be between 0.7 and 1 day. With a more relaxed definition of contact, i.e., by also including ner-contact binaries, we predict orbital periods to be up to 2 days. Such low orbital periods are reported for over half of the massive contact systems from the OGLE and MACHOsamples (Rucinski 1999; Pawlak 2016). However, as these samples also include late B-type binaries and the individual orbital parameters of the systems are unavailable, it is difficult to isolate their contributions to the sample.  \nOur predicted distributions indicate a strong correlation between the period and total mass of a contact binary, and also between its period and mass ratio, which we also find in the observed sample of contact binaries- wider binaries typically have higher masses and lower mass ratios (Figures 9 and A2). Despite this interesting agreement, which is further fortified by including near-contact binaries or by only retaining O-type stars in our synthetic populations, we still find that our predicted mass ratio distribution is far too skewed towards q=1 and, the correlation between the period and total mass of the observed systems (for total masses GLYPH<21> 40 M GLYPH<12> , especially in the Galaxy) is not as tight as our predictions. In order to evaluate these findings, let us examine the contribution of the uncertainties in the physics of our binary evolution models and in the observations of massive contact binaries.  \nAs in all massive star models, the main uncertainties relate to internal transport processes and mass loss. However, in contact binaries, some of these processes have a particular emphasis. One important issue may be that while our employed contact treatment does consider the mechanical component of a common envelope of both stars within the Roche approximation well, the thermal component is not dealt with. Despite their common envelope, no energy transport from one star to the other is considered in our models, and our contact scheme allows both stars to possess different surface temperatures. An argument that this may not cause large errors is that the internal structure of models of stars is largely independent of the outer boundary conditions (Kippenhahn & Weigert 1990). Besides, the two components of our binary models differ typically by less than 20% in their effective temperatures, and have a decreas-  \ning trend with time as the models tend towards a mass ratio of one. A temperature difference of this order between both components is confirmed by observations of contact binaries (cf., Tables 1, 2 and, Fig. 8) and more recently by Abdul-Masih et al. (2021) for the overcontact binaries, VFTS 352, V 382 Cyg and SMC-SC10 108086. They speculate however, that a heat exchange mechanism may be responsible for the near equalization of the surface temperatures of these stars.  \nInternal chemical mixing is also particularly uncertain in close binaries. It has been postulated that Eddington-Sweet circulations in tidally-locked binaries can lead to chemically homogeneous evolution (e.g., de Mink et al. 2009; Marchant et al. 2016; Hastings et al. 2020). Hastings et al. (2020) showed that the Eddington-Sweet circulations in tidally-locked binaries can be up to twice faster than in a comparable rotating single star. This effect is not included in our models. However, even with the enhanced Eddington-Sweet circulations, initial total masses above nearly 100 M GLYPH<12> are required for this channel to operate, which is beyond the upper limit of our model grids. From the observed contact systems, VFTS 217 is closest to this regime. However, Mahy et al. (2020b) find that neither of the components of this system is overluminous, although it should be borne in mind that this system is just approaching contact. On the other hand, Abdul-Masih et al. (2021) show that likely both components, but certainly the less massive one, of the Galactic system V382Cyg and SMC system OGLE SMC-SC10108086 are overluminous. This might argue for strong internal mixing, however, the surface abundances of the six stars in these three binaries are not strongly enriched (Abdul-Masih et al. 2019, 2021). While we can not draw any firm conclusion here, internal mixing remains only a weak candidate for explaining the discrepancies between our models and some of the observed binaries.  \nMass loss may also be an important process in some of the binary systems considered here. First, in particular while comparing with the most massive and metal-rich contact binaries, their strong stellar winds must be accounted for. While stellar winds are included in our models, the mass loss rates of Galactic O stars are thought to be about 1.6 times higher than of those in the LMC (Vink et al. 2001; Mokiem et al. 2007), which could be relevant for the most massive Galactic contact systems especially to explain their wide orbits. While we can only speculate how the orbital-parameter distributions of massive Galactic contact systems may behave, a separate grid of MWbinary models is required to ascertain the true nature of these distributions. We note that the widest known LMC (near) contact system, VFTS 217, which is the only clear outlier in the Pcontactqcontact diagram of the Magellanic Clouds sample (Fig. 9), has a total mass above 80 M GLYPH<12> . Among the MW contact binaries, LY Aur, XZ Cep, V348 Car and V729 Cyg have exceptionally wide periods ( GLYPH<21> 4 days) compared to the predicted period domain from our distributions. They have total masses ' 40 M GLYPH<12> , except for XZ Cep which although reported as a contact system, is in fact a semidetached binary.  \nA clue for understanding these wide and massive metal-rich systems can be found by examining the HRD positions of their components in Fig. 8. Our models predict that short-period contact systems should be found close to the ZAMS on the HRD, which is in agreement with the majority of the observed systems. This situation is only different for the Galactic systems with periods GLYPH<21> 4 days, which are much further evolved along the MS. Obviously, here both components keep overfilling their Roche lobes despite the large system dimensions. One reason they may do that could be their proximity to the Eddington limit, which inflates the envelopes of these stars (Sanyal et al. 2015). Since envelope inflation is metal-  \nlicity dependent (Sanyal et al. 2017, and see Fig. 8), this effect is most strongly expected in the MW, while in our lower metallicity models, it does not play a significant role.  \nMass-transfer related mass loss from binary systems may also be important here. Notably, our model calculations are conservative, in that all the mass transferred from the mass donor is assumed to end up on the mass gainer. It is well known that many mass transferring binaries can spill a considerable fraction of the transferred matter out of the binary system (de Mink et al. 2007; Langer 2012). In fact, one of the Galactic contact binaries, SV Cen does have an accretion disk detected around it (Linnell & Scheick 1991; Shematovich et al. 2017). While the dependence of this binary mass loss on the system parameters is not well known, and the slow, nuclear timescale mass exchange in contact systems may be close to conservative, the thermal timescale mass transfer which starts the contact binary evolution (cf., Section 3.1) may have well been non-conservative, especially for the systems with the rather large orbital periods (Pcontact ' 1 GLYPH<149> 2 days) as is the case for some of the most massive observed contact binaries. A possible mode of nonconservative mass transfer in our models may occur via isotropic re-emission from the mass gainer after it accretes mass (Soberman et al. 1997), which may increase the specific angular momentum of the binary and widen its orbit compared to what we will be obtained from a conservative assumption of mass transfer (Sen et al. in prep; Hastings et al. in prep).  \nIn our study we have not addressed the actual formation of massive binary systems with initial orbital periods within the range considered here, i.e., between 0.6 and 2.0 days. Though beyond the scope of this study, this is a relevant problem as pre-main-sequence (PMS) stars have larger radii relative to ZAMS stars, facilitating interaction during the PMS-phase, and ZAMS component masses are correlated. Studies of lower mass stars consider Kozai-Lidov oscillations and dynamical instability in triple systems as channels toward the formation of close binary systems (Sana et al. 2017; Moe & Kratter 2018). Indeed some of the MW contact binaries do have tertiary companions or binary systems detected, which may have also played a role in their current orbital period, inclination and eccentricity. For massive stars, interaction with primordial gas regulated by an external magnetic field (Lund & Bonnell 2018), along with disc fragmentation and the subsequent migration of protostars forming in these fragments (Meyer et al. 2018) have been proposed as well. That orbital hardening in general is an essential ingredient of close binary formation is supported by recent findings by Ramírez-Tannus et al. (2021), who find evidence for migration by comparing binary properties in young OB associations. Finally, we stress that our assumptions for the initial distribution of periods and mass ratios are subject to uncertainties. Observations of systems with masses and periods of interest are scarce. They are typically derived from OB associations and clusters that are estimated to be 0-4 Myrs old. Some systems may have already undergone some evolution and mass exchange by the time we can observe them (such as binaries that evolve like our example system 1). Our predictions for the observable P-q distribution are therefore likely somewhat sensitive to the precise choices made.  \nOverall, we can say that our contact binary evolution models appear to reproduce the observed properties of contact binaries in the lower total mass range (below about 45 M GLYPH<12> ) reasonably well, particularly for the metallicities of the Magellanic Clouds. Our models however, appear to overestimate the fraction of contact binaries with mass ratios close to one and fall short of explaining the scatter in the period-total mass distribution of the more massive contact binaries, especially in the MW.", '6 CONCLUSIONS': 'We have performed the first dedicated investigation of the evolution of massive contact binaries using our current theoretical understanding of this phase and made predictions for their observed population. We have computed two large grids of detailed binary evolution models for the LMC and SMC metallicities in MESA , by assuming mass transfer to be conservative and by using a contact prescription that considers the mechanical consequences of a common envelope within the Roche model, but neglects any energy exchange between the components of the contact binary.  \nFrom our models, we find that binaries that are in contact over nuclear timescales will inevitably merge. Hence we predict that all the observed contact binaries will eventually merge on the main sequence. From our models we find that a stellar merger is the outcome for most binaries in the Magellanic Clouds with initial periods less than 2 days and a definitive fate for all binaries with periods less than 1.5 days. This contact duration can last up to 10.7 Myrs in the lowest mass binaries in our grid, before they ultimately merge.  \nOur nuclear-timescale contact binary models evolve towards equal component masses, corresponding to the lowest energy state of binary systems (Adams et al. 2020). The longer a system spends time in contact, the likelier it attains a mass ratio of 1 on the main sequence. Nearly all the observed systems in deep contact and some in near contact have mass ratios between 0.9 and 1, and hence agree well with our result. From our models, we infer that these contact binaries would have had initial mass ratios GLYPH<21> 0 GLYPH<149> 8. We also expect massive binaries to spend less time in contact as their initial mass or initial period increases, and when their initial mass ratio deviates further from 1. None of our contact binary models undergo chemically homogeneous evolution in both grids.  \nWe expect nearly 4.5% of the O-star binaries in the LMC and SMC to be over-contact binaries (10 systems in total), assuming a binary fraction of GLYPH<25> 50% and constant star formation. With a binary fraction of 100%, the number of over-contact binaries will be 1.6 times as many and even more if we include binaries in near-contact (R/R RL GLYPH<21> 0 GLYPH<149> 9 for both stars).  \nWe find it essential to consider these near-contact binaries in our study due to the observational uncertainties in the degree of contact and in the configuration of reported contact binaries. We also find that their inclusion while computing the distribution of orbital period and mass ratio, reproduces the observed distribution more closely. In a population of O & early B-type binaries in the Magellanic Clouds, we expect that the contact binaries will have total mass / 45 M GLYPH<12> , periods GLYPH<159> 1 GLYPH<149> 5 days and mass ratios GLYPH<21> 0 GLYPH<149> 9. By including near-contact binaries, our distributions predict a wider domain of probable orbital periods of up to 2 days and of mass ratios as low as 0.4. From our modelled distributions, we predict that all reported contact binaries which have mass ratios / 0 GLYPH<149> 80 and periods greater than 1 day, are likely to be near-contact systems.  \nWe find little difference in the stellar and orbital-parameter distributions between our LMC and SMC contact binary populations, implying that metallicity does not play an important role in their orbital evolution at the values we have considered. Significant differences only appear in the orbital period distribution if we place mass cuts on our synthetic populations, by excluding systems with total masses less than 30 M GLYPH<12> .  \nOur binary models identify a distinct subspace in the orbital period-mass ratio diagram in which most contact systems are expected (Fig. 9), which we find indeed populated by most of the observed Magellanic Cloud contact binaries with total masses below 80 M GLYPH<12> . Like the general trend seen in the observed data, our  \nmodels also predict that wider binaries tend to be more massive and have small mass ratios. Overall, this argues for our treatment of the contact phase to capture the essential mechanics of the orbital evolution of the binaries during this phase, especially for systems with total masses GLYPH<21> 45 M GLYPH<12> . For higher mass systems we cannot reproduce the individual values of their mass ratios or orbital periods, especially when we consider the Galactic systems. This is mainly due to the fact that lower-mass systems are weighted higher in a distribution as they are favoured more by the initial mass function and their overall longer lifetimes. While imposing a mass cut on the distribution or including near-contact binaries does somewhat alleviate the problem, other physics phenomena may also be at play in the evolution of these massive systems such as stellar winds, non-conservative mass transfer and envelope inflation, which may be especially relevant for the Milky Way binaries.  \nTo test these ideas, we need to compute dedicated model grids for the Milky Way to compare with the Galactic massive contact binaries. In Paper II, we intend to explore these avenues and compare the stellar parameters of our models with the currently observed massive contact binaries, such as their luminosities, effective temperatures, rotational velocities and surface abundances.  \nWe conclude thus that the contact phase is an inevitable, nuclear-timescale evolutionary phase for most massive binary stars born with periods less than 2 days. We also identify it an observable phase preceding an imminent stellar merger, lending us the opportunity to probe the conditions that lead to one of the least understood evolutionary phases in stellar physics. While stellar mergers themselves produce exotic and peculiar stellar phenomena and supernovae, the contact phase in extremely metal-poor and massive binaries has also been flagged as a channel to produce the progenitors of merging black holes (Marchant et al. 2016). Therefore we hope that, ultimately, our results may not only lead to a better understanding of massive contact binaries in general, but also lead to better predictions of progenitors of stellar mergers and of massive compact object mergers, not only of those in our cosmic neighbourhood but also those which live so far out in the Universe that direct observations may remain difficult for the foreseeable time.', 'ACKNOWLEDGEMENTS': "This project was funded in part by the European Union's Horizon 2020 research and innovation program from the European Research Council (ERC, Grant agreement No. 715063, 772225: MULTIPLE), and by the Netherlands Organization for Scientific Research (NWO) as part of the Vidi research program BinWaves with project number 639.042.728 and the National Science Foundation under Grant No. (NSF grant number 2009131). D.Sz. is supported by the Alexander von Humboldt Foundation. PM acknowledges support from the FWO junior postdoctoral fellowship No. 12ZY520N. SdM is grateful to all members of the BinCosmos group for enlightening discussions and input, in particular to Lieke van Son and Katie Sharpe.", 'A1 Probability distributions of contact systems according to their initial parameters': 'In this section we elaborate how we compute the distribution of contact systems as a function of their initial parameters: P i GLYPH<150> q i and M T,i . To compute how contact systems are distributed as a function  \nof P i , we first compute the weighted contact duration of each system and then the contribution of that system towards a given P i :  \nHcontact \' P i \' = 𝛿 Pi 𝑤 s 𝜏 contact,s " n s=1 𝑤 s 𝜏 contact,s (A1)  \nwhere, P i 𝜖 » 0 GLYPH<149> 6 GLYPH<150> 0 GLYPH<149> 7 GLYPH<150> 0 GLYPH<149> 8 GLYPH<150> GLYPH<149> GLYPH<149> GLYPH<149> GLYPH<150> 2 GLYPH<149> 0 … days and 𝛿 Pi takes the value of either 0 or 1. For e.g., for systems with an initial period of 1 day, 𝛿 P1.0 ; for all other systems it takes the value of 0. 𝜏 contact,s is the net time spent in contact by the system and ws is the birth weight of the system as calculated in Eqs. 4.  \nWe compute similar functions for each q i and M T,i in our grid. These are later used to make 2D probability distributions of P i vs. q i and P i vs. M T,i (as will be seen in Fig. 9).', 'A2 Probability distributions of contact systems according to observed parameters': 'In this section we show how we compute the distribution of contact systems according to their present-day parameters.  \nWe divide the whole range of orbital periods for contact systems during their evolution in bins of 0.1 days. We take each system from the grid and determine how much time it spends in contact in each bin and weigh this time by its initial weight. To construct the distribution of the observed orbital period, we sum this value in each period bin and normalize its contribution to the total weighted contact duration of the entire grid (Hcontact \' Δ P \' ).  \nHcontact \' Δ P \' = 𝛿 Δ P " n s=1 𝑤 s dts, Δ P " n s=1 𝑤 s dts, Δ P (A2)  \nwhere, ws is the birth weight of the system, Δ P is the orbital period bin of systems while they are in contact such that Δ P 𝜖 » 0 GLYPH<149> 6 GLYPH<0> 0 GLYPH<149> 7 GLYPH<150> 0 GLYPH<149> 7 GLYPH<0> 0 GLYPH<149> 8 GLYPH<150> GLYPH<149> GLYPH<149> GLYPH<149> GLYPH<150> 3 GLYPH<149> 0 GLYPH<0> 3 GLYPH<149> 1 … days, dt s, Δ P is the time spent by the system while its in contact in each Δ P bin. 𝛿 Δ P is either 0 or 1; if a system spends time in contact in the bin Δ P = 0 GLYPH<149> 6 GLYPH<0> 0 GLYPH<149> 7 days, then 𝛿 0 GLYPH<149> 6 GLYPH<0> 0 GLYPH<149> 7 = 1 else it is 0.  \nFrom the distribution thus constructed, we can compute the probability of finding a contact system with a particular current orbital period and hence their most likely observed period. We similarly compute H contact, Δ q and H contact, Δ MT and the 2D probabilities of finding contact systems for combinations of mass ratio and total mass.', 'DATA AVAILABILITY': 'The template folders to run our single star and binary models used in this work are available at : https://github.com/athira11/ massive\\_contact\\_binaries.git . If they are used, please cite this work as the source of the data. If particular models are required, please send an email to the first author.  \nFigure A1. The contact duration ( 𝜏 contact) of each system of the entire LMC grid in a Pi-MT GLYPH<150> i for all q i . Background colors and symbols are as explained in Fig. 6.  \n<!-- image -->  \nglyph[circledot]  \nglyph[circledot]  \nglyph[circledot]', 'Models of massive contact binaries 21': 'Figure A2. The 2D distribution of Pcontact vs. qcontact for contact and over-contact binaries, with mass cuts on initial total mass: MT,i GLYPH<21> 30 M GLYPH<12> (left) and MT,i GLYPH<21> 50 M GLYPH<12> (right). Colors, histograms and markers are as explained in Fig. 9.  \n<!-- image -->', 'REFERENCES': 'This paper has been typeset from a T E X/L A T E X file prepared by the author.'}