Thesis work related of Pablo Marchant (both personal and supervised)

ARENBERG DOCTORAL SCHOOL Faculty of Science Detailed structure and evolution modeling of the tightest massive binary stars Matthias FABRY Examination committee: prof. dr. J. Magdaleni´c Zhukov, chair prof. dr. H. Sana, supervisor dr. P. Marchant, co-supervisor prof. dr. J. O. Sundqvist prof. dr. C. Aerts prof. dr. N. Langer (Universität Bonn, Germany) dr. E. Laplace (Heidelberg Institute for Theoretical Studies, Germany) prof. dr. O. R. Pols (Radboud Universiteit Nijmegen, the Netherlands) Dissertation presented in partial fulfillment of the requirements for the degree of Doctor of Science (PhD): Astronomy and Astrophysics May 2024 Cover Art: Inez Fabry Acknowledgments: M.F. received funding from the Flemish research foundation (FWO, Fonds voor Wetenshappelijk Onderzoek) under PhD fellowship No. 11H2421N. The research leading to these results has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement numbers 772225: MULTIPLES). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa. int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/ web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. Based on observations made with the Mercator Telescope, operated on the island of La Palma by the Flemish Community, at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. Figure 1.12 reproduced under the fair use policy of Art. XI.191/1 of Belgian Law. No copyright infringement is intended. ©2024 KU Leuven – Faculty of Science Uitgegeven in eigen beheer, Matthias Fabry, Celestijnenlaan 200D box 2401, B-3001 Leuven (Belgium) Alle rechten voorbehouden. Niets uit deze uitgave mag worden vermenigvuldigd en/of openbaar gemaakt worden door middel van druk, fotokopie, microfilm, elektronisch of op welke andere wijze ook zonder voorafgaande schriftelijke toestemming van de uitgever. All rights reserved. No part of the publication may be reproduced in any form by print, photoprint, microfilm, electronic or any other means without written permission from the publisher. Acknowledgments What a journey. In the first half of my PhD, the pandemic required all interactions with students and colleagues to move to Zoom (2/10; would not recommend). Despite the unbelievable situation, once it receded, I had the pleasure to meet many cool people abroad, both at conferences and during visits. Special thanks go to Silvia, Tomer and the API members for hosting me in Amsterdam, and to Norbert and the Bonn stellar group for letting me stay for four months! Also big thanks to the MESA devs, TAs and students at the Santa Barbara summer school for your inexhaustible enthusiasm for stars. Back home, I had a most wonderful time at the Institute of Astronomy here in Leuven. It really is a vibrant place, with lots of sports, receptions, and an especially tight-knit group of PhD students and post-docs. Special thanks go to: The volleyball group (Robin, Karan, Nico, Olivier, Luka, and many, many more...) for the many Friday nights on the sand (and at Spuye or Metafoor afterwards...). The 2023 Uitje team (too many to list here, THANKS everyone!) for organizing three days in the Ardennes full of activities. The UK/Ireland folks (Abi, Gareth, Calum, Emily) who introduced me to Six Nations Rugby at the Irish Pub. The SkIvS adventurers (Nich, Luc, JP, Silke, ...) whom I was able to join twice, taking in the sights, showing off the drip and shredding the gnar on the French alps. My fellow San Francisco travelers (Mathias, Joey, Vincent, Annachiara) for sitting through my very first live baseball game with me, walking through the Californian Redwoods, driving down highway 1, and otherwise making my first trip to the US one to never forget. Clio, Anita, Monique and the System team who all play invaluable support roles in the institute. Jesus and Saskia for the support on my three trips to La Palma. Observing on the mountain is truly an unforgettable experience! i ii ACKNOWLEDGMENTS And finally, massive thanks to the massive-star group. I consider myself blessed to have been able to meet every single one of you. Thanks Abi, Calum, Gareth, Laurent, Annachiara, Kunal, Soetkin, Tinne, Emma, Alan, Frank, Jaime, Julia, Michael, Cyprien, Philippe, Maddalena, Tomer, Reinhold, Dominic and Karan for the many fruitful discussions, interesting science and the fun on the team Uitjes! I deliberately left Hugues and Pablo out of this list. That’s because, of course, you two deserve separate paragraphs. Merci mille fois, Hugues, d’être un super mentor. I consider starting a PhD with you one of the best decisions I have ever made in my life (9/10; would recommend). Even before completing my work on your project, hell, even before finishing the first three months, you immediately let me write my own project, which resulted in the bulk of this thesis. You kick-started my academic career (which I hope to continue for many years to come) and supported me in my personal and professional development. I especially like your no-nonsense attitude, and making it extremely clear what to do and also (maybe more importantly) what not to do. If ever I need an opinion on anything even remotely related to astrophysics, I know who to turn to. Muchas gracias, Pablo, por mostrarme las estrellas. While I learned what stars are in my courses, it was you who taught me the fine details and showed me the inner workings of MESA , not to mention the myriad effects of binary processes. You have provided the majority of the scientific motivation for this thesis, which I will never take for granted. Alongside Hugues, you have supported me throughout the past years, and thanks to you, I have become the scientist I am today. You are a continuous inspiration to me and I will always look to you and your work as the gold standard for research in theoretical stellar astrophysics. I hope we can work together for many more years. Thanks also for being a cool guy and a good friend, and 3D printing some of the designs I made over the years. I’d be remiss not thanking all my other friends. Thanks to my pals at Golf Park Tervuren (and formerly DMGC), for all the fun rounds on the many Sunday mornings. I also have good memories from our Interclub competition campaigns! Y’all made it a really, really fun club I feel at home at. The little nine-hole course I learned to play on will always have a special place in my heart. I’ll definitely visit in the future! Thanks to the people of Guikui (who came up with this name anyway?), for dealing with me during our bachelor degree. You are a bit crazy sometimes, but that’s OK, and exactly why I like you. I was maybe not the most outgoing in this group during my PhD, but I still want to thank you since you helped me getting through some of the toughest times of my life. Big thanks to my three partners of the Lucky Basterds, Bart, Jordi and Lennart. By my estimation we played board games together for over a thousand hours, adventuring in Gloomhaven, Frosthaven, the 7th continent, Aeon Trespass, Scythe, ACKNOWLEDGMENTS iii and elsewhere. It takes a special kind of geek to be willing to do this, and I’m happy to have found them in you guys. I think it’s all worth it though, right? Remember tick tock, the (glorious) deus ex machina d10 crit roll in Aeon, me trying desperately to make expeditions work in Odin, and, of course, the terrible semi-cooperative shenanigans in Gloom/Frosthaven. On to 10000? For some more random acknowledgments: Thanks to deadmau5 for slapping some chords together that I enjoy, 1 Frank Herbert for writing Dune, the Miller brothers and Cyan Worlds for creating Myst and the D’ni universe, James SA Corey for writing the Expanse, Mike Pondsmith for creating Cyberpunk upon which both CDPR’s 2077 and Netrunner are based (not random: huge thanks to the Belgian Netrunner community, but especially Ruben for playing many games with me Tuesday nights in Demo-Spel!), Asobo studios for creating the best flight simulator to date, with which I can feel just a bit closer to my late grandfather, and Jon van Caneghem for developing Heroes of Might and Magic III, a simply awesome and timeless game. These people will probably never read this, but still, thanks, you kept me busy during most of my off-time in the past years, and that’s equally important I think. Ik heb de belangrijksten als laatst gehouden. Dank u meter, oma en opa om interesse te blijven tonen in mijn sterrenkundige avonturen. Wat ik bestudeer ligt misschien wat ver van jullie bed, maar toch laten jullie mij maar doen, met mijn hoofd in de ruimte. Dank u peter om mijn science-nonkel te zijn. Klinkt als iets uit een Amerikaanse sitcom, maar ik heb het geluk een échte te hebben! Dankzij jou weet ik wat mogelijk is als academicus, en ook bij wie ik terecht kan voor persoonlijk advies. Ten laatste, dank u mama, papa en zus voor de ongeziene steun. Dank u om ten minste alsof te doen geïnteresseerd te zijn in wat ik doe, en te blijven luisteren naar mijn gepalaver over sterren, pindanoten en zwarte gaten. De voorbije jaren waren niet altijd de gemakkelijkste, maar desondanks heb ik dit boekje kunnen schrijven, en dit was niet mogelijk zonder jullie hulp gedurende mijn hele leven. Ik ben hier extreem fier op, en ik hoop jullie ook. Matthias 1like “some chords” – deadmau5 Abstract Massive stars are responsible for producing heavy elements in the universe, and scatter them into their cosmic neighborhoods as they explode as supernovae. These processes pave the way for the formation of a new generation of stars that have the possibility to form planets out of the formed dust, potentially harboring life. With masses over eight times that of our Sun, massive stars lead bright and fast lives, taking only millions of years from birth to death, compared to billions for solar-type stars. Furthermore, much more so than the low-mass stars, massive stars come predominantly in pairs, called binary systems. Throughout their evolution, stars in binary systems interact, transferring mass between each other, which has major implications for their ultimate fate. This is one of the reasons that makes studying massive binary stars so interesting. For tight binary systems in particular, a whole list of extra effects come into play. Some binaries orbit so close that they deform through tidal interaction (the same process that sloshes the oceans on Earth), and in the extreme case, the stars come into contact, forming “peanut stars.” Massive contact binaries do not necessarily immediately merge into one, as they are rotating around each other very quickly, and they can be stable for millions of years. What is more is that, by observing samples of contact binaries, they seem to be stable while having unequal masses. The mass-ratio distribution is fairly uniform between mass ratios of around one half and unity. Unfortunately, no binary models that simulate the life of these stars have been able to reproduce this constraint, and predict that contact systems should overwhelmingly appear as equal-mass binaries. When modeling the life of stars, today, we use computers to solve the set of differential equations that govern their evolution. We include several effects that modify the structure of the stars, such as rotation, internal mixing, stellar winds, as well as the binary-star process of mass transfer. Up until today, tidal deformation and the process of energy transfer in contact binaries have not yet been accounted for in detailed models of binary evolution. Instead, the state-of-the-art models are spherically symmetric, or at most deformed cylindrically due to rotation. v Vulgariserende samenvatting Sterren zijn grote ballen gas en plasma dat onstaan door de samentrekking van stofen gaswolken in het heelal. Hun massa’s kunnen variëren van enkele honderdsten tot ettelijke tientallen maal de massa van onze Zon, dewelke op zijn beurt miljoenen keren zwaarder is dan de hele Aarde. Als een ster meer dan ongeveer acht keer zwaarder is dan de Zon, spreken we van een zware, of massieve ster. Massieve sterren zijn werkelijke kosmische machines, omdat zij zware elementen, zoals ijzer, in hun kernen maken en deze dan later doorheen hun melkwegstelsel slingeren wanneer zij exploderen in een supernova. Laagmassieve sterren zoals de Zon hebben deze eigenschappen niet. Hierdoor is het bestuderen van massieve sterren heel belangrijk, omdat zij verantwoordelijk zijn voor het creëren van elementen dat leven mogelijk maakt in het heelal. IJzer is bijvoorbeeld een essentieel ingrediënt in het menselijk bloed. Aan de hand van geavanceerde (computer)modellen kunnen astrofysici uitspraken doen over het leven van een ster. Zo kunnen we bijvoorbeeld de temperatuur of de levensduur van een ster berekenen, maar we kunnen ook hun binnenste structuur bestuderen. In de laatste honderd jaar zijn we zo te weten gekomen dat sterren hun energie uit nucleaire fusiereacties halen. De energie die in het centrum van de ster geproduceerd wordt, wordt dan naar de buitenste lagen getransporteerd oftewel door straling, oftewel door “te bubbelen,” net zoals water kookt in een pot op een gasvuur. Al deze processen hebben als einddoel de structuur van de ster stabiel te houden. Moest er immers geen energie geproduceerd worden, zouden sterren als een kaartenhuis in elkaar zakken ten gevolge van hun immens gewicht! Massieve sterren komen heel vaak voor met een partner, véél meer dan hun laagmassieve soortgenoten. Door hun wederzijdse aantrekkingskracht draaien zij rond elkaar, net als een koppel op de dansvloer. Zulke formatie noemen we een binaire ster, of ook een binair systeem. Soms vormen sterren zó dicht bij elkaar dat doorheen hun levens hun buitenste lagen in contact komen, een situatie waarvan u op de cover een artistieke illustratie kan beschouwen. Hun vorm lijkt wel op een pindanoot, en de naam “pindanoot-ster” blijft ook bij wetenschappers plakken. Dit is een zeer extreme fase in dewelke binaire sterren kunnen voorkomen, maar xiii xiv VULGARISERENDE SAMENVATTING toch gebeurt dit vrij frequent. Gedetailleerde modellen leren ons dat ongeveer 40% van alle binaire paren van massieve sterren een contactfase meemaken. Om goede voorspellingen te doen over het leven en de eigenschappen van binaire sterren, moet de contactfase dus nauw onder de loep genomen worden. In een computermodel van een (binaire) ster is het belangrijk dat alle relevante fysische processen mee worden berekend. In het bijzonder zijn er in een contactsysteem namelijk extra effecten die toegevoegd moeten worden, die tot nog toe enkel benaderd of zelfs helemaal verwaarloosd werden. Bijvoorbeeld worden sterren vaak voorgesteld als perfect ronde bollen, maar het moet duidelijk zijn door de illustratie op de cover dat dit niet klopt voor de sterren in een contactsysteem. De vervorming van de sterren is het gevolg van getijdenwerking, dezelfde kracht dat de zeeën op Aarde op en neer doet gaan. Deze thesis bouwt een methode om zulke vervorming voor te stellen in onze computermodellen. Een tweede effect in contactsystemen dat verwaarloosd werd is energietransfer. Omdat de buitenste lagen in zulke systemen overlappen, kan er energie vloeien van de ene ster naar de andere. Hoe de warmtestroom in een contactsysteem er precies uit ziet, is nog onduidelijk. De verwachting is echter wel dat energie vloeit van de meer heldere ster naar de minder heldere ster, wat vergeleken kan worden met hoe warmte vloeit van binnenin een hete radiator naar een koele kamer erbuiten. Dit verwachte resultaat modelleren we in deze thesis, om zo het effect van energietransfer op het leven van een contactsysteem te bestuderen. Uit de nieuwe computermodellen van contactsystemen, waarin we de effecten van getijdenvervorming en energietransfer in rekening brachten, moeten we besluiten dat zij nog steeds de eigenschappen van waargenomen contactsystemen niet helemaal kunnen verklaren. Ookal verbeteren zij enkele aspecten, moeten er dus nog verschillende andere effecten of processen zijn die we in de modellen over het hoofd zien. Het laatste woord over contactsystemen is dan ook nog niet geschreven. Desondanks zijn binaire systemen in contact unieke laboratoria waarmee astrofysici de theorieën van sterevolutie kunnen testen, om zo steeds betere en betere stermodellen te bouwen. AGB Asymptotic giant branch (star) BC Boundary condition BH Black hole CDF Cumulative distribution function CE Common envelope CHE Chemically homogeneous evolution ET Energy transfer GRB Gamma-ray burst GW Gravitational wave HMXB High-mass X-ray binary HRD Hertzsprung-Russell diagram MB Magnetic braking MC Monte Carlo MCMC Markov chain Monte Carlo (method) MESA Modules for Experiments in Stellar Astrophysics (code) MT Mass transfer NIR Near infrared PMS Pre-main sequence RGB Red giant branch (star) RL Roche lobe RLOF Roche-lobe overflow RV Radial velocity S/N Signal-to-noise ratio SB Spectroscopic binary SN Supernova WR Wolf-Rayet (star) ZAMS Zero-age main sequence ZKL von Zeipel-Kozai-Lidov (cycles) List of acronyms xv Contents Abstract v Beknopte samenvatting ix Vulgariserende samenvatting xiii xv List of acronyms xv Contents xvii 1 Introduction 1 1.1 Single stars ................................. 4 1.1.1 The equations ........................... 10 1.1.2 Star formation ........................... 19 1.1.3 Main-sequence evolution .................... 19 1.1.4 Post-main-sequence evolution ................. 22 1.1.5 Open questions in massive-star structure and evolution . . . 25 1.2 Binary stars ................................ 30 1.2.1 Orbital characterization and dynamical masses ........ 34 1.2.2 Stellar multiplicity ........................ 37 1.2.3 The Roche potential and mass transfer ............. 39 1.3 Contact binaries .............................. 44 1.3.1 Observed systems ........................ 45 1.3.2 Modeling efforts ......................... 46 1.3.3 Evolutionary pathways – Mergers ............... 48 1.4 Gravitational wave astronomy ..................... 50 1.4.1 Gravitational-wave progenitors ................. 53 1.5 This thesis ................................. 57 2 The dynamical mass of 9 Sagittarii 61 xvii xviii CONTENTS 2.1 Introduction ................................ 62 2.2 Observations ................................ 64 2.2.1 Optical spectra .......................... 64 2.2.2 Near-infrared Interferometry .................. 65 2.3 Orbital analysis .............................. 67 2.3.1 Astrometric orbit ......................... 67 2.3.2 Spectral disentangling ...................... 68 2.3.3 Distance .............................. 72 2.4 Atmosphere modeling .......................... 74 2.4.1 Setup ................................ 74 2.4.2 Results and discussion ...................... 75 2.5 Evolutionary modeling .......................... 76 2.6 Conclusions ................................ 78 3 Modeling tidal distortion in 1D stellar structure 81 3.1 Introduction ................................ 82 3.2 Methods .................................. 84 3.2.1 Modifications to the spherical stellar structure equations . . 84 3.2.2 Single star rotation ........................ 86 3.2.3 Synchronized Roche binaries .................. 87 3.2.4 Contact shells ........................... 88 3.2.5 Numerical calculations ...................... 89 3.2.6 Integration results ........................ 90 3.2.7 Comparison to literature ..................... 91 3.3 The 𝛺𝛤 limit ................................ 92 3.4 Atmospheric boundary conditions ................... 96 3.5 Stellar model comparison ........................ 98 3.5.1 Physical ingredients ....................... 99 3.5.2 Main sequence evolution of a single rotating star . . . . . . . 100 3.5.3 Evolution of a detached B-type star ..............102 3.5.4 Evolution of a twin contact binary ...............102 3.6 Discussion and conclusions .......................104 4 Modeling energy transfer in contact binaries 107 4.1 Introduction ................................108 4.2 Theory of energy transfer ........................109 4.2.1 Simple considerations ......................109 4.2.2 Models of energy transfer ....................111 4.2.3 Energy transfer in the Roche geometry ............112 4.3 Methods ..................................112 4.3.1 Physical assumptions in MESA ..................112 4.3.2 Shellularity and Roche lobe geometry .............113 4.3.3 Mass transfer ...........................115 4.3.4 Energy transfer ..........................115 CONTENTS xix 4.4 Stellar models ...............................117 4.4.1 Detailed example of ET evolution ...............117 4.4.2 Mass versus luminosity ratios ..................122 4.5 Conclusions ................................124 5 Energy transfer in a population of massive contact binaries 125 5.1 Introduction ................................125 5.2 Methodology ...............................126 5.2.1 Physical assumptions in MESA ..................126 5.2.2 Binary initialization .......................127 5.2.3 Outcomes and termination ...................128 5.2.4 Population synthesis computations ..............128 5.3 Results and discussion ..........................130 5.3.1 Roche-lobe overflow at the zero-age main sequence . . . . . 130 5.3.2 Interaction on the main sequence ................137 5.3.3 Luminosity ratios .........................141 5.3.4 Rejuvenation ...........................142 5.3.5 Uncertain population distributions ...............145 5.4 Conclusions ................................148 6 Conclusions and outlook 151 6.1 Future work and perspectives ......................153 A Additional figures and tables to Chapter 2 161 B spinOS 169 C Distortion in the single rotating star potential 173 C.1 Polynomial fits in the rotating shell potential .............173 C.2 Matching the tidal and single rotating star distortion models . . . . 175 D Additional considerations of energy transfer 177 D.1 Energy transfer location .........................177 D.2 Convergence study ............................178 D.2.1 Explicit evaluation ........................178 D.2.2 Implicit evaluation ........................178 E Additional figures and tables to Chapter 5 181 Bibliography 185 List of Publications 209 Chapter 1 Introduction “Space, the final frontier...” — Captain James T. Kirk When you look up on a clear night in or around Leuven, you can see a thousand or so objects on the sky. If above the horizon, the moon will be most prominent, and it presents itself as a crescent or a bright disk depending on its phase. The next brightest objects are the planets, such as Jupiter, Venus, Mars and Saturn. Even with a modest telescope or good binoculars, you can see Saturn’s rings, or the big red spot of Jupiter. Mercury can only be seen when it is reflecting a good amount of sunlight back to the Earth, while for Uranus you will probably need to squint your eyes in order to see it. Neptune is invisible to the naked eye in any but the very darkest skies, for which the Belgian ones anno 2024 definitely do not qualify. Next in line are objects outside our own solar system, like stars, galaxies and nebulae. In contrast to the planets, these objects are up to thousands of light years away in the case of stars or nebulae in our own Galaxy to millions or even billions of light years for other galaxies. The brightest star in the sky is Sirius, the ‘Dog star’, and can be seen toward the south in wintertime. It is so bright it bears significant importance in many societies and cultures throughout history. It was used by the Egyptians to predict the flooding of the Nile (Clagett,1989), while many Polynesian and other seafaring nations used it (and other bright stars) for navigation purposes (Bateson,1959). In modern day Leuven some of the bright stars you will be able to spot easily include Polaris, the pole star, Betelgeuse and Rigel in the constellation of Orion, Castor and Pollux in Gemini (twins), and others like Vega, Deneb and Altair. When it comes to seeing galaxies, nebulae and star 1 8 INTRODUCTION This number is around ten times bigger than the estimated age of the universe! If the Sun has an age of several billion years, only a few percent of the mass of the Sun needs to be actually available as fusible hydrogen. The stars thus have plenty of nuclear fuel in their tanks to explain their brightness and age. This also means that almost all energy sources on Earth are nuclear in origin. 4 With this final realization, it was Eddington who completed the set of equations that govern the life of a star, the general form of which has remained unchanged since. Further advances were made thanks to new quantum-mechanical insights. More accurate theories on the opacity of stellar material (i.e., the degree with which photons interact with atoms) were developed (Kramers,1923;Rosseland,1924). Additionally, the precise rates at which nuclear reactions take place and the amount of energy released as a result of them, depending on the temperature and density, could be computed from new nuclear theories (Atkinson & Houtermans, 1929;Atkinson,1936;Bethe & Critchfield,1938;Bethe,1939). As a consequence of the exclusion principle put forward by Wolfgang Pauli (Pauli,1925), which states that particles with half-integer spin (such as electrons) cannot occupy the same quantum state, electron degeneracy pressure was discovered and applied to astrophysics by Fowler (1926) in the form of an equation of state independent of temperature. Electron degeneracy would be instrumental in understanding the late stages of stellar evolution. Most notably, Subrahmanyan Chandrasekhar computed, using the relativistic version of electron degeneracy, that there is a maximal mass to an object having this equation of state (Chandrasekhar,1931a,b). Despite a rather public dispute between the well established Eddington, who rejected the theory of the up-and-coming student Chandrasekhar (Mena & Peres, 2018), this idea turned out to be correct, and the limiting mass of a white dwarf (see Sect. 1.1.4) is now known as the Chandrasekhar limit. Until the invention of the multipurpose, digital computer, all calculations involving stellar structure and evolution had to be done by hand, which necessitated that the equations were solvable analytically. This often required making multiple simplifications and approximations of the physical conditions, for example by stating that the gas of the whole star could be described by a single polytrope. A polytrope is the solution of the Lane-Emden equation (Lane,1870; Ritter,1898;Emden,1907), which was the then-prevailing equation describing stellar structure. It assumed a single mode of energy transport inside the star, and could not account for energy generation in the core. Furthermore, Eddington 4 Exceptions to this are tidal power, which originates from the gravitational energy of the Earth, Moon and Sun, and geothermal energy, coming from the heat of Earth’s interior. The latter, in turn, is split in two sources: gravitational from its formation and settling of the core, and nuclear, this time from the radioactive decay of naturally occurring radioisotopes within the Earth. The fraction that the nuclear component contributes to geothermal energy is highly debated, and ranges from 50% (Gando et al.,2011) to 80% (Turcotte & Schubert,2014), depending on the model. SINGLE STARS 9 had to initially consider fully homogeneous stars (with a constant gas-to-photon pressure, a source of debate) to keep the equations manageable, meaning he could only compute models for fully-mixed, uniform stars. The homogeneous stellar models of Eddington were not able to simultaneously explain the two populations in the HRD. You have the giant stars who were orders of magnitude more luminous than the main-sequence stars (also called dwarf stars) at roughly the same surface temperature. In 1938, Ernst Öpik first considered stellar models that were not homogeneous by entertaining the idea that the interior of stars is generally not well mixed (Öpik,1938). This meant that the core of a star could run out of hydrogen to burn and become inert. He predicted the core would rapidly contract and release gravitational energy instead, expanding the envelope that is still hydrogen rich. This process, i.e., the depletion of hydrogen in the core and the accompanying structural changes of the star, would naturally produce the giant stars in the HRD. Once computers entered the scene in the 1950s, stellar models could be computed using numerical techniques. In those days, machines with the computational power of today’s pocket calculators had the size of whole classrooms, but the physical processes occurring in stars could be programmed (with tape or punch cards! 5 ), and detailed stellar-evolution models were produced as early as 1952. As computers became more and more part of the standard equipment of the astrophysicist, innovative techniques in the implementation of the differential equation solver were developed, most notably with the contribution by Louis Henyey (Henyey et al.,1959). The so-called Henyey method is still used today for calculating stellar-evolution models, and uses a trial solution and a generalized Newton-Raphson root-finding algorithm to iteratively approach a better solution to the differential equations. As atomic and nuclear data kept improving, in part driven by research on nuclear weapons during World War II and the Cold War, and computers kept getting faster exponentially, decade after decade saw substantial progress in the calculation of stellar models. Rather than computing single models of one particular mass and composition, whole sets, or populations, of stellar models could be explored. This enabled to study the evolution of clusters of stars, where the assumption can be made that all stars are born more or less at the same time and more or less with the same composition. Alternatively, having a lot of computing power allows to study the impact of varying certain parameters that describe physical processes, such as the efficiency of internal mixing, the rate at which certain nuclear reactions occur, or the strength of mass loss of the star during its life (see Sect. 1.1.5). 5 Anecdotally, Rudolf Kippenhahn recalls that, from the patterns in the flashing lights and noises the computer made, he and his collaborators could infer which specific part of the star it was calculating (Gingerich,1978). 10 INTRODUCTION Today, although supercomputers are used to calculate grids of tens or even hundreds of thousands of stellar models, astrophysicists do not need supercomputers hosted at academic or governmental institutions to produce detailed models. Most stellar-evolution codes are designed to run on basically any computer, and it will take only on the order of hours to complete the integration of, for example, the life of the Sun, taking around ten thousand integration steps. This is in stark contrast with the situation of yesteryear, as Martin Schwarzschild noted in 1958: A person can perform more than twenty integration steps per day ... so that for a typical single integration of, say, forty steps, less that two days is needed. — Schwarzschild (1958). If one would attempt to compute the stellar models contained in this thesis without the help of modern computers, instead relying on manual integration like in the times of Schwarzschild, this thesis would require an estimated 8000 years to complete, and no funding agency would want to back a grant for that long. This goes to show how powerful computers have become in just seventy years or so. Progress in computer science does not stand still however, and, with machine learning technologies finding applications in the natural sciences, including astrophysics, maybe the field of stellar structure and evolution, in particular the computation of stellar models, will again be revolutionized soon. 1.1.1 The equations The study of stellar structure and evolution revolves around the solution of a set of partial differential equations, each one encoding a fundamental physical principle. They are: (i) Conservation of mass, (ii) Conservation of momentum, (iii) Conservation of energy, (iv) Transport of energy, and (v) Changes in composition. Under general conditions, the set of equations associated to these principles is not analytically solvable since we are dealing with coupled, non-linear, partial differential equations in three spatial dimensions plus a time dimension. However, by observing that stars are, generally speaking, round balls, we can make use of SINGLE STARS 11 spherical symmetry to cast the equations in the following form: 𝜕𝑀 𝜕𝑟 =4𝜋𝑟2𝜌, (1.5a) 𝜕𝑃 𝜕𝑟 =−𝜌𝜕𝛷 𝜕𝑟 =−𝜌𝐺𝑀 𝑟2,(1.5b) 𝜕𝐿 𝜕𝑟 =4𝜋𝑟2𝜀, (1.5c) 𝜕𝑇 𝜕𝑟 =−𝜌𝐺𝑀 𝑟2 𝑇 𝑃∇,(1.5d) 𝜕𝑋𝑖 𝜕𝑡 =∑︁ 𝑗 𝑟𝑗𝑖 −∑︁ 𝑘 𝑟𝑖𝑘 −1 𝜌𝑟2 𝜕 𝜕𝑟 𝜌𝑟2𝐷𝜕𝑋𝑖 𝜕𝑟 .(1.5e) We now have only one spatial dimension, 𝑟 , and the time, 𝑡 , left as independent variables. In these equations, several functions are defined. We have the set consisting of 𝑀=𝑀(𝑟, 𝑡) , which is the mass at time 𝑡 enclosed in a sphere of radius 𝑟 , and 𝑃=𝑃(𝑟, 𝑡), 𝑇 =𝑇(𝑟, 𝑡) and 𝐿=𝐿(𝑟, 𝑡) , which are, respectively, the pressure, temperature and luminosity of the stellar material at the radius 𝑟 and time 𝑡 . These four functions specify the structure of the star at each point in time. We have also 𝜌=𝜌(𝑃, 𝑇, 𝑋𝑖) , the density of the material, which is computed via an equation of state, and the potential 𝛷=𝛷(𝑀, 𝑟) , which is the gravitational potential associated with the mass configuration of the star, and can be solved from Poisson’s equation, ∇2𝛷= 4 𝜋𝐺𝜌 . The function 𝜀=𝜀(𝜌, 𝑇, 𝑋𝑖) encodes energy generation by nuclear reactions as well as thermodynamic heating or cooling by compression or expansion, and ∇ is the logarithmic derivative of temperature to pressure required to transport the luminosity 𝐿(𝑟) through the star. The last equation with 𝜕𝑋𝑖 𝜕𝑡 is actually a set of multiple equations that describes, for each element 𝑖 , the transmutation of 𝑖 into different elements 𝑘 at a rate of 𝑟𝑖𝑘 , while other elements 𝑗 are changed into 𝑖 at the rate 𝑟𝑗𝑖 . The third term represents internal mixing of the stellar material, occurring with diffusion coefficient 𝐷 . A full solution of all equations, giving the functions 𝑀(𝑟, 𝑡), 𝑃(𝑟, 𝑡), 𝑇 (𝑟, 𝑡), 𝐿(𝑟, 𝑡) and the elemental species 𝑋𝑖(𝑟, 𝑡) , with appropriate initial and boundary conditions is then called a stellar model. Equations of continuity and hydrostatic equilibrium The first two equations are results of classical physics and were well understood already at the start of the twentieth century. The first equation, Eq. (1.5a) , called the continuity equation, simply states that the mass, 𝑑𝑀 , contained in a thin shell of radius 𝑟 is the product of the density, 𝜌 , of the shell, its area, 4 𝜋𝑟2 , and its thickness, 12 INTRODUCTION Table 1.1: Final assembly reactions in the pp chain of hydrogen fusion into helium. pp-I pp-II pp-III 3He + 3He −−−→ 4He + 2 1H3He + 4He −−−→ 7Be + 𝛾 7Be + e–−−−→ 7Li + 𝜈e7Be + 1H−−−→ 8B + 𝛾 7Li + 1H−−−→ 24He 8B−−−→ 8Be + e++𝜈e 8Be −−−→ 24He 𝑑𝑟 . Equation (1.5b) asserts that forces arising from the gravitational potential, 𝛷 , are balanced by the internal pressure gradient coming from the gas making up the star and the radiation it produces. We call this situation hydrostatic equilibrium. In most of the life of a star, hydrostatic equilibrium is satisfied, or very nearly so. However, for example, when a massive star collapses at the end of its life, we speak of hydrodynamic evolution, and an acceleration term −𝜌 𝜕2𝑟𝜕𝑡2 must be added to the right hand side of Eq. (1.5b) . Additionally, Eq. (1.5b) ignores the fact that, at the surface of the star, radiative forces can launch a stellar wind, but this effect is usually taken into account by incorporating mass loss in the models (see Sect. 1.1.5). Equations of energy conservation & transmutation The third equation, Eq. (1.5c) , is the equation of conservation of energy and is intimately connected with Eq. (1.5e) specifying the transmutation of elements into others. As discussed in the introductory paragraphs, it was found in the early twentieth century that the main source of stellar energy is nuclear. On the main sequence, stars convert hydrogen into helium. In low-mass stars ( 𝑀≲ 2 𝑀⊙ ), this process is achieved by direct assembly, also called the proton-proton or pp chain. Two hydrogen nuclei, 1H , which consist of just a single proton, combine while one of them changes into a neutron, forming deuterium, and emit a positron e+ and a neutrino 𝜈e:1H+1H−−−→ 2H+e++𝜈e(1.6) Adding then another proton forms 3He and emits a photon, 𝛾: 2H+1H−−−→ 3He +𝛾(1.7) Final assembly to 4He then occurs by either combining two 3He nuclei (pp-I) or, with 4He already present, going via beryllium and lithium (pp-II) or boron (pp-III), as presented in Table 1.1. Each of these reactions releases an amount of energy, 𝜀 , to be used to hold up the star against its own gravity. SINGLE STARS 13 Figure 1.4: CNO cycle of hydrogen burning into helium. In massive main-sequence stars this is the dominant process from which nuclear energy is liberated. Adapted from Marchant (2018). The second way to fuse hydrogen is through the CNO cycle, which is the dominant process in intermediate-mass (2 𝑀⊙≲𝑀≲ 8 𝑀⊙) and massive stars ( 𝑀≳ 8 𝑀⊙ ). Instead of assembling a helium nucleus from scratch by starting from single protons, in this process we start with already present carbon and nitrogen nuclei, and add protons one at a time until a helium nucleus splits off. This happens in two distinct cycles, going through both stable and unstable isotopes of carbon, nitrogen, oxygen and fluorine as shown in Fig. 1.4. An important property of this cycle is that the proton capture of nitrogen, 14N+1H−−−→ 15O+𝛾(1.8) occurs the slowest, meaning that the CNO cycle, while burning hydrogen to helium, converts a fraction of the carbon and oxygen into nitrogen. Observing an enhanced nitrogen abundance on the surface of a star is therefore a good indication that nuclearly processed material of the core has been brought to the surface by internal mixing processes (see Sect. 1.1.5). When a star leaves the main sequence as its hydrogen runs out in the core, the next nuclear fusion reaction available is the combination of three helium nuclei to a carbon nucleus, which is called the triple alpha process: 34He −−−→ 12C+𝛾. (1.9) 14 INTRODUCTION This reaction is crucially dependent on the stability of the intermediary product beryllium: 4He +4He −−−→ ←−−− 8Be.(1.10) Beryllium-8 is unstable, and decays back to two 4He nuclei with a half-life of only about 8 . 2 × 10 −17 s. Luckily, the density and temperature in a stellar core where these reactions take place is high enough for the third helium nucleus to encounter the 8Be particle before it decays again, and fuse to 12C. The final reaction: 8Be +4He −−−→ 12C∗−−−→ 12C+𝛾, (1.11) however passes through an excited resonance state of carbon, 12C* . Predicted by Fred Hoyle in the early fifties (Hoyle,1954), it is now called the Hoyle state. The state was experimentally found not much later at Caltech (Dunbar et al.,1953; Cook et al.,1957), but it took until 2011 to theoretically calculate its existence (Epelbaum et al.,2011). The fact that stellar nuclear burning has to go through 12C* and that this state exists almost exactly at the sum of energies of 8Be and 4He is of considerable interest to philosophers in context of the anthropic principle and the fine-tuning of the universe (Kragh,2010). If this state did not exist, or had its properties altered by slightly different fundamental physical constants, the production of elements heavier than helium could have been radically different, or stalled altogether, and life as we know it might not have emerged. In massive stars heavier than about 8 𝑀⊙ , thermonuclear burning proceeds from 12C to heavier elements like 16O , 20Ne , 28Si , and finally 56Fe . All these steps use multiple intermediary nuclei and interaction channels, much too large to list in this thesis. Instead we refer to works on stellar nucleosynthesis, such as Clayton (1983). We do however note the reaction 12C+4He −−−→ 16O+𝛾, (1.12) which has a particularly hard to determine reaction rate (deBoer et al.,2017). Because this rate establishes the precise carbon to oxygen ratio after helium burning, it has important implications for the final fate of the massive star (e.g., on the masses of black holes, Farag et al.,2022). Figure 1.5 shows the binding energy per nucleon of the most stable isotope of each mass number 𝐴 . It shows that the first step in the fusion chain, hydrogen burning to helium, is the most energetically favorable reaction since the jump in binding energy is by far the largest. This translates in the fact that a star will spend the vast majority of its life burning hydrogen in its core. Each progressive burning stage releases less energy, which means these stages will take less time to complete. As soon as 56Fe is formed in the core of a massive star, nuclear fusion stops as it is no longer energetically favorable for the star. The whole machinery of producing energy to maintain a pressure gradient which counteracts gravity falls apart like a house of cards. This announces the dramatic end of the life of the massive star with a supernova (see Sect. 1.1.4). SINGLE STARS 15 Figure 1.5: Binding energy per nucleon of the most stable isotope per mass number 𝐴 . Indicated are several important isotopes in stellar nucleosynthesis, as well as silver-108, gold-198 and uranium-235, the latter being the prime fissile material in modern nuclear power plants. The vertical line indicates iron-56, which is the most efficiently bound nucleus. Nuclear data from Huang et al. (2021); Wang et al. (2021). In many of the nuclear reactions of the pp chain and the CNO cycle referred to above, neutrinos are produced as a result of interaction with the weak nuclear force. Furthermore, in burning stages of carbon and later, the conditions in the star are such that photons are energetic enough to create an electron-positron pair, which occasionally turn into a neutrino-antineutrino pair. However, once produced, neutrinos are notoriously non-reactive with protons and electrons, and therefore do not deposit their energy into the gas of the star. Stars thus lose energy to neutrino production, which can be substantial in the later phases of massive-star evolution. Additionally, gas could be heating, cooling, contracting or expanding, all of which affects the energy balance in the stellar layers. Rather than hiding all possible energy terms in the general function 𝜀 , the energy equation is more explicitly written as: 𝜕𝐿 𝜕𝑟 =4𝜋𝑟2𝜀nuc −𝜀neu −𝑐𝑃 𝜕𝑇 𝜕𝑡 +𝛿 𝜌 𝜕𝑃 𝜕𝑡 ,(1.13) where the four terms on the right hand side are, in order, the nuclear energy released from fusion, energy losses by neutrinos, thermodynamic heating or cooling with specific heat capacity 𝑐𝑃 and volumetric expansion or compression with 𝛿=−𝜕ln 𝜌/𝜕ln𝑇. 16 INTRODUCTION Equation of energy transport The transport of energy within the star, specified by the fourth equation, Eq. (1.5d) , was first thought to be mediated solely by convective blobs of gas. Convection happens all around us, take for example boiling water on a stove. The intense heat near the bottom of the pot creates bubbles that rise to the surface. The same effect also powers ocean currents, weather patterns in Earth’s atmosphere, and causes tectonic activity of Earth’s crust as the mantle below convects heat from the core to the surface. However, it was put forward by Karl Schwarzschild that radiation can also play a significant role in the transportation of energy from the core to the outer layers of a star (Schwarzschild,1906). This idea was put on solid theoretical grounds some time later by Eddington (1916), where he also discussed the sources of opacity inside stellar material. The type of energy transport, radiative or convective, is an important factor in determining the structure of a star. Considering first radiative transport, which is a diffusive process, Fick’s law applies (Fick,1855), which, when applied to radiation, gives the radiative flux: Frad =−4𝑎𝑐 3 𝑇3 𝜅𝜌 ∇𝑻.(1.14) Here the opacity, 𝜅 , of the material is used, and 𝑎= 8 𝜋5𝑘4 B/( 15 𝑐3ℎ3) is the radiation density constant, defined using the fundamental constants of Boltzmann, 𝑘B , and Planck, ℎ . Using then that the luminosity is 𝐿(𝑟)= 4 𝜋𝑟2Frad , from Eq. (1.14) we have for a spherically symmetric configuration: 𝜕𝑇 𝜕𝑟 rad =−3 16𝜋𝑎𝑐 𝜅𝜌𝐿 𝑟2𝑇3,(1.15) which, dividing by the equation of hydrostatic equilibrium, defines a radiative temperature gradient ∇rad: ∇rad ≡dln𝑇 dln 𝑃rad =3 16𝜋𝑎𝑐𝐺 𝜅𝐿𝑃 𝑀𝑇4.(1.16) If energy is transported fully by radiation, the temperature within the star will thus vary with pressure as given by Eq. (1.16). Convective energy transport will occur if gas motion is unstable against adiabatic displacement. Put simply, if a blob of gas rises (or sinks) in the star and finds itself in an environment that makes it rise (or sink) even more, large-scale convective bubbles will be created. Defining the adiabatic temperature gradient at any point in the star: ∇ad ≡dln𝑇 dln 𝑃𝑠 ,(1.17) SINGLE STARS 17 the subscript meaning to take the derivative with constant specific entropy, 𝑠 , a stellar layer will be unstable to convection if: ∇ad <∇rad.(1.18) This is called the Schwarzschild criterion of convection (Schwarzschild,1906). The Belgian astrophysicist from Liège, Paul Ledoux, modified this criterion by also considering stabilizing molecular gradients (Ledoux,1947). He realized that when a convective blob moves to a region where the chemical composition is different, when accounting for this in the stability criterion, an extra term needs to be added that depends on the composition gradient. The Ledoux criterion for convection is: ∇ad <∇rad −𝜑 𝛿∇𝜇,(1.19) with the definitions 𝜑≡𝜕ln 𝜌/𝜕ln 𝜇 , 𝛿≡ − 𝜕ln 𝜌/𝜕ln𝑇 and ∇𝜇≡𝜕ln 𝜇/𝜕ln 𝑃 to write the stability condition concisely. The mean molecular weight, 6𝜇 , is defined as the amount of mass the average free particle has in units of 𝑚u: 𝜇=𝜌 𝑛𝑚u,(1.20) with 𝑛 the number of free particles per unit volume in the mixture. A higher 𝜇 thus means that particles have more nucleons on average. Conversely, a lower 𝜇 means that particles are lighter, or that the mixture is more ionized, which increases the number of free electrons. Looking at Eqs. (1.16) and (1.19) , we see that convection will set in if either the luminosity to mass ratio, 𝐿/𝑀 , or the opacity, 𝜅 , is particularly high, while steep composition gradients, ∇𝜇 (e.g., left behind by nuclear burning), have stabilizing effects. Computing the actual temperature gradient, ∇ , of the star in a convective zone requires a theory of turbulent motion. To emphasize, ∇ for a convective region is not a result of the stellar-structure equations, but is rather an input required to describe the temperature gradient of a convective region. Describing convection is hard, as it is a turbulent flow that is chaotic. In principle, the full Navier-Stokes equations need to be solved in three dimensions, after which parameterizations to one dimension can be attempted. We can however make the approximation that convection will occur on a typical length scale, and consider that all convective bubbles travel that length and then promptly deposit their excess heat into the environment. This set of assumptions defines the mixing length theory (MLT) of convection, and was pioneered by Erika Böhm-Vitense as she first described convection in the Sun (Vitense,1953;Böhm-Vitense,1958). MLT allows to compute a one-dimensional efficiency of convection, 𝜁 , by comparing the convective heat 6 This is, in my opinion, another confusing misnomer, as, in the context of stellar structure, this has nothing to do with weight nor with molecules. I prefer mean particle mass. 24 INTRODUCTION Figure 1.10 : Crab Nebula as imaged by the James Webb Space Telescope. It is the remnant of a supernova that erupted in 1054. Image credit: NASA, ESA, CSA, STScI, T. Temim. collapses under its own weight, causing a violent supernova accompanied by a neutrino outburst (however, a supernova can be avoided altogether, see below). In a supernova, part of the mass of the massive star progenitor (including its nuclearly processed material as well as heavy elements produced in the supernova itself) is ejected into the interstellar medium, making massive stars the primary engine of chemical enrichment in galaxies (Burbidge et al.,1957). Some of the best studied supernovae are the ones that occurred in our galactic neighborhood. In 1054, Chinese and Japanese astronomers detected a “guest star” (Stephenson & Green,2003), which evolved into the Crab Nebula we see today (Fig. 1.10). More recently, SN 1987A erupted in February 1987 in the Large Magellanic Cloud (LMC), one of the satellite galaxies of our Milky Way. It was the first nearby supernova studied with modern telescopes and instruments, including neutrino detectors, and thus it provided excellent observational data about the physics of core-collapse supernovae. The remnant of a core-collapse supernova is either a neutron star or black hole (BH), the former being a dense ball of nuclear material held up by neutron degeneracy pressure and strong nuclear interactions, the latter being a singularity that can currently only be described by Einstein’s theory of general relativity (Einstein, 1915;Schwarzschild,1916). The main parameters that determine what remnant (if any) will be left behind after a supernova is the mass of the stellar core, and the metallicity of the material, i.e., the fraction of mass in metals. Metallicity is important as it drives mass loss (see Sect. 1.1.5) and so it influences the mass of the carbon-oxygen core by the time helium is depleted. Furthermore, the explodability of the core depends on the complex interaction between neutrinos and the core structure itself that is beyond the scope of this thesis (see Heger SINGLE STARS 25 et al. 2003,2023, for summaries). Finally, if the core is too hard to explode, a supernova may not occur at all, and the star collapses directly into a BH without launching much material into interstellar space. The main result from all these considerations is that it is harder to form BHs at solar metallicity than at low metallicity, where carbon-oxygen cores are heavier. Those that do form are due to supernova fallback, where material that is launched in the supernova shock drops back onto the proto-neutron star and tips it over the mass limit to become a BH. However, the orbital characteristics of Cygnus X-1 and VFTS 243 indicate that BH formation through direct collapse is likely what occurred in these systems, and both are located in the high-metallicity environments of the LMC and the Milky Way (Mirabel & Rodrigues,2003;Shenar et al.,2022). We draw this conclusion because the eccentricity of the orbit is very low, which suggests that there was very little mass loss during the formation of the BH, thus favoring a direct collapse over a supernova. 1.1.5 Open questions in massive-star structure and evolution The physical processes outlined above give, generally, a good picture of massivestar structure and evolution. We can explain the general trends in the life of a massive star, and understand the chemical enrichment by heavy elements resulting from their winds and supernovae. We have detected both neutron stars and black holes left behind by massive stars. The first neutron star was discovered by Jocelyn Bell using radio observations (Hewish et al.,1968), while the first convincing black hole candidate was Cygnus X-1, identified independently by Bolton (1972) and Webster & Murdin (1972) using a combination of X-ray observations from the UHURU space telescope and measurements from radio antennas. Despite these successes, there are still multiple unsolved discrepancies between theory and observation regarding stellar evolution, particularly so in the massivestar regime, some of which we outline below. Mass discrepancy A major problem present in the literature is the so-called mass discrepancy problem. For single stars, there are two main ways of calculating its mass. One involves measuring the luminosity and effective temperature of the star, and matching them to the position of evolutionary models in the HRD; the best fit giving the evolutionary mass. The other requires a measurement of the surface gravity 𝑔=𝐺𝑀/𝑅2 through spectroscopic data, which, if the radius can be estimated from the luminosity and effective temperature, gives the spectroscopic mass of the star. Many studies however find a systematic mismatch of both mass calculations (Herrero et al.,1992;Mokiem et al.,2007;Martins et al.,2012; Tkachenko et al.,2020;Mahy et al.,2020). At this moment it is unclear whether 26 INTRODUCTION the problem lies in the accuracy of determining surface gravities, which can be challenging for O stars, or if the evolutionary tracks suffer from uncertain physical modeling. Mass loss Massive stars have a strong stellar wind on the main sequence (of order 10 −6𝑀⊙yr−1 or more). This is in contrast with low-mass stars, which lose negligible mass on the main sequence (e.g., the Sun has a mass loss rate of merely 10−14 𝑀⊙yr−1). Since stellar-structure and evolution theory does not (a priori) predict the mass-loss rates of stars, other theories or experiments need to be used. The first evidence of stellar winds in massive stars came from Morton (1967), who used measurements in the ultraviolet to infer outflows from three stars. 9 Later, it was proposed by Lucy & Solomon (1970) and Castor et al. (1975) that such outflows are driven by the absorption and scattering of radiation by elements (in particular metals) in the outer layers of the star. The radiation coming from the interior therefore deposits its momentum in these atoms and drives them out. Actually quantifying how strong a radiation field can accelerate atoms this way is a difficult problem. Vink et al. (2001) calculated the interactions of a large number of photons with over 10 5 different atomic lines of elements ranging from hydrogen to zinc. More recent calculations however predict lower mass-loss rates (Sundqvist et al.,2019;Björklund et al., 2021,2022). For example, a 15 𝑀⊙ star with luminosity of log(𝐿/𝐿⊙)= 4 . 5 and a surface temperature of 30000K has a predicted mass-loss rate by Vink et al. (2001) of ¤ 𝑀≈ 1 × 10 −8𝑀⊙yr−1 , while Björklund et al. (2022) predict ¤ 𝑀≈ 2 × 10 −9𝑀⊙yr−1 , almost an order of magnitude lower. It has become apparent that mass loss, both on the main sequence and beyond, can significantly change the evolution of a massive star. For example, if mass loss during the red supergiant phase is particularly strong and blows away the whole hydrogen rich envelope, the hot helium core can become exposed to outside observers. If it does, such objects are called Wolf-Rayet (WR) stars, after Charles Wolf and Georges Rayet discovered they had strong and broad emission lines in their spectra, which are direct consequences of an optically thick wind. It is however unclear if WR stars can be produced from a supergiant wind alone (e.g., Langer et al.,1994;Meynet & Maeder,2005), or if binary interaction is needed to take away the hydrogen layers (e.g.,Pauli et al.,2022). Therefore, understanding the evolutionary pathway and observable quantities of a star depends on knowing the precise mass-loss rates that massive stars have during the various phases of their lives. 9 Since ultraviolet light is largely absorbed by the atmosphere, he needed a rocket-mounted spectrograph to get to altitudes of 200km to get sufficient signal. SINGLE STARS 27 Rotation and internal mixing In stating the equations of stellar structure in Eqs. (1.5) , we assumed no forces other than self-gravity and internal pressure gradients act on the stellar material. However, as soon as a ball of gas is spinning about an axis, another term should be considered in the balance of forces of Eq. (1.5b) . To sustain rotational motion, a centripetal acceleration is needed, which will be supplied by the gravitational attraction. Conversely, when viewed in the frame of reference of the rotating star, a centrifugal force appears that points away from the rotation axis. In the equation of hydrostatic equilibrium, we write it as: 𝜕𝑃 𝜕𝑟 =−𝜌𝜕𝛷 𝜕𝑟 −𝜌𝑟 𝛺2,(1.24) where 𝛺 is the angular velocity about the axis of rotation. A major issue with Eq. (1.24) is that it is only valid for in the equatorial plane of a star. By introducing rotation, we necessarily added an independent spatial coordinate, typically chosen to the the polar angle, 𝜃 . In 2D then, the equation of hydrostatic equilibrium, including rotation, reads: ∇𝑃=−𝜌∇𝛷−1 2𝜌𝛺2∇(𝑟sin 𝜃)2.(1.25) Note that the gravitational potential 𝛷 is also no longer spherically symmetric, since the mass distribution is deformed. In principle this means all equations of stellar structure need to be solved in two dimensions, 𝑟 and 𝜃 . While possible, this usually means sacrificing solving the time dependent problem, so no evolution is simulated, as is done by Roxburgh (2004) and Espinosa Lara & Rieutord (2013). Deupree (1990,1995) and recently, Mombarg et al. (2023) succeeded in computing two-dimensional stellar-structure models while also computing the evolution along the main sequence. An alternative treatment of rotation exist however, as Kippenhahn & Thomas (1970) provided a method to include rotation into 1D stellar-structure calculations, which we will introduce in Sect. 3.2. Another uncertain physical process associated to rotation is rotationally induced mixing. Very early on, it was identified by Edvard Hugo von Zeipel that a rotating star cannot be in hydrostatic and radiative thermal equilibrium at the same time (von Zeipel,1924). This arises since, in hydrostatic equilibrium, the amount of energy that radiation can transport is proportional to the magnitude of the local effective gravity, Frad ∼𝒈eff ≡ −∇𝛷−1/2𝛺2∇(𝑟sin𝜃)2 . This relation is called von Zeipel’s theorem. In radiative thermal equilibrium, the equation of conservation of energy can be written as: ∇·Frad =𝜀nuc (1.26) If a star is in hydrostatic equilibrium, layers of constant density, temperature, pressure and composition align with surfaces of constant effective potential 𝛹= 𝛷+1 2𝑟2𝛺2 (von Zeipel,1924). Von Zeipel then showed that, on an equipotential surface, the left-hand side of Eq. (1.26) is a function of the polar angle 𝜃 , while the 28 INTRODUCTION Figure 1.11: Circulation currents in a slowly, uniformly rotating star. These currents balance out the non-uniform radiative energy transport as a result of von Zeipel’s theorem. Figure from Sweet (1950). right-hand side is constant, since 𝜀nuc =𝜀nuc(𝜌, 𝑇, 𝑋𝑖) is a function of microphysical properties. In hydrostatic equilibrium then, under general conditions, radiative equilibrium cannot be satisfied. It is hard to imagine that the nuclear fusion processes conspire with stellar rotation to keep isobars and surfaces of constant temperature aligned with the effective potential. Instead, Eddington (1925) and Vogt (1925) independently proposed a solution to this problem by considering meridional circulation currents. These currents operate to carry away the excess energy generated in one region and deposit them in another so as to maintain radiative thermal equilibrium and hydrostatic equilibrium. Later, Sweet (1950) calculated the circulation velocities in the case of a slowly, uniformly rotating star, and produced the picture in Fig. 1.11. Such flows are now called Eddington-Sweet circulations. While the origin of meridional circulation is reasonably well understood, the mixing of stellar material associated with these currents is not. From Fig. 1.11, it is clear that material near the core, which is altered by nuclear processes, can be transported up to the surface. However, this picture is a rather simplified view of the detailed physics that happens. Zahn (1992) shows that turbulent viscosity is at the root of the mixing process, rather than the meridional circulation, and derives it can only be efficient for very fast rotators. Furthermore, since fused elements near the core (like helium, carbon, etc.) are heavier than untouched hydrogen in the envelope of the star, there is a molecular weight gradient that inhibits these SINGLE STARS 29 circulation currents. Different groups implement this interaction in different ways (compare, e.g.,Heger et al.,2000, and Meynet & Maeder,2000), and produce different results. While both models find elements like nitrogen are enhanced near the surface, only the models of Meynet & Maeder (2000) produce significant helium enrichment for fast-rotating, massive main-sequence stars (unless the stars rotate so fast they evolve chemically homogeneously, see Sect. 1.4.1). Using the 2D treatment, Mombarg et al. (2023) implemented the theory of Zahn (1992) and showed in their models that meridional flows are indeed ineffective at transporting chemical elements from the core to the surface. Convective-core overshooting In massive stars, the core is unstable against convective motion of material, while in the envelope, radiative transport of energy occurs (see Fig. 1.6). At the edge of the core is thus a convective boundary, beyond which convective blobs should not exist anymore. It is however not unreasonable to think that the inertia of such blobs carry them past the boundary. This process is called convective-core overshooting. The question is how far convective elements can penetrate into the radiative envelope (if at all), and in so doing mix the material that would otherwise not participate in the nuclear burning of the core. The net effect of overshooting is the lengthening of the life of the star, since the more nuclear fuel a star has available, the longer it can maintain equilibrium against gravity. Furthermore, stronger overshooting increases the mass of the helium core, which has important ramifications for the further evolution (and fate) of the massive star. Just like describing convection itself (see Sect. 1.1.1), theoretically deriving the properties of core overshooting is very hard. We note here the work of Canuto (1992), who derived a theory of convection that naturally includes a description of overshooting, but it has not found wide-spread use in stellar-structure calculations. Instead, overshooting is commonly parameterized by a dimensionless number, 𝛼ov , which signifies how many pressure scale heights, 𝐻𝑃 , convection extends beyond the convective boundary set by the Ledoux criterion, with 𝐻𝑃=−dln 𝑃 d𝑟−1 .(1.27) Using observations of clusters of stars, Brott et al. (2011) calibrated the overshoot parameter to 𝛼ov = 0 . 335 for massive stars in the LMC, while Schroder et al. (1997) constrained it to 𝛼ov ∈ [ 0 . 24 , 0 . 32 ] for intermediate-mass stars. More recently, Tkachenko et al. (2020) found values consistent with 𝛼ov = 0 . 4 for most stars in a sample of galactic eclipsing binaries. Another way to constrain overshooting is by employing the theory of asteroseismology, the study of stellar oscillations (Aerts et al.,2010). Just as with Earth’s seismology, the surface oscillations of stars allow us to peer inside them and gauge 30 INTRODUCTION the mass of the convective core (see e.g.,Johnston,2021 for a list of such convective core measurements), and thus constrain the strength of the overshoot (e.g.,Aerts et al.,2003;Moravveji et al.,2016;Burssens et al.,2023). Finally, with increasing efficiency of hardware and innovation in computational algorithms, in particular the use of CPU parallelization and GPU computing, more and more three-dimensional simulations of convection in stellar cores are performed (e.g.,Browning et al.,2004;Cristini et al.,2017;Edelmann et al.,2019; Vanon et al.,2023;Andrassy et al.,2024). Such simulations could provide a direct mass and luminosity dependent measure of the overshoot parameter, as well as the temperature gradient in the core. 1.2 Binary stars The historical introduction in this section is mainly based on Aitken (1935). While our own Sun is a lonely star, it is not impossible for stars to live and evolve near each other. One famous example featured in pop culture is the stars hosting the planet Tatooine in the Star Wars franchise. As Luke looks over the desert at sunset (suns-set? Or just sunsets?), two reddish stars are sinking toward the horizon, see Fig. 1.12. Thanks to Newton’s synthesis of motion and gravitation (Newton,1687,1999), we know that the forces that keep such stars revolving around each other are the same ones that make the Earth move around the Sun, or that make apples fall from trees. The study of binary stars, like regular single stars, is not a recent affair. As early as the Greek and Roman period, Ptolemy of Alexandria used the term double star to describe 𝜈 Sagittarii, which are two stars separated by only about half of the Moon’s apparent diameter (which is around thirty arcminutes or half a degree). The first modern observations of double stars are usually attributed to the Italian Jean Baptiste Riccioli in 1650, who found the star Mizar in Ursa Major had two bright components. Some historians challenge this credit, and letters from Galileo to Benedetto Castelli indicate it was Castelli who first resolved Mizar into two components as early as 1617 (Ondra,1999). Christiaan Huygens and Robert Hooke, among many others, made similar such observations in the seventeenth century. Initially, there was no specific interest in these objects, as it was thought it was mere coincidence that two stars happen to align when viewed from Earth. It was in the mid-eighteenth century that this perspective changed. Based upon general probability arguments, John Mitchell argued in 1767 that stars that appear close to each other are indeed very likely to physically be near each other too (as opposed to scattered randomly throughout the galaxy and aligning by chance), and thus compose a stellar system. This realization sparked William Herschel and BINARY STARS 31 Figure 1.12: Binary star sunset over the planet Tatooine in the movie Star Wars Episode IV: A New Hope. © Lucasfilm Ltd., Walt Disney Studios. Reproduced under the Fair Use Act. Christian Mayer to systematically look for binary stars, and by 1782, hundreds of double stars were cataloged (Herschel & Watson,1782). At the start of the nineteenth century, it was Herschel who made the formal distinction between double stars and proper binary systems, by requiring for the latter that the two stars are “united by the bond of their own mutual gravitation toward each other” (Herschel,1802). One year later, he published the first set of measurements that irrevocably proved that Newton’s laws determine the motion of some of the observed double stars (Herschel,1803). Up to this point, only visual binaries were discovered. These are binary systems that are resolvable by direct imaging with telescopes, where the size of the main mirror (and the quality of the atmosphere) determines how close stars can be distinguished. However, with the advent of spectroscopy, in 1887, Antonia Maury at Harvard noticed that the dark lines in the spectra from Mizar A appeared single sometimes, and double other times, see Fig. 1.13. Together with Pickering, they concluded that one component of the visual binary is in itself another binary system (Pickering,1890). Based upon the Doppler-Fizeau principle, in aspectroscopic binary (SB) then, as a star alternately moves toward and away from the observer, a spectral line gets blueshifted to shorter wavelengths as the star is approaching, or redshifted to longer wavelengths when the star is receding. The effect we experience in daily life when hearing an ambulance, car or airplane pass by is analogous. The position of the line correlates with the radial velocity (RV) of the star, i.e., the speed with which it approaches or recedes from the observer. When the spectral line doubled, the Harvard team thus saw one star receding and the other star approaching the Earth, while when the line was single, both 32 INTRODUCTION Figure 1.13: Two photographic plates showing the line splitting of Mizar A’s spectra. On the left panel, a double line can be distinguished, while on the right, the line appears single. Images reproduced from Ondra (1999). stars were moving tangentially to the plane of the sky (relative to the motion of the binary system as a whole). Only months later, Hermann Carl Vogel saw the spectra of Algol move in correlation with the periodic dimming of the star (Vogel, 1890). The theory by John Goodricke that a dark(er) object passed in front of the star (Goodricke,1783) thus turned out to be correct, and the class of eclipsing binaries had its first discovery. In the late nineteenth century, methods for characterizing orbits were developed. Thorvald N. Thiele published a method for computing the orbital elements of a visual binary by describing the apparent ellipse the orbit draws on the sky (Thiele, 1883). This was later developed on by Robert T.A. Innes (van den Bos,1926), resulting in what is now called the set of Thiele-Innes constants describing visual binaries. While the equation describing the RV of spectroscopic binaries can be derived from elementary geometry, actually fitting observed RVs was a challenge in the past. Today, we can use numerical methods and computing power to find a well fitting orbit, but in the early twentieth century several algorithms involving iteration were composed, e.g., from Lehmann-Filhés (1894), Schwarzschild (1900) and Russell (1902), among others. From the mid-twentieth century, observations with X-ray telescopes uncovered another way to detect binaries. When a star is donating material to a nearby companion, an accretion disc can form around the accretor where the in-falling matter heats up tremendously. This hot gas then emits high-energy radiation in the form of X-rays, for which a specific type of telescope is needed to detect them. The very first such instruments were rocket mounted to get above the Earth’s atmosphere, but subsequent designs were full space telescopes like UHURU (Giacconi et al.,1971), Chandra (Weisskopf et al.,2000), XMM-Newton (Jansen et al.,2001), and NuSTAR (Harrison et al.,2013). In order to reach the temperatures required to generate X-rays, the accretor needs to be very compact so the in-falling matter gains significant energy. Elevated X-ray emission from an object is thus a telltale sign of a binary system that is exchanging mass, and it also provided the first indirect evidence of black holes existing in nature. BINARY STARS 33 One of the main benefits of studying binary systems is that they provide a direct measurement of the masses of the stars through Kepler’s laws. Such dynamical mass estimates usually suffer from much less uncertainty than gauging the mass from evolutionary models or having to rely on atmospheric modeling as is done for single stars. It is still far from a trivial task, as, for visual binaries, computing the masses requires knowing the distance to the system as well as the motion of both components, rather than their relative motion. For spectroscopic binaries, a major unknown is the inclination of the orbit (i.e., the angle of the orbital plane to that of the sky), and allows only to estimate a lower limit on the mass of the stars. Given these problems, eclipsing binaries are particularly useful, as, in order to be eclipsing, the inclination must be very close to 90 degrees. Furthermore, the duration of the transits relates to the relative sizes of the stars to the orbit and so mass and radius can be measured simultaneously. Fast-forwarding to today, we have catalogs of hundreds of thousands of binary systems, coming from large all-sky surveys such as OGLE 10 (Graczyk et al.,2011), ASAS 11 (Paczynski et al.,2006) and space missions like Kepler (Borucki et al., 2010), TESS12 (Ricker et al.,2015) and Gaia (Gaia Collaboration et al.,2023). With such catalogs, statistical studies on the period, eccentricity and mass distributions can be undertaken to gain a better understanding of stellar and binary evolution. It would be naive to assume that a star that is part of a binary would evolve in precisely the same way as if it was single. Several key differences are immediately clear: Ι n close binary systems, tidal interactions can directly influence the rotation rate of either star. There is a parallel with the Earth-Moon system, as tidal interactions between the Earth and Moon have made the rotation period of the Moon equal its orbital period around the Earth (so that it always shows the same face to us, we call this process tidal locking). In a single star however, rotation is unaffected by external forces. But more important is the proximity of the companion itself. As stellar evolution drives stars to become bigger, at some point the more massive star of the two, which evolves faster, grows to the point where the gravitational pull of the companion on its outer atmospheric layers supersedes its own. At that moment, a phase of mass transfer will be initiated, where material of the donor star falls toward the accretor star. Constructing detailed models of this form of binary interaction, and the impact on the evolution of the stars, is the subject of many current efforts in binary-star astrophysics. 10Optical Gravitational Lensing Experiment 11All Sky Automated Survey 12Transiting Exoplanet Survey Satellite 40 INTRODUCTION giant. If we then make the approximation that the dense cores of the stars are point masses, and that the total gravitational potential is unaffected by the tenuous outer layers of the stars, we get: 𝑉(𝒙)=−𝐺𝑀1 |𝒙|−𝐺𝑀2 |𝒙−𝒂|,(1.34) where 𝒂 is the vector separating the centers of mass of the stars. From this potential it follows that the centers of mass of the stars move around one another in ellipses, just as Johannes Kepler found for the planets around the Sun. If we make one final simplification, namely that the orbits of the stars are circular, and then move our reference frame into the rotating binary, we arrive at the Roche potential of the binary star: 𝑉(𝒙)=−𝐺𝑀1 |𝒙|−𝐺𝑀2 |𝒙−𝒂|−1 2(𝜴×𝒙)2,(1.35) where 𝜴 is the angular velocity vector pointing along the axis of rotation. Its magnitude, 𝛺2=𝐺(𝑀1+𝑀2) 𝑎3,(1.36) follows from Kepler’s third law. Equation (1.35) is named after nineteenth-century French mathematician Édouard Roche, who studied the restricted three-body problem using this potential. Figure 1.16 shows an impression of the Roche potential in the equatorial plane. We can see the two wells of the point masses that represent the heavy cores of the stars, as well as the outer regions that are dominated by the centrifugal force which appeared when changing into an accelerating reference frame. Without any other force, the slopes in Fig. 1.16 determine how a marble would start to roll if placed in this potential. If we now consider surfaces of constant potential, i.e., equipotentials, in the equatorial plane, we arrive at the contours pictured in Fig. 1.17. Marked there also are the Lagrangian points, L 1 through L 5 . These are special locations where the gradient (i.e., the slope) of the potential is zero. A particle that is exactly at any of the Lagrangian points would thus stay there indefinitely if not disturbed by an outside force. However, points L 1 , L 2 and L 3 are unstable, meaning that the tiniest perturbation from them will amplify, and the particle will leave that point. In contrast, L 4 and L 5 can be stable thanks to the Coriolis force acting on a moving particle (which cannot be represented in a potential graph).14 The existence of the Lagrangian points has important ramifications for the life of stars that have a nearby companion. Generally, a star will grow in radius throughout its life. This is most apparent when a star exhausts its hydrogen in the 14 This is not a trivial calculation. Only for mass ratios bigger than 𝑀1/𝑀2=25+3√69 2≈ 25 are L 4 and L 5 stable. Furthermore, they are never asymptotically stable, meaning a particle will never exactly return to L 4 or L 5 , but only stay in their vicinity. This nice result can be obtained from the restricted three-body problem by solving the linearized equations of motion around the Lagrangian points (Greenspan,2014). BINARY STARS 41 Figure 1.16 : Representation of the Roche potential in the equatorial plane for 𝑞=𝑀2 𝑀1= 0 . 3. The deep gravitational wells are created by the massive cores of the stars, while the outer regions are dominated by the centrifugal force. Figure 1.17 : Roche equipotential contours in the equatorial plane (for 𝑞=𝑀2 𝑀1= 0 . 3), along with the five Lagrangian points marked as red dots. 42 INTRODUCTION core and becomes a giant, where the stellar size changes typically by two orders of magnitude. But, even during the main sequence itself, as the chemical mixture in the core is slowly changing from hydrogen to helium, the star needs to adapt and grow in size to accommodate for the increasing luminosity generated (see Fig. 1.8). The first Lagrangian point, L 1 , is a showstopper however, as, when a star grows so its surface crosses L 1 , mass will start funneling toward the companion. Figure 1.18 shows another view of the Roche potential, this time along the axis connecting the centers of the stars. In the left panel, both stars sit nicely inside their own wells, the so-called Roche lobes (RLs), while in the right panel, the more massive star has grown to fill its RL. We give these binary configurations the names ‘detached’ and ‘semi-detached’, respectively. The size of the RL depends on the separation, 𝑎 , and the mass ratio, 𝑞 , of the binary, and can be computed by performing a numerical integration of the Roche potential of Eq. (1.35) . A widely used analytical fit for the volume-equivalent radius, i.e., the radius of a sphere having the same volume as the teardrop-shaped RL, exists (Eggleton,1983): 𝑅RL,1=0.49𝑞−2/3 0.6𝑞−2/3+ln 1+𝑞−1/3𝑎, (1.37) where we defined 𝑞=𝑀2/𝑀1 . The main feature of this equation is that the RL size increases as the separation (or orbital period) increases, and vice versa. The process of mass spilling over the L 1 point is called mass transfer (MT), or Roche-lobe overflow (RLOF), and material that was once part of the donor star can end up on the accretor star. For main-sequence donor stars, we expect its radius to never significantly increase beyond the RL, because the mass-transfer rate, i.e., the amount of mass transferred per second, increases with the degree of overflow (Kolb & Ritter,1990). For giant donors however, whose outer layers are rarefied and tenuous, the overflow can become quite large even for moderate mass-transfer rates. Mass-transfer phases are categorized depending on the evolutionary status of the donor star, coined by Kippenhahn & Weigert (1967). For main-sequence donors, we use the case A designation. For post-main-sequence donors, we use case B and case C for stars that have exhausted hydrogen and helium in their cores, respectively.15 It is important to consider also the effects of mass transfer on the orbital parameters themselves, and here we focus on the case where mass transfer is fully conservative, i.e., all mass lost from the donor is accreted onto the companion. When mass is transferred from a more massive star to a less massive star, conservation of (orbital) angular momentum says the separation between the stars must shrink. Conversely, mass transfer from a lower-mass star to a higher-mass star will widen the orbit. This means that a mass-transfer phase of a main-sequence donor (case A) can be 15 Although Lauterborn (1970) states Kippenhahn & Weigert (1967) categorized all three cases, to our knowledge, Kippenhahn & Weigert (1967) made no mention of case C. Instead, we find Lauterborn (1970) makes the first mention of case C mass transfer. BINARY STARS 43 Figure 1.18: Configurations in the Roche potential. On the left, both stars reside within their own Roche lobes, while on the right, mass transfer through Roche-lobe overflow is occurring. split into two sub-phases. First, as the more massive star reaches its RL, it will start mass transfer to a less massive companion, meaning the orbit will shrink, resulting in a decreasing RL size. The donor star cannot thermally adjust quick enough to this new RL size, and will overflow it even more, increasing the mass-transfer rate, which exacerbates the problem. This scenario is called fast or thermal timescale mass-transfer, as it takes (roughly) one thermal timescale for the donor to adjust its radius to its new mass (Pols,1994). In most binary systems, the accretor star will be gaining so much mass during fast case A it becomes the more massive star in the system. From that moment on, the binary separation widens as a result of further mass transfer, which increases the RL size of both components. Once the donor has regained thermal equilibrium, and the star can immediately adjust its radius to its new mass, the mass-transfer rates remain rather benign as it is the nuclear expansion of the donor that drives mass transfer. This scenario is called slow case A mass-transfer, and occurs on the nuclear timescale of the donor. In very close binary evolution, both phases occur one after the other (see, e.g., Sect. 4.4.1). The transfer of mass between components in a binary is the main interaction process by which binary-star evolution differs significantly from single-star evolution. Binary interaction has been very successful in explaining various astrophysical phenomena that could not be understood in terms of single-star structure and evolution. One famous such case is the so-called Algol paradox. Algol is both a spectroscopic and eclipsing binary, 16 and from its analysis the masses and radii of the two stars could be determined with good precision. The primary star, Algol A, has 𝑀𝐴= 3 . 17 ± 0 . 21 𝑀⊙ and 𝑅𝐴= 2 . 73 ± 0 . 20 𝑅⊙ , while the secondary, Algol B, has 𝑀𝐵= 0 . 70 ± 0 . 08 𝑀⊙ and 𝑅𝐵= 3 . 48 ± 0 . 28 𝑅⊙ (Baron et al., 2012). If we make the assumption that these two stars were formed at the same time, then conventional, single-star evolution arguments should make the more 16It is actually a triple star system, but we focus on the inner binary here. 44 INTRODUCTION Figure 1.19: Contact configuration of a binary star. Both stars are overflowing their respective RLs, and share common layers above L 1 . The surface layer of the contact binary is a Roche equipotential. massive star evolve faster and become larger before the less massive star. In Algol, this is clearly not the case as 𝑀𝐴> 𝑀𝐵 , but 𝑅𝐵> 𝑅𝐴 , so it seems the less massive star has grown faster than the more massive one! In principle it could be that both stars were not born simultaneously, and later paired up to become the binary star we see today, but that is impossible as a 0 . 7 𝑀⊙ star has not had enough time yet to expand to 3 . 5 𝑅⊙ , the universe is too young. Instead, mass transfer can fairly easily resolve the paradox: In the past, Algol B was more massive than Algol A, and evolved more rapidly, until mass was transferred from Algol B to Algol A, which has flipped the mass ratio. Algol A and B are thus a binary interaction products. Algol A appears younger than its age suggests because for much of its life it was much less massive. Algol B on the other hand appears much older for its mass, for the opposite reason. 1.3 Contact binaries When a star grows to fill its RL, it will start a phase of mass transfer to its companion in a semi-detached configuration. However, it is possible that both stars simultaneously want to overfill their respective RLs. In this case, both stars will grow to have their surface layers on the same equipotential somewhere between L 1 and L 2 , see Fig. 1.19. As a result, the stars enter into a contact configuration, and their shape merits the colloquialism “peanut stars.” In the literature, both the terms overcontact and contact are used for the same configuration. “Overcontact” most likely originates from the portmanteau of “overflow” and “contact,” and some authors distinguish “overcontact” from “contact” from the amount of overflow. There is however no clear, physical boundary separating contact binaries from deep “overcontact” binaries. Therefore, in this thesis, we use the term “contact” for any configuration where both components overflow their RLs. CONTACT BINARIES 45 Figure 1.20: Surface model of VFTS 352, a massive contact binary in the LMC. The color indicates the local effective temperature, which varies across the surface because of the loss of a spherically symmetric gravitational field. Reproduced from Abdul-Masih et al. (2021). 1.3.1 Observed systems The characterization of massive binaries is not an easy task. It requires using several observational techniques to accurately determine the orbital and stellar parameters. In Sect. 1.2.1, we introduced spectroscopic and astrometric constraints, but a combined spectroscopic and photometric analysis can also fully determine a system. For contact binaries, this is the most widely used approach (e.g.,Lorenzo et al.,2014;Abdul-Masih et al.,2021;Kobulnicky et al.,2022). One of the best characterized massive contact binaries is VFTS 352, in the Tarantula star-forming region of the LMC. By simultaneously fitting RV and light curve data, Almeida et al. (2015) found it is a 28 . 6 ± 0 . 3 𝑀⊙+ 28 . 8 ± 0 . 3 𝑀⊙ binary, i.e., a nearly equal-mass system. They also found both stars are overluminous and hotter than what single-star theory would predict. Abdul-Masih et al. (2021) characterized two more contact binaries, and similarly found overluminous components in both of them. A potential source of this discrepancy is that 1D, spherical atmosphere models were fitted to (clearly) non-spherical, highly distorted stars. Abdul-Masih et al. (2020) developed a spectroscopic patch model that more accurately fits the spectra than the 1D models. This method assigns a “local” spectrum to patches across the surface, thus taking into account temperature variations due to tidal distortion (see Fig. 1.20), and integrates the contributions of all patches to come to a composite spectrum. It is quite challenging to confidently say an observed contact binary star is indeed experiencing a contact phase, or whether it is nearly so in a (semi-)detached configuration. Usually, contact binary candidates are identified through their 46 INTRODUCTION photometric light curves obtained from stellar variability surveys like OGLE, ASAS, MACHO, Kepler, Gaia and TESS (Szymanski et al.,2001;Paczynski et al., 2006;Rucinski et al.,2007;Prša et al.,2011a;Ricker et al.,2015;Gaia Collaboration et al.,2016,2023). Such surveys provide hundreds of thousands of (primarily lowmass) candidates. Unfortunately, a light curve of a contact binary looks very much like a “nearly-in-contact-but-still-detached” binary. To confirm the status of the binary, in particular to have an independent constraint on mass ratio, spectroscopic measurements need to be taken, whose analysis has its own challenges alluded to in the previous paragraph. Even with homogeneous datasets, Mahy et al. (2020) could only confirm the contact configuration of VFTS 352. Three other massive contact binary candidates (out of a sample of thirteen targets) have measured parameters whose uncertainty ranges make them consistent with a detached configuration. Despite the observational challenges, over the past decade, the following picture of observed, massive contact binaries emerged. A good fraction of contact binaries are observed with unequal-mass components, meaning their mass ratio, 𝑞=𝑀2/𝑀1 , is significantly away from unity. In Fig. 1.21, we plot the cumulative distribution of observed, massive (near-)contact binaries (Tables E.1 and E.2), and deduce, under the assumption this is a representative sample, that the probability of finding contact binaries with mass ratios more extreme than, say, 𝑞= 0 . 85, is around 50%. Furthermore, from studying the evolution of the orbital period over the past thirty years, Abdul-Masih et al. (2022) determined that the period derivative in six contact binaries is rather low, on the order of the nuclear timescale, meaning their mass-ratio evolution also occurs on the nuclear timescale. It thus seems that contact binaries with unequal-mass components are stable in that configuration. This is an important observational constraint that models of contact binaries would ideally want to reproduce. 1.3.2 Modeling efforts Modeling of contact binaries got well underway in the late 1960s, thanks to growing surveys of variable stars. Many of these variable stars turned out to be partially eclipsing, low-mass contact binaries, called W UMa stars (Eggen,1967; Binnendijk,1970). The surveys revealed many systems contained unequal-mass components. Theoretically, this was not expected, as Kuiper (1941) identified that the mass-radius relationship of single stars in equilibrium, combined with the geometric constraint from the Roche geometry, did not allow for contact binaries with mass ratios different from unity. This discrepancy between theory and observation is called Kuiper’s paradox. To resolve it, Lucy (1967b) first considered the process of energy transfer in the common layers of the contact binary. Also Biermann & Thomas (1972), Vilhu (1973), Flannery (1976), Shu et al. (1976) and Webbink (1977) all worked on advancing the models of W UMa stars. CONTACT BINARIES 47 Figure 1.21: Cumulative probability distribution function of the observed sample of massive (near-)contact binaries (see Tables E.1 and E.2). It follows nearly a 𝑞3 distribution. About a half of the systems are observed at mass ratios more extreme than 𝑞=0.85. Unfortunately, all of these models contained inconsistencies, or were not applicable to the massive contact binaries. For example, Hazlehurst (1993) showed that the model of Shu et al. (1976) violated the second law of thermodynamics, while the energy-transfer model of Lucy (1967b) is not applicable to the radiative envelopes of massive stars. Since the initial efforts to model W UMa stars, little progress has been made at obtaining a complete model of contact binary stars, especially the massive ones. As a result, much of the theoretical work from fifty years ago is ignored in recent modeling. Now that it is clear massive stars undergo contact phases frequently, many studies focus on close binary evolution with specific considerations of the contact phase, e.g. Marchant et al. (2016), Menon et al. (2021), Sen et al. (2022), and Henneco et al. (2024). However, they treat contact phases in a rudimentary way only. Mass transfer is computed consistently, by requiring that the masstransfer rate is adjusted so that the surfaces of both stars lie on a common Roche equipotential. Unfortunately, other effects are ignored in all models and calculations. For example, the spherically symmetric stellar-structure equations, Eqs. (1.5) , or their equivalents for single rotating stars, are solved, so that any effect of tidal deformation cannot be taken into account. Furthermore, Hastings et al. (2020) calculated that internal mixing currents in tidally deformed stars are stronger than the original prescription of Eddington (1925) and Sweet (1950) by 48 INTRODUCTION about a factor of two. Tidally distorted stars could thus live longer and have larger cores as they mix more fresh fuel into their stellar cores. Finally, the process of energy transfer, which was heavily discussed for W UMa stars in the seventies, and taken up again by a few papers in the 2000s (e.g.,Kähler,2004;Yakut & Eggleton, 2005), has never been considered in massive binary evolution. 1.3.3 Evolutionary pathways – Mergers Contact binaries represent a unique, but not-uncommon phase in binary evolution. Recent calculations estimate that at least 40% of all massive stars undergo at least one contact phase throughout their lives (Henneco et al.,2024). Generally speaking, there are two ways a contact configuration might form two with main-sequence components (see Henneco et al.,2024, for a complete view of contact phases in all stars). First, when the accretor in a semi-detached configuration swells as a consequence of it gaining mass, it can overflow its Roche lobe too and engage in contact with its donor (case AR as introduced in Nelson & Eggleton,2001). This happens in vigorous, thermal-timescale mass-transfer phases. Second, in cases of slow mass transfer, the accretor can “overtake” the evolution of the donor (cases AS/AE of Nelson & Eggleton,2001). In this case, the accretor, growing on its nuclear timescale, enters contact with the donor, and the mass-transfer rate will soon reverse. Both these scenarios fall under the “accretor expansion” outcome of Henneco et al. (2024). Around half of all case A binary systems that enter a contact phase are expected to merge (Henneco et al.,2024), and this number rises to essentially all systems if one considers nuclear-timescale contact phases only (Menon et al.,2021). The merger process of two stars is a violent and dynamical occurrence, resulting in matter being flung outward, and a temporary brightening of the system (e.g.,Pejcha et al., 2016). Because the probability of observing a phenomenon is proportional to its duration, stellar mergers are hard to catch, as they happen rather quickly on astronomical timescales. Nevertheless, thanks to long-term monitoring programs and transient surveys, which use dedicated telescopes that scan the whole sky every few days, we have caught several events called Luminous Red Novae (LRNe; see Kasliwal et al.,2017;Pastorello et al.,2019). Such events are dimmer than supernovae, but more luminous than classical novae coming from white dwarfs. The eruptions of V838 Monocerotis and V1309 Scorpii are the canonical examples of LRNe. Both were rather unremarkable stars that quite suddenly brightened by several magnitudes. V838 Mon erupted in January 2002, and was first noticed by Brown et al. (2002) as a “peculiar variable.” In the following years, various models for the event were put forward, like an expanding red giant swallowing its close planets (Retter & Marom,2003), or the merger of a low-mass star with an intermediate-mass one (Soker & Tylenda,2003). Because the distance to this CONTACT BINARIES 49 source was adjusted by analysis of the light echo (i.e., the delayed illumination of the surrounding gas and dust; Tylenda,2004), Tylenda et al. (2005) modeled the progenitor as a massive, 𝑀∼ 8 𝑀⊙ , star. More recently, mass loss from the outer Lagrangian point is also considered as the origin for the eruption of V838 Mon (Pejcha et al.,2016). Apart from its similarity to V1309 Sco, the precise nature of the eruption of V838 Mon has not been uncovered yet. By comparing the various proposed models, Tylenda & Soker (2006) conclude a merger of a 8 𝑀⊙ star with a 0 . 4 𝑀⊙ pre-main-sequence star is the most promising (although it is unsure whether a contact binary formed). V1309 Sco is perhaps even better studied, because there is well-resolved OGLE data of the progenitor, which was undoubtedly a contact binary. Its outburst in 2008 was discovered by Nakano et al. (2008), after which Mason et al. (2010) established that it was not a classical nova. Tylenda et al. (2011) analyzed both the outburst and the progenitor data, and conclusively showed that this object was indeed a contact binary that experienced a dynamical merger event. The OGLE data revealed that, in the years preceding the merger, the orbital period was decreasing exponentially, which hints at a runaway, self-strengthening process. During the outburst, the system brightened by about seven magnitudes in around six months, finishing with a sharp, ten-day rise of three more magnitudes. Soon after, St e¸ pie´n (2011) computed evolutionary models for the progenitor, and estimated initial masses of 𝑀1≈ 1 . 2 𝑀⊙ and 𝑀2≈ 0 . 6 𝑀⊙ for the components, respectively. These stars first reverse their mass ratio by case A mass transfer, after which the second phase of mass transfer becomes unstable, resulting in the observed inspiral and merger. Pejcha (2014) focused on the exponential rise of the outburst, and connected it to rapid mass loss from the outer Lagrangian point L 2 . Several more LRNe have been detected and connected to mergers or common-envelope events (Smith et al.,2016; Pastorello et al.,2019;Blagorodnova et al.,2017,2020,2021). In common-envelope evolution, like contact binaries, both components share their outer layers, but a crucial difference is that there is loss of co-rotation (Ivanova et al.,2013;Röpke & De Marco,2023). While a contact binary can be described as a (quasi-)hydrostatic object, a common-envelope phase is always a hydrodynamical affair, where the core of a star, or a compact object, is moving through the shared envelope of the system. Furthermore, in common-envelope events, the Roche potential is no longer a suitable approximation of the geometry of the system. The result of a merger, once it has thermally relaxed, is often a rather peculiar star. If it is part of a stellar cluster, it will appear younger than its fellow cluster members. It will also appear hotter (“bluer”) on the HRD, giving these stars the name blue stragglers, first identified in the globular cluster M3 by Sandage (1953). While there is agreement that stellar mergers contribute to the number of blue stragglers (Hills & Day,1976), the precise amount is unclear (Leonard,1989). Other possibilities are that they are stars with particularly well mixed material (see Sect. 1.1.5;Wheeler,1979), that they have accreted a significant amount of mass 56 INTRODUCTION a black hole as its nuclear fuel is expended. At this point in the evolution, the orbital period is still on the order of hundreds of days, meaning no gravitational wave merger is expected if the second star also became a black hole at this period. We thus need a mechanism that shrinks the orbit by at least two orders of magnitude. As time moves on, the second star starts evolving significantly, causing it to fill its Roche lobe. If the mass ratio is low, 𝑀BH ≪𝑀2 , the second phase of mass transfer is unstable, forming a common envelope (CE). In common-envelope evolution, the black hole spirals in through the envelope of the giant secondary, and experiences drag. It is then postulated that the gravitational energy liberated from the inspiral is converted to unbind the outer layers of the giant. In this process, the orbit of the black hole can shrink by orders of magnitude, enabling a gravitationalwave merger event to happen in a Hubble time (e.g.,Tutukov & Yungelson, 1993;Stevenson et al.,2017;Giacobbo & Mapelli,2018). While we understand the qualitative picture of a common-envelope inspiral, there are still many uncertainties in its quantitative evolution, in particular on the efficiency of energy conversion, and on the endpoint of the common-envelope phase (see, e.g.,Ivanova et al.,2013, for a review). One of the possible outcomes is that the envelope is not ejected, and that the black hole merges with the core of the giant star, preempting any gravitational-wave merger happening from the system (e.g.,Klencki et al.,2021). 3. Stable mass transfer. A sister of the common-envelope channel, this scenario requires that the second phase of mass transfer is stable. Under certain conditions, if mass lost from the donor is not accreted by the companion, the period of the binary can shrink without the mass-transfer rate increasing uncontrollably (van den Heuvel et al.,2017). In essence, the non-conservative mass transfer acts as an energy and an angular momentum sink, with which the binary can also shrink in period by several orders of magnitude. Characterizing exactly when mass transfer is stable is crucial to understand what parameter space of binaries might go through this channel (Marchant et al.,2021;Klencki et al.,2021;Olejak et al.,2021;Gallegos-Garcia et al., 2021). 4. Population III stars As early as 1926 by Jan Hendrik Oort, and later in 1944 by Walter Baade, stars were categorized in two classes depending on their measured velocity. Population I stars are the slowly moving, metalrich stars of the galactic disk, while population II stars were faster moving, metal-poorer stars located in the galactic bulk and globular clusters in the halo. Since globular clusters are very old, population II stars are older than population I stars. Later, the third class, population III, was added for the stars that were formed very soon after the big bang (in the first billion years or so, Abel et al.,2002). Because they are the very first generation, and one requires stars to form heavy elements, population III stars are thought to be essentially metal-free (Bromm et al.,2002). Population III stars should be THIS THESIS 57 more compact for the same mass as population I stars, so they more easily avoid interaction during their lives. Kinugawa et al. (2014) modeled a large set of population III stars, and computed the compact object merger rates from it. They find that, even though the merger rate from population III stars is lower than for population I stars, they are more easily observable because the average masses of the merging black holes are much higher for the population III stars. This conclusion however depends on cosmological parameters and the star-formation-rate history. The community has not yet reached consensus on which of the channels produces the observed gravitational-wave mergers. Most probably, all of them will contribute to a certain degree, but uncertainties in the modeling prevent us from making precise predictions. While the increasing number of gravitational-wave detections can help constrain the physical processes of the progenitors (Mandel & Farmer,2022), efforts to improve the modeling itself, and understanding its associated physics, are the main way to progress in this growing field. 1.5 This thesis The main goal of this thesis is to advance on the theoretical modeling of massive contact binary stars. Given the frequency of contact phases, their pivotal role in stellar mergers and associated exotic phenomena, and occurrence in the chemically homogeneous evolution channel of gravitational-wave progenitors, it is timely to concentrate efforts on the accurate modeling of contact binaries. However, we will first make a brief but important excursion into the observational characterization of binary stars. Most model-independent masses of (very) massive stars come from short-period eclipsing binaries. However, such measurements might be affected by additional physical effects, such as tidal deformation and unconstrained internal mixing. In this sense, young, long-period binaries, which are unaffected by close binary interaction, offer an excellent alternative to calibrate massive-star models. In Chap. 2, we will analyze the long-period, massive binary star 9 Sagittarii (9 Sgr). The central question we want to solve is whether the spectroscopic analysis of Rauw et al. (2012) suffered from biases using line-profile fitting of the highly blended lines of 9 Sgr. We will use relative astrometry and spectral disentangling to determine all the orbital parameters and in so doing calculate the dynamical masses of both stars. Afterwards, we will analyze the atmospheres of the stars, deriving their temperature and abundances of elements like carbon, nitrogen and oxygen, and compare the stellar parameters with evolutionary models. 58 INTRODUCTION Figure 1.25: Pathway of a binary black hole merger through the common-envelope and stable mass-transfer channels. Two phases of interaction are required to produce merging black holes in this scenario. First a case A or B mass-transfer phase occurs, later followed by either a common-envelope phase or a stable masstransfer phase. Representative periods are indicated at various phases in the evolution. Acronyms are explained in the preface.Adapted from Marchant et al. (2016). THIS THESIS 59 In Chap. 3, we develop the methodology to accurately represent contact binaries in a one-dimensional stellar-evolution code. We will use the method of Kippenhahn & Thomas (1970) and Endal & Sofia (1976) to represent tidally distorted stars under the shellular approximation, and adapt the equations of stellar structure accordingly. We will then quantify how the surface properties of a stellar model change when going from a spherically symmetric geometry to a rotating or tidally deformed geometry. Chapter 4deals with the process of energy transfer in contact binaries. We will give a brief overview of the theory of energy transfer, after which we implement an energy-transfer scheme into our stellar-evolution models. The effects on the shared layers of the stars is carefully studied, and differences in the evolutionary pathways of a binary model with and without energy transfer are identified. In Chap. 5, we will compute a large grid of binary-evolution models, where, in one grid, we include the process of energy transfer, while in another, we exclude it. In this way, we will get a measure of the effects of energy transfer across the entire population of massive binary stars. In particular, we aim to answer whether including energy transfer in theoretical modeling can resolve the discrepancy on the mass-ratio distribution introduced in Sect. 1.3.2. Finally, Chap. 6concludes this thesis with a summary and outlook. Chapter 2 The dynamical mass of 9 Sagittarii This chapter is mainly based on: Resolving the dynamical mass tension of the massive binary 9 Sgr M. Fabry, C. Hawcroft, A. J. Frost, L. Mahy, P. Marchant, J-B. Le Bouquin, H. Sana ASTRONOMY & ASTROPHYSICS, 651, A119 (2021) Authors contributions: M. Fabry did the majority of the work appearing in this chapter. M. Fabry, L. Mahy and H. Sana jointly developed the methodology of the spectral-disentangling approach. J-B. Le Bouquin and A. J. Frost reduced the astrometric data, while C. Hawcroft performed the atmosphere analyses. P. Marchant provided critical insights on the evolutionary models. The text was written by M. Fabry, and subsequently improved on by all co-authors. With respect to the original publication, some of the text was altered to fit the chapter format of this thesis. Original abstract: Context. Direct dynamical mass measurements of stars with masses above 30 𝑀⊙ are rare. This is the result of the low yield of the upper initial-mass function and the limited number of such systems in eclipsing binaries. Long-period double-lined spectroscopic binaries that are also resolved astrometrically offer an alternative to eclipsing binaries for obtaining absolute masses of stellar objects. 9 Sgr (HD 164794) is such long-period high-mass binary. Unfortunately, there exists large 61 62 THE DYNAMICAL MASS OF 9 SAGITTARII tension between its total dynamical mass inferred spectroscopically from radialvelocity measurements and that from astrometric data. Aims. Our goal is to resolve the mass tension of 9 Sgr that exists in literature, to characterize the fundamental parameters and surface abundances of both stars as well as to determine the evolutionary status of the binary system, henceforth providing a reference calibration point to confront evolutionary models at high masses. Methods. We obtain the astrometric orbit from existing and new multi-epoch VLTI/PIONIER and VLTI/GRAVITY interferometric measurements. Using archival and new spectroscopy, we perform a grid-based spectral disentangling search to constrain the semi-amplitudes of the radial-velocity curves. We compute atmospheric parameters and surface abundances by adjusting FASTWIND atmosphere models and we compare our results with evolutionary tracks computed with the Bonn Evolutionary Code (BEC). Results. Grid spectral disentangling of 9 Sgr supports the presence of a 53 𝑀⊙ primary and a 39 𝑀⊙ secondary, in excellent agreement with their observed spectral types. In combination with the size of the apparent orbit, this puts 9 Sgr at a distance of 1 . 31 ± 0 . 06 kpc . Our best-fit models reveal a large mass discrepancy between the dynamical and spectroscopic masses, which we attribute to artifacts from repeated spectral normalization before and after the disentangling process. Comparison with BEC evolutionary tracks shows the components of 9 Sgr are most likely coeval with an age of roughly 1Myr. Conclusions. Our analysis clears up the contradiction between mass and orbital inclination estimates reported in previous studies. We detect the presence of significant CNO-processed material at the surface of the primary, suggesting enhanced internal mixing compared to currently implemented levels in the BEC models. The present measurements provide a high-quality high-mass anchor to validate stellar-evolution models, and to test the efficiency of internal mixing processes. 2.1 Introduction Massive stars drive the chemical enrichment of heavy elements and inject large amounts of kinetic energy into their neighborhoods through their strong, linedriven winds and final explosion as supernovae and gamma-ray bursts. Obtaining accurate mass measurements of stars in the upper part of the HRD has been a challenge, however, as evolutionary models of high-mass stars are riddled with physical uncertainties. Furthermore, spectroscopic masses, obtained through INTRODUCTION 63 atmospheric model fitting, are intrinsically inaccurate. Therefore, direct mass measurements that are independent of atmosphere or evolutionary models offer valuable constraints to gauge the quality of the models. In Sect. 1.2.1, we introduced how model-independent masses can be found through Kepler’s laws, which provide so-called dynamical masses. Unfortunately, no single observational technique can fully characterize the orbit and dynamical masses of the two components. Either a double-lined spectroscopic binary (SB2) has to be eclipsing, or, astrometric data, either absolute or relative, must be available. In both cases, multi-epoch observations are required. Traditionally, eclipsing SB2s are considered to provide the best constraints on the orbital parameters. Yet, they are rare and uncertainties about the effects of tidal deformation, mutual illumination and/or binary interaction may pollute the obtained results, making it challenging to confront these objects to single-star models. This is particularly the case in the realm of massive stars. In this context, astrometric binaries are a valuable alternative to eclipsing SB2s. Recent advances in optical long-baseline interferometry have identified a number of such systems (e.g.,Sana et al.,2013b,2014;Mayer et al.,2014;Maíz Apellániz et al.,2017;Mahy et al.,2018;Lanthermann et al.,2023), which offer new opportunities to obtain dynamical mass constraints of stars in the high-mass regime. 9 Sgr (HD164794) is such a long-period, astrometric SB2 system in the Lagoon Nebula. It was first studied by Abbott et al. (1984) in the context of its variable synchrotron emission, which is interpreted to be a result of wind-wind collisions of binaries (Pittard & Dougherty,2006), hinting that the then presumed single star 9 Sgr was in fact a binary. Subsequent studies by Rauw et al. (2002a,b) and Nazé et al. (2008) confirmed the presence of elevated X-ray emission, a typical indication of colliding winds in O + O binaries (Rauw et al.,2002c;Sana et al.,2004,2006;Rauw & Nazé,2016). Rauw et al. (2012) first confirmed the long-period binary nature of 9 Sgr through radial velocity (RV) measurements and classified its components as O3.5V((f+)) and O5.5V((f)). Rauw et al. (2016) studied the periastron passage of 9 Sgr and reported a maximum in the X-ray emission coming from shocked gas in the interaction zone of the stellar winds, as expected from a wind-wind collision in a wide binary where the shocked material cools adiabatically (Stevens et al., 1992). The system was resolved for the first time in 2009 using the Astronomical MultiBEam combineR (AMBER) and in 2013 with the Precision Integrated-Optics Nearinfrared Imaging ExpeRiment (PIONIER) by Sana et al. (2014). Le Bouquin et al. (2017) constrained the astrometric orbit using multi-epoch AMBER and PIONIER interferometric measurements, uncovering a discrepancy with the spectroscopic analyses of Rauw et al. (2012,2016). While the interferometric measurements of Le Bouquin et al. (2017) firmly excluded inclinations below 80 degrees, the 64 THE DYNAMICAL MASS OF 9 SAGITTARII RV curve semi-amplitudes of Rauw et al. (2012,2016) resulted in an estimated inclination of about 50 degrees if the stars were to have masses representative of their spectral types. Another long standing issue is the distance to 9 Sgr, and whether or not it is a member of the young open cluster NGC 6530. Prisinzano et al. (2005) and Kharchenko et al. (2005) measured distances to the cluster of around 1 . 25 kpc , while earlier measurements indicated a distance close to 1 . 8 kpc (van den Ancker et al.,1997;Sung et al.,2000) as a result of differences in the adopted reddening laws. The combination of the astrometric and spectroscopic measurements provide a direct constraint on the distance, allowing us to confirm the current Gaia eDR3 measurement of 1 . 21 kpc . Gaia does not suffer from systematics owing to an assumed reddening law, and thus provides a powerful constraint, but its measurement can still be impacted by the multiplicity of the system. In this work, we aim to resolve the existing discrepant results that cast doubt on the current mass estimates, evolutionary status, and cluster membership of 9 Sgr. To do so, we leverage the accuracy of relative astrometry and we perform a grid spectral disentangling analysis on spectroscopic data to fully constrain the orbit of 9 Sgr. We further use the atmospheric properties of the stars in the system to derive its evolutionary status and test stellar-evolution models at high masses. This chapter is organized as follows. Section 2.2 covers the observational data that were used. We present and discuss the orbital analysis and spectral disentangling in Sect. 2.3, the atmosphere modeling of both components in Sect. 2.4, and the evolutionary status of 9 Sgr is discussed in Sect. 2.5. Section 2.6 presents our conclusions and final remarks. 2.2 Observations We combine archival and new optical spectroscopy with near-infrared (NIR) interferometry of 9 Sgr. Most of these measurements were part of long-term monitoring programs. 2.2.1 Optical spectra Archival data consist of 57 spectra from the High Efficiency and Resolution Mercator Echelle Spectrograph (HERMES, Raskin et al. 2011, used in Rauw et al. 2016), 20 Fiber-fed Extended Range Optical Spectrograph (FEROS) spectra (Kaufer et al.,1999) and 49 Ultraviolet and Visual Echelle Spectrograph (UVES) spectra (Dekker et al.,2000) in the blue and red arms (used in Rauw et al.,2012). Observations that are not previously analyzed are listed in Table A.1 and consist OBSERVATIONS 65 of five additional HERMES spectra and one spectrum from the High Accuracy Radial velocity Planet Searcher (HARPS) spectrograph (Mayor et al.,2003). The FEROS spectra cover the spectral domain between 3700 and 9000 Å and have a resolving power of 𝑅≈ 48000. The UVES spectra cover the wavelength range 𝜆𝜆 = 3500 – 5000 Å with its blue arm and 𝜆𝜆 = 5000 – 7000 Å with its red arm, and each have a resolving power of 𝑅≈ 40000. The HERMES and HARPS spectra both have 𝑅≈ 85000 with a coverage of 𝜆𝜆 = 3800 – 9000 Å and 𝜆𝜆 = 3750 – 6900 Å , respectively. The FEROS, UVES, and HARPS spectra were obtained through the ESO archive science portal and were pre-reduced with their respective pipelines. The HERMES spectra are reduced using the HERMES Data Reduction Software (DRS) pipeline. Finally, all spectra were normalized over their whole spectral domain by fitting a cubic spline function through selected continuum regions. 2.2.2 Near-infrared Interferometry We used the previously published AMBER and PIONIER dataset obtained from Jun 2009 to Aug 2016 (Le Bouquin et al.,2017), along with two new PIONIER observations obtained in May and August 2017. These new data show for the first time the system turning back on the apparent ellipse and provide an almost complete coverage of the nine-year orbit. They were obtained with PIONIER (Le Bouquin et al.,2011) at the Very Large Telescope Interferometer (VLTI) using the four 1.8 meter Auxiliary Telescopes (ATs) in configurations B2-K0-D0-I3 and A0-G1-J2-J3, offering a maximal projected baseline of 120 and 130 meter respectively. The PIONIER data were reduced and analyzed, as described in Le Bouquin et al. (2017), using the pndrs package. 1 Each observation produces six visibilities and four closure phases, delivering relative astrometry with submilliarcsecond precision and an H-band flux ratio of 𝑓H= 0 . 62 ± 0 . 02. The full journal of interferometric observations is given in Table A.2 along with the measured astrometric properties of the system. Additionally, three VLTI/GRAVITY (GRAVITY Collaboration et al.,2017) measurements were taken, of which two were obtained in June 2016 and the other in September 2016. These observations were part of the science verification (SV) program and used the ATs in configurations A0-G1-D0-C1 and A0-G1-J2-K0. As GRAVITY is a spectro-interferometer, it provides six visibilities and four closure phases for each wavelength bin in the 2 . 0 – 2 . 4 µm NIR K band with a spectral resolving power of 𝑅≈ 4000. These SV data were reduced with the standard GRAVITY pipeline (Lapeyrere et al. 2014, version 1.0.11) and fitted parametrically to a binary model with PMOIRED. 2 The uncertainties on the relative astrometry were estimated by adding a bootstrapped error and a systematic error of 0 . 1 mas in quadrature. In the bootstrapping procedure, data were drawn randomly to create 1http://www.jmmc.fr/pndrs 2https://github.com/amerand/PMOIRED 72 THE DYNAMICAL MASS OF 9 SAGITTARII Figure 2.3: Two-dimensional histogram of the results of the MC sampling in (𝐾1, 𝐾2) space. The color indicates the number of samples 𝑁 (out of a total of 1500) at the corresponding (𝐾1, 𝐾2) pairs. The result from the original data is encircled in red. 2.3.3 Distance Because the semi-amplitudes 𝐾1 and 𝐾2 set the absolute scale of the orbit, we can infer the distance to 9 Sgr by measuring against the apparent orbit found through the interferometry (Fig. 2.1), via the total mass of the system and Kepler’s third law, as follows: 𝑑=1 𝑎app 3 √︄𝐺(𝑀1+𝑀2)𝑃2 4𝜋2,(2.6) where 𝑎app is the semimajor axis of the apparent orbit in angular units. The values (𝐾1, 𝐾2)=( 36 , 49 )km s −1 correspond to a total dynamical mass of 92 𝑀⊙ , which results in a distance of 𝑑=1310 ±60pc,(2.7) where we propagated the MC errors. This value lies within the Gaia eDR3 geometric distance reported by Bailer-Jones et al. (2021) of 1218 +100 −108 pc . With these results, it seems likely that 9 Sgr is a member of NGC 6530, confirming the measurements of Prisinzano et al. (2005) and especially Kharchenko et al. (2005), who quoted a distance of 𝑑= 1322 pc . ORBITAL ANALYSIS 73 Figure 2.4: Disentangled, renormalized spectra of both components using 𝐾1= 36 km s −1 and 𝐾2= 49 km s −1 . The spectrum of the primary is shifted up by 0.4 units with respect to the secondary for clarity. In H 𝛼 , the core is contaminated by nebular emission present in the FEROS and HERMES spectra. The model spectra are those obtained with the best-fit atmospheric parameters from the FASTWIND analysis (Sect. 2.4). Conversely, we can firmly exclude distances to 9 Sgr of over 1500 pc because that would require a total mass of over 130 𝑀⊙ , which our grid-disentangling results do not support. Similarly, the distance of 1780 ± 80 pc of Sung et al. (2000), adopted by Rauw et al. (2012), is incompatible with the interferometric and Gaia measurements of the apparent size of the orbit and the derived component masses. 74 THE DYNAMICAL MASS OF 9 SAGITTARII 2.4 Atmosphere modeling 2.4.1 Setup Using the disentangled spectra, we adjust theoretical line profiles computed with the FASTWIND non-local thermal equilibrium atmosphere code suitable for the expanding atmosphere of O-type stars (Santolaya-Rey et al.,1997;Puls et al.,2005;Carneiro et al.,2016;Puls,2017;Sundqvist & Puls,2018). To reduce the dimensionality of the parameter space, the rotational and macroturbulent velocities are estimated using the iacob-broad tool (Simón-Díaz & Herrero, 2014). Analyzing (with the goodness-of-fit method) the OIII 𝜆 5592 and NV 𝜆 4603 lines for the primary and OIII 𝜆 5592, He I𝜆 4713 and He I𝜆 5876 for the secondary, we find a projected rotational velocity 𝑣sin𝑖= 102 +8 −12 km s −1 and a macroturbulent velocity 𝑣mac = 77 +23 −20 km s −1 for the primary star, while we take 𝑣sin𝑖= 67 +6 −13 km s −1 and 𝑣mac = 48 +21 −14 km s −1 for the secondary. Similar values are obtained when using the Fourier transform method. The stellar atmosphere models are then iterated using a genetic algorithm (Charbonneau,1995;Mokiem et al.,2005) within a predefined parameter space to optimize a 𝜒2 fitness metric until a convergence to the best fit with the spectrum is reached. The version of the genetic algorithm used is detailed in Abdul-Masih et al. (2019). We set the 𝛽 exponent of the wind-acceleration law to 0.85, as appropriate for main-sequence stars (Muijres et al.,2012). The microturbulence velocity is fixed to 𝑣mic = 10 km s −1 in the computation by FASTWIND; in the formal integral it is selected on the criteria 𝑣mic =max( 10 km s −1, 0 . 1 𝑣wind) . We include optically thin wind clumping, with a near-constant clumping factor 𝑓cl throughout the wind. Lastly, we opted to clip the core of the H 𝛼 line to avoid fitting the nebular emission remnant (visible in the bottom-rightmost panel of Fig. 2.4). The full list of fitted spectral lines is shown in Table A.3. To provide an absolute magnitude anchor point for the atmospheric model, from which the component luminosities will be computed, we adopt the photometric data from the Two Micron All-Sky Survey (2MASS, Skrutskie et al. 2006), giving an apparent magnitude in the NIR K S band of 𝑚KS= 5 . 731 ± 0 . 024 for the total system. Using the interstellar absorption coefficient 𝐴𝑉= 1 . 338 ± 0 . 021 from Maíz Apellániz & Barbá (2018), the color correction 𝐴KS/𝐴𝑉= 0 . 116 from Fitzpatrick (1999) and a distance of 1 . 31 ± 0 . 06 kpc calculated in Sect. 2.3.3, the absolute magnitude in the K S band comes to 𝑀KS=− 5 . 01 ± 0 . 10. We correct for the measured PIONIER flux ratio 𝑓= 0 . 62 between the components, where we assume that it remains unchanged between the H and K S band (as expected for such hot objects). This yields absolute component magnitudes of 𝑀KS=− 4 . 49 ± 0 . 10 and − 3 . 96 ± 0 . 10 for the primary and secondary, respectively. We note that these values correspond to fainter stars than their derived spectral types suggest. Synthetic photometry of ATMOSPHERE MODELING 75 Martins & Plez (2006) give 𝑀K=− 4 . 98 and 𝑀K=− 4 . 39 for the O3V primary and O5V secondary, respectively, which suggests the stars are slightly more compact. Finally, using the bolometric correction in the K-band from Martins & Plez (2006), 𝐵𝐶K=28.80 −7.24log(𝑇eff), and the magnitude-luminosity relation log(𝐿/𝐿⊙)=−0.4(𝑀K+𝐵𝐶K−4.75),(2.8) we compute the component luminosities from the K-band absolute magnitudes and the modeled effective temperatures. 2.4.2 Results and discussion We list the spectroscopic parameters of the resulting best-fit atmospheric models and the resulting inferred parameters in the leftmost column of Table 2.4. The corresponding theoretical spectra are plotted in Fig. 2.4 with dashed lines. We note the obvious nebular contamination of H 𝛼 , as well as the general trend that the disentangled spectra are slightly shallower than the model spectra in the deep and broad lines (like H 𝛿 , H 𝛾 and He II 𝜆 4686), suggesting issues in the normalization of these broad lines. The best-fit parameters depend on which line features were considered in the fit and with what weights. For example, giving more weight to He I lines would result in lower inferred effective temperatures and vice versa. Correspondingly, the inferred surface gravities would be lower for lower 𝑇eff and vice versa. Therefore, conservative errors of 1 kK and 0 . 2 dex on 𝑇eff and log𝑔 , respectively, are adopted. Furthermore, the determination of the quality of the CNO abundance measurements is challenging. The carbon abundance for the secondary, for example, is fitted to 𝜀C=log(𝑁C/𝑁H) + 12 = 9 . 12; this is an unusually high measurement that has (to our knowledge) never before been observed. We note that this measurement is driven by the CIII 𝜆 5696 line, which is in emission. At 𝑇eff = 42 kK ,FASTWIND can only reconcile this line in emission by boosting the carbon abundance. Keeping the issues presented by this line in mind, however (see Martins & Hillier,2012), we adopt the minimum of a formal 0 . 1 dex error and the statistical error of the grid of models. Since this formal error is somewhat arbitrary, even these uncertainties should be interpreted with great care. Therefore we can only argue for qualitative enrichment of nitrogen in the primary and enrichment of carbon in the secondary. For added justification, we show in Figs. A.5 and A.6 the comparison of the model spectra, their error ranges along with spectra using the Brott et al. (2011) CNO baseline abundances in various diagnostic lines of the CNO elements. The best-fit log𝑔 values then provide the spectroscopic masses of the stars, which are found to be 32 ± 16 𝑀⊙ and 19 ± 10 𝑀⊙ for the primary and secondary, respectively. These masses are significantly lower than their dynamical 76 THE DYNAMICAL MASS OF 9 SAGITTARII counterparts, albeit with large error bars, and are not representative of dwarf stars of that luminosity. The main reason for this discrepancy is the low inferred log𝑔 , which should be raised by about 0 . 25 dex for both stars, that is, slightly beyond the adopted uncertainty, to match the dynamical masses. The mass discrepancy problem (Herrero et al.,1992) is still an open issue in massive-star spectroscopy, and while more recent studies (e.g.,Mahy et al.,2020) show that for stars above ∼ 35 𝑀⊙ , the discrepancy largely disappears, in this analysis, it is still present. The repeated normalization of the spectra before and after disentangling could be the root cause of this fact, as the reconstruction plots in Fig. A.3 and A.4 and the comparison to the model spectra (Fig. 2.4) hint towards. 2.5 Evolutionary modeling We compare our previous results with the Milky Way evolutionary tracks of Brott et al. (2011), using the Bayesian search tool BONNSAI (Schneider et al.,2014). 5 The BONNSAI tool allows us to search the rotating single-star evolution tracks of Brott et al. (2011) for the highest-likelihood stellar model that corresponds to measured quantities. We input the observed log 𝐿, 𝑇eff,𝑌He , and 𝑣sin𝑖 from the left column of Table 2.4 and request the highest likelihood models of both stars in the grid. To avoid biasing the Bayesian search, we refrain from using the log𝑔 due to the uncertainties posed by the mass discrepancy. In a first search, we do not input the CNO abundances obtained from FASTWIND. The parameters from the highest likelihood models replicated from our spectroscopic and photometric observables of this search are given in the middle column of Table 2.4. The comparison with evolutionary tracks point toward relatively compact and coeval stars with an age of about 1 Myr . We find that the evolutionary masses are within error of the dynamical masses, which provides a further indication that the spectroscopic mass likely suffers from systematic errors. Additionally, the evolutionary models favor lower CNO surface abundances than are spectroscopically inferred, especially nitrogen in the primary and carbon in the secondary. The rather modest rotational velocities and young ages do not allow for rotational mixing to modify the surface composition; the CNO abundances returned by BONNSAI correspond to the baseline value of the Brott et al. (2011) models with very small uncertainties. From the atmosphere models in Table 2.4, it is clear that the CNO composition between the primary and the secondary is different. This is not reflected in the evolutionary tracks because both models prefer the baseline values of the Brott et al. (2011) tracks, namely 𝜀C= 8 . 13 , 𝜀N= 7 . 64 and 𝜀O= 8 . 55 (middle column of Table 2.4). Furthermore, if we assume the abundances of the secondary are 5The BONNSAI web-service is available at www.astro.uni-bonn.de/stars/bonnsai. EVOLUTIONARY MODELING 77 Table 2.4: Parameters of the best-fit, genetically evolved FASTWIND atmospheric model (described in Sect. 2.4), along with replicated observables from the Brott et al. (2011) models using BONNSAI (Sect. 2.5). The errors correspond to the 1 𝜎 confidence level. Empty entries are indeterminable for that parameter. FASTWIND BONNSAI Primary Secondary Primary Secondary Primary, Scaled CNO 𝑇eff/kK 46.0±1.0 42.0±1.0 45.9+0.6 −0.941.9±0.9 46.0+0.6 −1.0 log(𝑔/[cgs]) 3.87 ±0.20 3.87 ±0.20 4.11 ±0.05 4.12+0.06 −0.07 4.10+0.02 −0.06 log ¤ 𝑀 𝑀⊙yr−1−6.6±0.2−6.6±0.2 ... ... ... 𝑓cl 29 ±5 22 ±3 ... ... ... 𝑣sin𝑖/kms−1a 102+8 −12 67+6 −13 ... ... ... 𝑣rot/kms−1... ... 110+59 −26 70+8 −15 330+26 −30 𝑌He 0.25 ±0.04 0.24 ±0.03 0.26b0.26b0.28+0.08 −0.02 𝜀C8.17+0.60 −0.55 9.12 ±0.10*8.14+0.01 −0.03 8.13b7.12+0.55 −0.05 𝜀N8.45+0.10 −0.29 7.42 ±0.10 7.63+0.09 −0.01 7.64b8.72+0.10 −0.27 𝜀O8.63+0.10 −0.70 8.64+0.10 −0.13 8.55+0.01 −0.02 8.55b8.55+0.01 −0.61 log(𝐿/𝐿⊙)5.68 ±0.08 5.35 ±0.08 5.64+0.07 −0.06 5.33+0.08 −0.06 5.67+0.06 −0.07 𝑅/𝑅⊙10.8±1.0 8.9±1.2 10.45+0.88 −0.59 8.73+0.75 −0.67 10.73+0.79 −0.61 𝑀spec/𝑀⊙32.1±16.0 18.9±10.1 ... ... ... 𝑀evol/𝑀⊙... ... 53.4+3.2 −3.337.0+2.0 −2.353.8±4.7 Age (Myr)... ... 0.52+0.32 −0.33 1.00+0.48 −0.58 1.00+0.80 −0.41 a Determined using IACOB-BROAD, not FASTWIND. b Very small error, see Sect. 2.5. * Highly uncertain measurement, this formal error is likely not representative, see Sect. 2.4.2. baseline for 9 Sgr, this source has a different CNO baseline than the Brott et al. (2011) tracks. It is hard to justify observationally that the observed abundances of the secondary are baseline for 9 Sgr or the cluster NGC 6530. But we expect the least massive star with the lowest rotational velocity to be least contaminated by surface enrichment from a theoretical standpoint. We thus test if BONNSAI finds different models for the primary if we scale down the observed CNO abundances to the Brott et al. (2011) baseline. For the abundances of the primary to maintain the same fractional difference versus the secondary, this amounts to calculating 𝜀′ X=𝜀X,base,Brott −𝜀X,base,obs +𝜀X,prim,obs . Keeping the doubtful C abundance measurement of the secondary in mind (Sect. 2.4.2), we refrain from scaling C and compute 𝜀′ N= 8 . 67 +0.14 −0.31 and 𝜀′ O= 8 . 54 +0.14 −0.71 . Using then again log 𝐿, 𝑇eff,𝑌He , and 𝑣sin𝑖from the FASTWIND models, along with these new N and O abundances as input for BONNSAI, we obtain other highest likelihood evolutionary parameters; these are listed in the rightmost column of Table 2.4. These results point to a different scenario. Here the rotational velocity is significantly higher, allowing significant rotational mixing to occur. The primary age has increased to match that of the secondary as well, while the evolutionary mass, log𝑔and 𝑇eff are only slightly changed when comparing to the non-scaled results (middle column of 78 THE DYNAMICAL MASS OF 9 SAGITTARII Table 2.4). The major implication of this is that either the rotational axis and the normal to the orbital plane are heavily inclined, up to an estimated 68 degrees to explain the observed projected rotational velocity, or the effect of rotational mixing is underestimated in the stellar-evolution models. We cannot exclude either that a mixing mechanism weakly dependent on rotation is operating on stars in this mass regime. Distinguishing between these scenarios however requires greater confidence in the quality of the CNO abundance measurements. We summarize our results of Sects. 2.4 and 2.5 in an HRD and a Kiel diagram that is overplotted by several of the evolutionary tracks and isochrones of Brott et al. (2011). We note that while the location of the models on the HR diagram matches well, there is a mismatch of the spectroscopic mass inferred from the FASTWIND models and the evolutionary masses. In the Kiel diagram, there is a poorer match as expected from the mass discrepancy discussed in Sect. 2.4.2. 2.6 Conclusions In this chapter, we have obtained disentangled spectra of 9 Sgr using a combination of high angular resolution astrometry and spectral grid disentangling with the fd3 code. The astrometric measurements solidify the long period of 8.9 yr and have a near edge-on inclination of 86.5 degrees. Our results confirm the presence of an O3V+O5V massive binary, which has inferred dynamical masses of about 53 and 39 𝑀⊙ , making 9 Sgr one of the most massive galactic O+O binaries ever resolved. Furthermore, to our knowledge, this is only the second instance of a dynamical mass estimate of a galactic O3V star (the other from Mahy et al. 2018). By re-deriving the semi-amplitudes of the RV curves through a spectral disentangling analysis, we clear up the contradictory results between the previous RV measurements of Rauw et al. (2012) and the high inclination of the interferometric orbit by Le Bouquin et al. (2017). Furthermore, the results show 9 Sgr is a member of the young open cluster NGC 6530. The combined dynamical, atmospheric, and evolutionary modeling shows 9 Sgr contains massive stars of roughly 53 𝑀⊙ and 37 𝑀⊙ for the primary and secondary, respectively. 9 Sgr is a unique system in the far top left corner of the HRD, and therefore provides an equally unique opportunity to use its stellar and systemic parameters to compare with massive star evolutionary models as well as binary formation scenarios. CONCLUSIONS 79 Figure 2.5: Top: HRD showing the location of the best-fit FASTWIND models (FW) and highest-likelihood evolutionary models from Brott et al. (2011) (BONNSAI). Overplotted are evolutionary tracks for different masses and initial rotation velocities (grey lines), along with isochrones for the 100 km s −1 initial rotation velocity models (blue lines). Bottom: Kiel diagram with an equivalent legend as the HRD above. Using the scaled CNO abundances does not appreciably move the model of the primary in either diagram. Chapter 3 Modeling tidal distortion in 1D stellar structure “A beginning is the time for taking the most delicate care that the balances are correct.” — Princess Irulan This chapter is mainly based on: Modeling overcontact binaries, I. The effect of tidal distortion M. Fabry, P. Marchant and H. Sana ASTRONOMY & ASTROPHYSICS, 661, A123 (2022) Authors contributions: M. Fabry did the majority of the work appearing in this chapter. P. Marchant assisted in developing the methodology, while H. Sana provided critical feedback throughout the project. The text was written by M. Fabry and improved on by all co-authors. With respect to the original publication, some of the text was altered to fit the chapter format of this thesis. Original Abstract: Context. In the realm of massive stars, strong binary interactions are commonplace. One extreme case is that of contact systems, which are expected to be part of the evolution of all stars evolving towards a merger and hypothesized as playing a 81 88 MODELING TIDAL DISTORTION IN 1D STELLAR STRUCTURE Figure 3.1: Illustration of the Roche geometry in the equatorial plane. It shows the L 1, L 2 , and L 3 equipotentials (gray), as well as the splitting surfaces (blue dashed). The splitting surfaces have their normals always perpendicular to the potential gradient. 3.2.4 Contact shells In the context of 1D evolution models, we have to distinguish the two stellar components within the geometry of a Roche binary. For equipotential shells lying within their respective RLs, this is trivial, as such layers are physically separated, but shells in contact are shared between the two stars. We therefore construct three “splitting surfaces,” one through each of the colinear Lagrangian points L 1, L 2 , and L 3 , separating the Roche geometry in three main parts. This will ensure distinction between the two components, as well as either one from regions beyond the outer Lagrangian points L 2 and L 3 , from which outflows are expected if any component overflows them. We construct them in such a way that the condition: ∇𝜳·d𝑺=0 (3.18) holds on all points of these splitting surfaces (see Fig. 3.1, in blue dashed lines). In other words, the surface normal is always perpendicular to the gradient of the potential of Eq. (3.17) . We note that by this construction of the splitting surfaces, the scalar product d 𝒏· d 𝑺 is zero. This ensures that Eq. (3.4) still holds, and if the contact layers are shellular, there is no net flux of radiative energy between the components following Eq. (3.11) . However, since the condition of shellularity cannot be held exactly, there will exist a flow 𝐿𝛹 ,trans =∫𝑑𝛹 F· d 𝑺split across equipotentials, representing energy transfer from one component to the other. In this work, we assume this contribution is zero and defer the modeling of energy METHODS 89 transfer to later work. We note also that in the case 𝑞≠ 1, the L 1 splitting surface is not a vertical plane through L 1 , which differentiates our results of equipotential volumes and surface areas from those of Mochnacki (1984) and Marchant et al. (2021). 3.2.5 Numerical calculations In order to use the modified stellar-structure equations in a stellar evolution code, we require for each shell in the geometry the quantities 𝑉𝛹 , 𝑆𝛹 , ⟨𝑔⟩ and 𝑔−1 which specify 𝑟𝛹, 𝑓𝑃 and 𝑓𝑇 of Eqs. (3.1) , (3.9) , and (3.14) . Additionally, the specific moment of inertia, 𝑖rot , is needed to calculate the rotational velocity of the equipotential shell from the angular momentum: 𝛺𝛹=𝑗rot/𝑖rot . We approach this problem using numerical methods. First of all we note that the problem is symmetric under 𝑦→ −𝑦 and 𝑧→ −𝑧 , so we need only explicitly do calculations in the quadrant of positive 𝑦and 𝑧. Given a set of mass ratios, 𝑞 , and target equipotential values, 𝛹0 , we calculate, in a spherical coordinate system on rays of constant (𝜃, 𝜙) , the points 𝑟=𝑟(𝜃, 𝜙) that lie on the corresponding equipotential surface, such that: 𝛹0−𝛹(𝑟, 𝜃, 𝜙)=0.(3.19) When solving Eq. (3.19) , care has to be taken in order to avoid roots beyond the splitting surfaces (Sect. 3.2.4), as these belong to a part the geometry not associated to the considered component star, and would instead belong to the other component or to the region beyond the outer Lagrangian points. We therefore restrict the solver to values 𝑟 smaller than the value 𝑟split(𝜃, 𝜙) on the splitting surface that the ray (𝜃, 𝜙) crosses. If it finds no roots in the interval ( 0 , 𝑟split) , 𝑟split is recorded instead as the splitting surfaces mark the maximal extent of the components. 90 MODELING TIDAL DISTORTION IN 1D STELLAR STRUCTURE Next, for these points 𝑟(𝜃, 𝜙), we numerically integrate the following quantities: 𝑉𝛹=∫𝛹 𝑟3 3sin𝜃d𝜃d𝜙, (3.20) 𝑆𝛹=∫𝛹 1 ˆ𝒏·ˆ𝒓𝑟2sin𝜃d𝜃d𝜙, (3.21) ⟨𝑔⟩=1 𝑆𝛹∫𝛹 𝑔 ˆ𝒏·ˆ𝒓𝑟2sin𝜃d𝜃d𝜙, (3.22) 𝑔−1=1 𝑆𝛹∫𝛹 1 𝑔ˆ𝒏·ˆ𝒓𝑟2sin𝜃d𝜃d𝜙, (3.23) 𝑖rot =∫𝛹(𝑥2+𝑦2)d𝑚𝛹 ∫𝛹d𝑚𝛹 =∫𝛹𝑟2sin2𝜃𝑔−1d𝛹d𝑆 𝑆𝛹𝑔−1d𝛹, =1 𝑆𝛹𝑔−1∫𝛹 𝑟2sin2𝜃 𝑔ˆ𝒏·ˆ𝒓𝑟2sin𝜃d𝜃d𝜙, (3.24) where ˆ𝒏·ˆ𝒓=d𝛹 d𝑟.|∇𝜳| and we used the notation d 𝑚𝛹=𝜌 d 𝑛 d 𝑆 to signify a small mass element on an equipotential shell, such that d 𝑀𝛹=∫𝛹 d 𝑚𝛹 . From these, 𝑓𝑃 and 𝑓𝑇can be computed for each equipotential shell with Eqs. (3.9) and (3.14). 3.2.6 Integration results We consider the component of a Roche binary at the origin with mass 𝑀1 and a companion with mass 𝑀2 , and we computed all the necessary integrals of Eqs. (3.20) - (3.24) in the two-dimensional parameter space of mass ratio 𝑞≡𝑀2/𝑀1 and Roche equipotentials, 𝛹0 . We emphasize here that we calculate solely the properties of the component with mass 𝑀1 and that the properties of the companion 𝑀2 can be obtained by inverting the mass ratio. We sampled 280 equally spaced mass ratios in logarithmic space in the interval log𝑞∈ [− 7 , 7 ] , as well as the equal mass ratio log𝑞= 0 case while in potential space, we sampled 158 values in several logarithmic ranges from 𝛹0= 50 𝛹L1 to 𝛹0=1 2(𝛹Lout +𝛹𝐿4) of different densities. With L out we denote the outer Lagrangian point of the considered component of the binary (always having 𝑥Lout < 0, which is L 3 if 𝑞≤ 1 and L 2 if 𝑞≥ 1). As quantities exhibit important variations around the overflow potentials 𝛹L1 and 𝛹Lout , those will be the regions we sample more densely. We note that𝛹L1< 𝛹L2≤𝛹L3< 𝛹L4=𝛹L5, as in our convention𝛹is strictly negative. METHODS 91 Figure 3.2: Structure correction factor 𝑓𝑃 as a function of fractional RL radius 𝑟𝛹/𝑅RL for various mass ratios. The discontinuities in all graphs correspond to the crossing of the Lagrangian points L 1 (at 𝑟𝛹/𝑅RL = 1, by construction) and L out at larger radii. For reference, the equivalent single rotating star inferred 𝑓𝑃 is overplotted in dashed lines. Figure 3.2 shows the integration results of 𝑓𝑃 for various mass ratios as a function of their radius relative to the RL radius 𝑟𝛹/𝑅RL , where we defined 𝑟𝛹L1≡𝑅RL . We observe discontinuities at the Lagrangian overflow points, which is a result of the splitting surfaces limiting the considered volume to one component (Fig. 3.1). For a comparison with the single rotating star deformation model, we overplot the equivalent 𝑓𝑃 values as if the binary star was interpreted as a single rotating star with the same volume and rotational velocity (see Appendix C.2 for further details). 3.2.7 Comparison to literature As verification of our integrations, we compare our results to calculations performed by Mochnacki (1984). As this study considered vertical splitting surfaces through L 1 at all mass ratio’s, for contact layers we are limited to comparing the mass ratio unity case. Figure 3.3 shows excellent agreement between our 𝑓𝑃 values and those computed via Eq. (3.9) using the results from 92 MODELING TIDAL DISTORTION IN 1D STELLAR STRUCTURE Figure 3.3: Comparison of the calculations of the structure correction factor 𝑓𝑃 from Mochnacki (1984) and those in this work, as a function of the fractional RL radius 𝑟𝛹/𝑅RL . We performed a resolution convergence test with double, quadruple and octuple the default spatial resolution, and conclude our integrations have well converged at an eight-fold resolution increase. Mochnacki (1984), with the relative difference being less than about 1 part in 10 3 . A resolution convergence test of our integrations show our calculations converge to the 10 −5 level at eight times our original resolution. The final computations of all integrals for use in stellar evolution instruments are done with this eight-fold resolution increase, and corresponds to dividing the interval 𝜙=[ 0 , 𝜋] in 𝑛𝜙= 5280 parts and cos𝜃=[0,1]in 𝑛𝜃=𝑛𝜙/2 parts so that 𝛥𝜙 =𝜋/𝑛𝜙and 𝛥cos 𝜃=2/𝑛𝜙. 3.3 The 𝛺𝛤 limit Thanks to the radiation pressure of the photon flux generated in the core, a star gains radiative support against its self-gravity. In the situation where these effects are canceled out exactly, the star has reached the Eddington limit, which is given as (for a single, non-rotating and thus spherical star): 𝐿=𝐿Edd ≡4𝜋𝑐𝐺𝑀 𝜅.(3.25) THE 𝛺𝛤 LIMIT 93 Defining then the Eddington factor 𝛤as: 𝛤≡𝐿rad 𝐿Edd ,(3.26) with 𝐿rad the radiative luminosity, equal to 𝐿 in a radiative envelope, 𝛤= 1 is then an equivalent statement of the Eddington limit. If 𝛤 > 1, the sum of radiative and gravitational acceleration is directed outward. In the interior of the star this can in principle be compensated by an inversion of the gas pressure (Joss et al.,1973), but if a star approaches 𝛤= 1 at the surface, strong outflows are expected to develop (Gräfener et al.,2011). After adding a centrifugal contribution to the force balance, we have: 𝒈tot =𝒈grav +𝒈rot +𝒈rad =𝒈+𝒈rad =0,(3.27) where we write 𝒈=∇𝜳 as the effective gravity. Langer (1997) then considered the Eddington limit on the equator of a rotating (but non-deformed) star, and, assuming that the radiation field is isotropic, the radiative flux is F=𝐿/ 4 𝜋𝑅2 , so that the force balance reduces to 1−𝜔2−𝛤=0,(3.28) where we have defined 𝜔=𝛺/𝛺crit,class as the fractional classical critical rotation rate 𝛺crit,class =√︃𝐺𝑀 𝑅3 . It is important to make the distinction here as the formal critical velocity is lowered to: 𝑣crit =√︂𝐺𝑀 𝑅(1−𝛤)=√︁𝑣crit,class(1−𝛤).(3.29) Therefore, the closer a star is to the Eddington limit, the lower the critical velocity will be. Maeder & Meynet (2000) later noted that this treatment does not take von Zeipel’s theorem into account, namely that the radiative flux is dependent on local effective gravity, Eq. (3.11) . Expressing this for a rotating star to lowest order in rotational velocity 𝛺gives 𝒈tot =𝒈(1−𝛤𝛺(𝜗)),(3.30) where 𝜗 signifies a polar dependence and 𝛤𝛺≡𝛤×𝑓(𝛺, 𝜗) is now a rotationdependent Eddington factor. Since the strength of the radiative flux is a function of the effective gravity through the von Zeipel effect, Maeder & Meynet (2000) showed that, neglecting the 𝜗 dependence, Eq. (3.30) bifurcates when 𝛤= 0 . 639. Below this, 𝒈= 0 is the only solution, so that, independent of luminosity, the break up velocity is the classical velocity 𝑣crit,class =√︁𝐺𝑀/𝑟e , with 𝑟e the equatorial radius of the rotating star. Above the bifurcation point, both 𝒈= 0 and, importantly, 1−𝛤𝛺=0,(3.31) 94 MODELING TIDAL DISTORTION IN 1D STELLAR STRUCTURE determine the critical velocities, the second of which is lower than the classical critical velocity. The lowest of these two velocities then determines the physical break-up velocity of a rotating star. We further expand on the notion of break-up limits by considering the case of deformed stars due to a conservative potential. As a starting point, we consider the same criterion of zero net acceleration, 𝒈tot = 0 , as the stability boundary. The contributions to this acceleration are the effective gravity, 𝒈 , as a result of the potential, and the radiative acceleration. Using von Zeipel’s expression for radiative flux (Eq. 3.11), combined with the equation of hydrostatic equilibrium results in: 𝒈rad =−4𝑎𝑇3 3 d𝑇 d𝑃𝒈,(3.32) which, upon using Eq. (3.13) as the radiative gradient in the non-spherical geometry, gives the stability criterion: 𝒈1−𝛤𝑓𝑇 𝑓𝑃=0.(3.33) We recognize however that von Zeipel’s law of radiation can break down close to critical velocity due to significant baroclinicity, that is, departure from shellularity (see, e.g.,Espinosa Lara & Rieutord,2011). Continuing on our assumption of shellularity, Eq. (3.33) is of a very similar form as that of Maeder & Meynet (2000), namely, 𝒈( 1 −𝛤 𝛹)= 0, where 𝛤 𝛹=𝛤×𝑓 is again a product of the classical Eddington factor with a function dependent on the geometry. For the single rotating star deformation, 𝑓=𝑓(𝜔) , dependent on the fractional critical rotation rate, while in the tidal deformation case the function 𝑓 is of the form 𝑓=𝑓(𝑟𝛹/𝑅RL, 𝑞) , dependent on the mass ratio and degree of Roche filling. In Eq. (3.33) , the classical solution 𝒈= 0 is present of course, where the effective gravity vanishes as a result of rotational (and tidal) support only. As above, in this case there is no corresponding Eddington luminosity as it is independent of any radiative acceleration due to the von Zeipel effect. The other solution is given by: 1−𝛤𝑓𝑇 𝑓𝑃 =0,(3.34) which translates to a modified expression of the Eddington luminosity: 𝐿Edd,𝛹 =𝐿Edd 𝑓𝑃 𝑓𝑇 ,(3.35) or the Eddington factor (see also Sanyal et al.,2015): 𝛤 𝛹=𝐿𝛹 𝐿Edd,𝛹 =𝜅𝐿𝛹 4𝜋𝑐𝐺𝑀𝛹 𝑓𝑇 𝑓𝑃 .(3.36) THE 𝛺𝛤 LIMIT 95 Figure 3.4: Maximal Eddington factor 𝛤max that a rotating star can have as a function of fractional classical critical rotation rate of 𝜔=𝛺/𝛺crit,class . Our model of shellular stars has a maximal reduction to 𝛤max =0.639. These expressions, as in the case of the equations derived in Sect. 3.2.1, are applicable in any conservative potential𝛹. As the criteria for instability, Eqs. (3.28) , (3.31) , and (3.34) all express a departure of the Eddington limit of a star from the classical Eddington limit of 𝛤= 1. In the model developed here, it amounts to a reduction of the maximal Eddington factor by the ratio 𝑓𝑃/𝑓𝑇 , which is smaller than one in both the single rotating star or synchronized binary case. In Fig. 3.4, we show the comparison of this reduction from the results of Langer (1997) to this work in the case of a single rotating star. Moreover, consistent to Maeder & Meynet (2000), we find a bifurcation in the critical rotation velocity at at 𝛤= 0 . 639. Below this number, only the classical critical velocity, 𝜔= 1, will make the star unstable, while above it, this velocity is reduced. The case of a synchronized Roche binary is shown in Fig. 3.5. For example, a binary with a mass ratio of unity, overflowing to its outer Lagrangian point, has its Eddington limit reduced to about 84% its classical value. 96 MODELING TIDAL DISTORTION IN 1D STELLAR STRUCTURE Figure 3.5: Reduction, 𝑓𝑃/𝑓𝑇 , of the maximal Eddington factor as a function of mass ratio, 𝑞 , and fractional overflow radius, 𝑟𝛹/𝑅RL . The equivalent radius of a star filling up to its outer Lagrangian point is marked with the green line. 3.4 Atmospheric boundary conditions The equations of stellar structure must be supplied with appropriate boundary conditions (BCs) at the center and atmosphere of the star. In the center, the conditions are such that 𝑟=𝑀=𝐿= 0, which transform to 𝑟𝛹=𝑀𝛹=𝐿𝛹= 0 in our equipotential shell model. At the surface, an atmosphere model is integrated to supply the pressure 𝑃surf and temperature 𝑇surf . A common choice is the so-called gray, Eddington atmosphere, which is a plane-parallel model integrated from 𝜏=0 to some predetermined atmospheric optical depth 𝜏surf, and is presented in, for example, Cox & Giuli (1968, henceforth CG). In spherical stars, the surface optical depth 𝜏surf is reached at the same physical perpendicular depth, 𝑧 , at all points in the star. For deformed shellular stars, however, the optical depth 𝜏 across an equipotential depends on local effective gravity, as follows: d𝜏=𝜅𝜌 d𝑧=−𝜅𝜌 d𝑛 d𝛹d𝛹=−𝜅𝜌𝑔−1d𝛹, (3.37) with 𝜅 and 𝜌 as the opacity and density of the equipotential, respectively. Consequently, equipotential surfaces do not coincide with surfaces of constant optical depth. For example, in single rotating stars, a fixed 𝜏surf is reached at a deeper equipotential on the poles compared to on the equator, resulting in a ATMOSPHERIC BOUNDARY CONDITIONS 97 different observed effective temperature between those regions (which can further be modulated by a gravity darkening law). The same is true for highly deformed binary stars like contact systems. The atmosphere properties output by stellar evolution codes should therefore also take into account this effect of position dependent atmospheric depth. Within the treatment of deformed stars in a conservative potential, we proceeded as follows. We start by defining the global effective temperature 𝑇eff,𝛹 of an equipotential by the law of Stefan-Boltzmann: 𝐿𝛹=∫𝛹 F·d𝑺≡𝑆𝛹𝜎𝑇4 eff,𝛹 ,(3.38) that is, the global effective temperature of an equipotential shell is that of a black body radiating a luminosity 𝐿𝛹 over an area 𝑆𝛹 , which need not be spherical. Next, using the equation of radiative flux (Eq. 3.11), we can define a local effective temperature 𝑇eff,ℓ that varies across the equipotential: F=𝑔 ⟨𝑔⟩𝜎𝑇4 eff,𝛹 =𝜎𝑇4 eff,ℓ .(3.39) We now consider a gray, plane-parallel, Eddington-approximated atmosphere to construct appropriate BCs. At each point on the surface equipotential, the temperature profile can be written as (CG): 𝑇4(𝜏)=1 2𝑇4 eff,ℓ 1+3 2𝜏.(3.40) This expression is derived using two standard assumptions that can be similarly applied to the case of deformed stars. First, it is assumed that the radiation pressure at all optical depths can be written as 𝑃𝑟(𝜏)=1 3𝑎𝑇4(𝜏) . Second, the intensity of radiation on top of the atmosphere is assumed to be isotropic in the outward directions, meaning that we ignore limb darkening, so that the pressure at the 𝜏=0 surface is computed to be (CG): 𝑃𝑟(0)=2F 3𝑐.(3.41) The first BC we consider constrains the temperature 𝑇surf of the outermost boundary by requiring that 𝑇surf =𝑇eff,𝛹 . To find the second BC on the surface pressure 𝑃surf, we consider the equation of hydrostatic equilibrium: 𝑑𝑃 𝑑𝜏 =𝑔 𝜅.(3.42) Taking a point on the surface equipotential where 𝑇eff,ℓ =𝑇eff,𝛹 , then by Eq. (3.39) , we have that 𝑔=⟨𝑔⟩ . As we also require that 𝑇surf =𝑇eff,𝛹 , Eq. (3.40) implies that the surface optical depth is 𝜏surf = 2 / 3. The surface pressure can then be integrated 104 MODELING TIDAL DISTORTION IN 1D STELLAR STRUCTURE Figure 3.9: Effective temperature evolution of the 32 𝑀⊙ contact model with the different assumptions on BC and rotation model. We notice that the models including rotation are up to 5% cooler during contact. Between the different rotation models and BCs, we find only minor differences. The top panel shows the RL evolution, which explains the difference in age at onset of overflow of the Lagrangian points through spin-orbit coupling. of the stars, this process causes the so-called Darwin instability (Darwin,1879; Hut,1980), where the stars catastrophically spiral in on each other on the orbital timescale. 3.6 Discussion and conclusions We developed the methodology needed to account for tidal deformation in binary stars in 1D stellar evolution codes. They are represented in the structure-correction factors of 𝑓𝑃 and 𝑓𝑇 that multiply the equations of hydrostatic equilibrium and radiative energy transport, respectively. Additionally, modified expressions of the Eddington limit and atmospheric boundary conditions as a result of the deformation have been calculated. We show that the radii predicted by stellar-evolution models can shift by ∼ 5% if centrifugal deformation is included, while also including tidal deformation leads to a smaller ∼ 1% effect. This means comparing observed stellar parameters to non-rotating models could result in a mismatch and this could have implications DISCUSSION AND CONCLUSIONS 105 for high-precision astrophysics. In a study where observed surface properties with uncertainties on the order of 5% are used, centrifugal deformation has to be included in the stellar models. In binaries, if the precision of measurement is better than 1%, then the tidal deformation has to be taken into account. We note that the effect observed here is only due to geometrical considerations, that is, the equator bulges out due to rapid rotation or the presence of a companion. Mixing as a result of rapid rotation was not included and will further affect the radius evolution (and age). Generally, the effect of rotational mixing is to keep the stars compact as fresh fuel is introduced into the core, so this may somewhat cancel the geometrical effect. However, the exact interplay and outcome of these competing effects is far from established. The methods developed here have a wide range of applicability as they can be used for any (semi-)detached binary of arbitrary mass ratio. This can be applied, for instance, in studying the effect of tidal forces in mass transfer stability, which can have a potential impact on the formation of gravitational-wave sources through stable mass transfer (van den Heuvel et al.,2017;Marchant et al.,2021). There are, however, limitations to our model, most notably that the structure-correction factors computed from the binary Roche potential only apply to fully tidally locked systems. The general case, where the rotation of the components is free, cannot be described with a conservative potential, and is thus not suitable to the treatment of Kippenhahn & Thomas (1970) and Endal & Sofia (1976) used here. Developing a global model of free rotation including tides is a complex problem, but a potential workaround would be to consider free rotation for stars well within their RLs and use the centrifugal deformation corrections of single stars, and then implementing a switch to the synchronized binary corrections once a star fills an appreciable amount of its RL. For contact configurations, currently only the mass ratio of unity can be considered consistently, as energy transfer in the contact layers, is zero. When we include the processes of mass and energy transfer in contact layers, which is the subject of the next chapter, the methods developed here enable the consistent modeling of all contact binary systems. Chapter 4 Modeling energy transfer in contact binaries “E.T. phone home...” — E.T., the Extra-Terrestrial This chapter is mainly based on: Modeling contact binaries, II. The effect of energy transfer M. Fabry, P. Marchant, N. Langer and H. Sana ASTRONOMY & ASTROPHYSICS, 672, A175 (2023) Authors contributions: M. Fabry did the majority of the work appearing in this chapter. M. Fabry, P. Marchant and N. Langer jointly developed the methodology, while H. Sana provided critical feedback throughout the project. The text was written by M. Fabry, with additions from N. Langer, and subsequently improved on by all co-authors. With respect to the original publication, some of the text was altered to fit the chapter format of this thesis. Original Abstract: Context. It is common for massive stars to engage in binary interactions. In close binaries, the components can enter a contact phase, when both stars simultaneously overflow their respective Roche lobes. While observational constraints on the stellar properties of such systems exist, the most detailed stellar 107 108 MODELING ENERGY TRANSFER IN CONTACT BINARIES evolution models that feature a contact phase are not fully reconcilable with those measurements. Aims. We aim to consistently model the contact phases of binary stars in a 1D stellar evolution code. To this end, we have developed a methodology to account for energy transfer in the common contact layers. Methods. We implemented an approximative model for energy transfer between the components of a contact binary based on the von Zeipel theorem in the stellar evolution code MESA . We compared structure and evolution models both with and without this transfer. We then analyzed the implications for the observable properties of the contact phase. Results. Implementing energy transfer helps in eliminating baroclinicity in the common envelope between the components of a contact binary, which (if present) would drive strong thermal flows. We find that accounting for energy transfer in massive contact binaries significantly alters the mass-ratio evolution and can extend the lifetime of an unequal mass ratio contact system. 4.1 Introduction The modeling of energy transfer (ET) goes back to the ’60s in an effort to resolve Kuiper’s paradox. Lucy (1968) first considered ET in common convective envelopes, and computed the first approximated structure models for W UMa stars. Later, Lucy (1976); Flannery (1976); Hazlehurst (1985) and Kähler (1989) calculated models that are out of thermal equilibrium and show cyclic behavior. Shu et al. (1976,1979) and Lubow & Shu (1977,1979) (collectively SLA) constructed models of contact binaries where they dropped the requirement of a continuous structure at the layer coinciding with the equipotential surface of the first Lagrangian point ( 𝐿1 ). These models resolve Kuiper’s paradox regardless of the thermal structure of the envelopes, namely, whether they are radiative or convective. Despite much criticism (see e.g., Hazlehurst,1993;Kähler,1989), this is the simplest model of ET in radiative envelopes available in the literature. Accurate modeling of massive, long-lived contact systems has been attempted in the past. Marchant et al. (2016) computed detailed evolution model grids of massive binaries with initial periods down to 0 . 5d, which includes the regime of contact binaries at the zero-age main sequence. Those models however did not include ET between contact components, as the models concerned binaries of mass ratio close to unity, 𝑀2/𝑀1= 0 . 8 – 1, and the effect was thought to be minimal. Sen et al. (2022) used the models of Marchant (2018) to study the semi-detached Algol systems, although these models also included contact phases. Menon et al. (2021) computed models with initial mass ratios down to 𝑀2/𝑀1= 0 . 6 with the THEORY OF ENERGY TRANSFER 109 study of massive, long-lived contact binaries in mind. They found, also without including ET in contact phases, a strong correlation between observed mass ratio and period in contact systems that is broadly in agreement with observations. However, the mass-ratio distribution they derive is heavily skewed toward values close to unity, which is not supported observationally. These authors suggested that including ET in contact phases of unequal mass components could alleviate this discrepancy. In this chapter, we apply the ET model of SLA in common stellar layers to modern stellar structure models to further advance our evolutionary modeling of massive contact binaries. In Sect. 4.2, we describe and discuss the theory of ET. Section 4.3 describes the physical setup of the stellar evolution code, along with our ET implementation. In Sect. 4.4, we compare models computed with and without ET in contact layers and discuss the difference in the observable properties of the models. Lastly, in Sect. 4.5, we provide our concluding remarks. 4.2 Theory of energy transfer 4.2.1 Simple considerations Theoretical modeling of contact stars started with Kuiper (1941), who stated that stable contact systems of uniform composition with unequal masses cannot exist as a result of differing mass-radius relationships. For the galactic ZAMS models of Brott et al. (2011), we derive the mass-radius relation of single stars to be approximately: 𝑅2 𝑅1ZAMS =𝑀2 𝑀10.57 =𝑞0.57.(4.1) However, the condition of contact in a binary following the Roche geometry constrains the surface of the stars to the same equipotential, which leads to: 𝑅2 𝑅1Roche ≈𝑀2 𝑀10.46 =𝑞0.46,(4.2) for stars not overflowing their RL too much. Clearly, Eqs. (4.1) and (4.2) cannot be satisfied simultaneously unless 𝑀1=𝑀2 . However, while the contact condition needs to be satisfied from dynamical arguments, the stellar structure does not, a priori, need to follow a single-star model. Even though the dense stellar core will be largely unperturbed due to the companion, the outer layers of contact components are highly distorted, causing variations in the total radii with respect to spherical models, as seen explicitly in Chap. 3. This result had not yet taken the possibility of ET into account and we expect further changes under the 110 MODELING ENERGY TRANSFER IN CONTACT BINARIES consideration that a hotter gas under similar pressure takes up more volume. Therefore, the ZAMS mass-radius relation in Eq. (4.1) is not expected to be satisfied under general contact conditions, and so Kuiper’s paradox can be resolved by providing alternative stellar models that have a mass-radius relation closer to the contact condition in Eq. (4.2). Following the Roche-lobe geometry, Lucy (1968) finds the ratio of the surface areas 𝑆2/𝑆1 of two stars in contact to be proportional to (𝑀2/𝑀1)𝛽 , with 𝛽= 0 . 96. Using the approximation 𝛽≃ 1, combined with von Zeipel’s theorem of gravity darkening, 𝑇4 eff ∝𝑔 (von Zeipel,1924), and the Stefan-Boltzmann law, 𝐿∝𝑆𝑇4 eff , results in the simple expectation that in contact binaries, the luminosity ratio follows the mass ratio (Lucy,1968;Tassoul,2000): 𝐿2 𝐿1≃𝑀2 𝑀1 =𝑞. (4.3) Single main-sequence stars, on the other hand, follow the well-known massluminosity relation: 𝐿s,2 𝐿s,1≃𝑀2 𝑀1𝛼 =𝑞𝛼,(4.4) with 𝛼≃ 2 – 3 for the upper main sequence (Gräfener et al.,2011;Köhler et al.,2015). We use the symbol “ ≃ ” to denote these relations are approximations to simple power laws. We see that this leads to a difference between the luminosity of a single star 𝐿s,1 and the luminosity 𝐿1 of a star of the same mass in a contact binary of 𝛥𝐿1=𝐿1−𝐿s,1=−𝑓 𝐿1,(4.5) and for the companion 𝛥𝐿2=𝐿2−𝐿s,2=𝑓 𝐿1≃𝑓 𝑞𝐿2,(4.6) since 𝐿2≃𝑞𝐿1 and we require 𝛥𝐿1+𝛥𝐿2= 0 to conserve energy. This defines 𝑓 as 𝑓≃𝑞−𝑞𝛼 1+𝑞𝛼.(4.7) Therefore, the two stars in a contact binary can fulfill the single-star massluminosity relation in Eq. (4.4) in their cores and the contact binary massluminosity condition in Eq. (4.3) at their surfaces if the amount of energy per time, given by Eqs. (4.5) and (4.6) , is transferred from the more massive to the less massive star in their common envelope. THEORY OF ENERGY TRANSFER 111 4.2.2 Models of energy transfer The general solution to Kuiper’s paradox is to consider detailed stellar models with the inclusion of ET between the binary components. Several models of ET are given in the literature. Lucy (1968) and Biermann & Thomas (1972) provided a first solution by adjusting the adiabatic constants of convective envelopes in contact components. However, this is an unsatisfactory solution, since this setup requires the stars to be burning hydrogen through different nuclear chains or cycles in the case of Lucy (1968) or that the models exhibit inaccurate light curves, as in Biermann & Thomas (1972). Other models, such as those of Lucy (1976) or Flannery (1976), relaxed the requirement of thermal equilibrium and constructed models of W UMa stars that exhibited thermal cycles. Kähler (1989) presented a detailed model that required turbulent motions in the common envelope to explain early-type (radiative) W UMa binaries. Meanwhile, SLA presented the contact discontinuity model of contact binaries, by relaxing the requirement of continuous structural quantities across the RL. This is the only model that treats the common envelope as a single volume of the binary structure, at the price of hiding a heat engine in a very thin region around the RL. One peculiar feature is that this model necessitated a temperature inversion at the RL layer in one of the components as otherwise they would not be able to construct thermally stable contact models of uniform composition (in order to model binaries at zero age). This feature has received criticism in that the proposed heat engine violates the second law of thermodynamics and cannot be stable over thermal timescales (see Hazlehurst,1993;Kähler,1989, and references therein). Eventually, Kähler (2004) drew the following conclusion based on all collected theoretical arguments: internal circulation currents must exist in the less luminous component to reduce the radiative temperature gradient, since the luminosity carried by radiation is reduced by the circulation luminosity. Given the complexity of the theoretical problem of the structure of contact binaries, especially with radiative envelopes, it is beyond the scope of this work to further develop the analytic theory. Instead, we apply an ET model in modern stellarstructure calculations. Using shellularity as our base assumption of the stellar structure (see Sect. 4.3.2 for the precise definition), the model of SLA is a natural choice, although we recognize that this comes with the apparent thermodynamical problems stated above. However, we believe we avoid the most fundamental one, as we do not explicitly require a temperature inversion in our models. We only used the model of SLA to compute the amount of energy transferred (see Sect. 4.2.3). 112 MODELING ENERGY TRANSFER IN CONTACT BINARIES 4.2.3 Energy transfer in the Roche geometry The work of SLA provides a general model of ET by introducing the notion of an energy flow at the base of the common envelope. If the contact layers are shellular, and they satisfy von Zeipel’s gravity darkening, conservation of energy at the RL implies: 𝐿′ 1=(𝐿1+𝐿2)𝑆1⟨𝑔⟩1 𝑆1⟨𝑔⟩1+𝑆2⟨𝑔⟩2 ,(4.8a) 𝐿′ 2=(𝐿1+𝐿2)𝑆2⟨𝑔⟩2 𝑆1⟨𝑔⟩1+𝑆2⟨𝑔⟩2 .(4.8b) Here, 𝑆 is the surface area of the RL, ⟨𝑔⟩ is the surface averaged effective gravity at the RL, and the primed quantities specify the state just above the ET layer, while unprimed those just below. This equation specifies that the fraction of the total luminosity that each component radiates is proportional to 𝑆⟨𝑔⟩ . The transferred luminosity then equals 𝐿trans =𝐿1−𝐿′ 1=𝐿1𝑆2⟨𝑔⟩2−𝐿2𝑆1⟨𝑔⟩1 𝑆1⟨𝑔⟩1+𝑆2⟨𝑔⟩2 .(4.9) Comparing Eqs. (4.8) against Eqs. (4.5)-(4.6), we find for the fraction 𝑓: 𝑓=𝐿1𝑆2⟨𝑔⟩2−𝐿2𝑆1⟨𝑔⟩1 (𝐿1+𝐿2)𝑆1⟨𝑔⟩1 ,(4.10) which is consistent with Eq. (4.7) as 𝑆2⟨𝑔⟩2 𝑆1⟨𝑔⟩1≃𝑞and 𝐿2 𝐿1≃𝑞𝛼. 4.3 Methods To investigate the effect of ET on the evolution of contact binaries, we computed binary evolution models using the stellar evolution code MESA (Paxton et al., 2011,2013,2015,2018,2019;Jermyn et al.,2023), version r22.11.1. These models are the first massive binary evolution models that include ET in contact layers. We followed the evolution of binaries from the ZAMS until the least massive component overflows the second Lagrangian point. 4.3.1 Physical assumptions in MESA The microphysical setup of our stellar evolution models remains mostly the same as in Chap. 3, with some additions from the newer MESA version. We use the same nuclear microphysics as in Chap. 3. The EOS has an additional table available from METHODS 113 the Skye project (Jermyn et al.,2021), its blend specified in Jermyn et al. (2023). Radiative opacities are also blended from CO-enhanced tables of OPAL opacities (Iglesias & Rogers,1993,1996) and Ferguson et al. (2005) for lower temperatures. At high temperatures, Compton scattering opacity is from Poutanen (2017), and electron conduction opacities are taken from Cassisi et al. (2007) and Blouin et al. (2020). Like Chap. 3, we set the metallicity of stars to the solar value, 𝑍⊙ , where 𝑍⊙=0.0142, with metal fractions as determined from Asplund et al. (2009). Mass loss through winds is accounted for by following the prescription of Brott et al. (2011). If the surface hydrogen fraction is 𝑋 > 0 . 7, the mass loss rate is taken either from Vink et al. (2001), for temperatures above the iron bi-stability jump (calibrated also by Vink et al.,2001), or the maximum of the rates from Vink et al. (2001) and Nieuwenhuijzen & de Jager (1995), below this temperature. For 𝑋 < 0 . 4, the wind prescription of Hamann et al. (1995) is used, albeit decreased by a factor of ten. When 0 . 4 < 𝑋 < 0 . 7, the wind is linearly interpolated between the above results. We allow part of the wind launched by a star to be accreted by its companion using the Bondi-Hoyle mechanism (Bondi & Hoyle,1944) as implemented by Hurley et al. (2002). For internal mixing and rotation, we use the same setup as in Sect. 3.2. The treatment of mass transfer (MT) is explained in Sect. 4.3.3, while the implementation of energy transfer (ET) is shown in Sect. 4.3.4. 4.3.2 Shellularity and Roche lobe geometry Since we deal with highly tidally deformed stars, we use the modifications to the stellar structure equations from Sect. 3.2.1 to incorporate the RL geometry into a one-dimensional (1D) stellar evolution code. The stars are therefore modeled as hydrostatic structures living in the Roche potential 𝛹 of a fully synchronized binary. Furthermore, we assume stellar layers to be shellular. Shellularity is reached when all intensive quantities (in particular the temperature, pressure, and mass density) are constant along a stellar layer and when such a layer coincides with a unique equipotential surface. It should be emphasized that in the context of 1D models, the full shellularity of layers of a (single) star is an assumption – not a self-consistently modeled feature as, of course, there is no 3D structure to explicitly test the shellularity of a stellar layer. However, overflowing layers of components in a contact binary can be tested whether they are shellular with respect to each other if the temperature, density, and so on at the 1D cells are equal for equal values of Roche potential. Given the method we used to split the common envelope of the contact binary (Sect. 3.2.4), the computation of the tidal deformation corrections based on equipotential surfaces, and the usage of the corresponding outer boundary condition (Sect. 3.4), 120 MODELING ENERGY TRANSFER IN CONTACT BINARIES Figure 4.3: Contact duration distributions as function of observed mass ratio ¯𝑞=min 𝑀1 𝑀2,𝑀2 𝑀1 .Top: Differential duration (in bins of 𝛥¯𝑞= 0 . 005). Bottom: Cumulative duration in configurations more extreme than ¯𝑞. seen in the left column of Fig. 4.4 by the joining of the temperature and density profiles of the primary and secondary. In particular, at the surface, the density of the ET components agree to within 0 . 2%, while the temperature to within 0 . 5%, whereas the surface properties of the no-ET components vary more than 10% (as expected for models of differing mass). While the former differences are hard to constrain from measurement, current data-analysis techniques are able to measure surface temperatures with accuracies better than 10%, so that the no-ET model can be ruled out observationally. A significant shortcoming of our model is the assumed thickness of the ET layer. Shu et al. (1979) argued that the thickness 𝑑of the ET layer is on the order of: 𝑑 𝑎∼𝛿0.4, 𝛿 =𝐿trans/𝑎2 𝜌ℎ𝑐𝑠 ,(4.17) with 𝑎 as the binary separation, and 𝜌 , ℎ, and 𝑐𝑠 as the density, specific enthalpy and local sound speed evaluated at the RL, respectively. Lubow & Shu (1979) computed this number to be 𝑑/𝑎∼ 10 −2 for stars of masses around 4 – 8 𝑀⊙ which justifies SLA in modeling the layer as a discontinuity in the stellar profile, located at the RL radius. STELLAR MODELS 121 Figure 4.4: Profiles of outer layers of the binary components at the onset of inverse MT. It shows (top to bottom) the temperature, density, and luminosity profiles for both components in the ET (red, left column) and no ET (blue, right column) cases, all as a function of the scaled radius coordinate ˜𝑟=𝑟−𝑅RL 𝑅−𝑅RL . The gray vertical lines show the location of the RL. However, direct computation of Eq. (4.17) for our 25 𝑀⊙+ 20 𝑀⊙ model show that this estimation breaks down for higher masses, see the red line in Fig. 4.5. This suggests that the energy redistribution flow, modeled as a discontinuity by SLA, is not sufficient in these higher mass stars. Except when the binary has reached considerable overflow, where 𝑑/𝑎≲ 10 −2 , we find that the thickness needed can be a significant fraction of the binary separation, even surpassing it in the early stages of the contact phase. Another way to compute the thickness of the ET layer is to use Bernouilli’s equation. If we consider fluid motion from far away from L 1 on the primary star (location 𝑖) toward L1(location 𝑓), we have 1 2𝑣2 𝑓+∫𝑓 𝑖 𝑑𝑃 𝜌=0,(4.18) 122 MODELING ENERGY TRANSFER IN CONTACT BINARIES where we have already canceled the potential terms 𝛹𝑖,𝛹𝑓 since we move along an equipotential surface, and the initial velocity, 𝑣𝑖 , is assumed to be negligible. Making the estimation: ∫𝑓 𝑖 𝑑𝑃 𝜌≈1 𝜌𝑓+1 𝜌𝑖𝑃𝑓−𝑃𝑖,(4.19) and assuming that the transferred energy through the binary neck of width, 𝑏 , by a mass flow of ¤ 𝑀=𝜌𝑖𝑣𝑓𝑏2is: 𝐿trans =𝜌𝑖𝑣𝑓𝑏2(𝑐𝑝,1+𝑐𝑝,2)(𝑇𝑓−𝑇𝑖),(4.20) we compute for the minimal thickness of the ET layer: 𝑏≈v u u u t𝐿trans 𝜌𝑖(𝑐𝑝,1+𝑐𝑝,2)(𝑇𝑓−𝑇𝑖)√︂1 𝜌𝑓+1 𝜌𝑖𝑃𝑓−𝑃𝑖.(4.21) Finally the thickness of the layer at the neck is related to the thickness far away from L1via: 𝑑 𝑎≈𝑏 𝑎2 ,(4.22) since the Roche potential varies quadratically near L 1 and linearly elsewhere. Equation (4.21) gives a lower limit on the width of ET layer so that a balanced mass flow, ¤ 𝑀 , in the contact binary can carry a to be transferred luminosity, 𝐿trans . Conversely, it can be interpreted as 𝐿trans being the maximal luminosity the mass flow can carry in a layer of fixed width, 𝑏. We plot the Bernouilli-computed thickness of Eq. (4.22) in Fig. 4.5. Similarly, we see that the required thickness, 𝑑/𝑎 , is much larger than what the radius of the primary star allows room for. It is only at later times, when the mass ratio has equilibrated and deeper contact is engaged, that the estimated width becomes smaller than the overflow rate, 𝑅−𝑅RL, of the primary. 4.4.2 Mass versus luminosity ratios During the nuclear-timescale, inverted MT phase, contact is engaged so that our ET scheme acts to move luminosity from one component to the other. As mentioned in Sect. 4.2, we expect the luminosity ratio of contact binaries to follow the mass ratio, 𝐿∝𝑀 , as opposed to detached stars following a single-star mass-luminosity relation, 𝐿∝𝑀𝛼 , with 𝛼≃ 2 – 3. Figure 4.6 shows the evolution of the luminosity ratio as function of the mass ratio during the slow MT phase of the 25 𝑀⊙+ 20 𝑀⊙ explored in Sect. 4.4.1. In this graph, the models evolve from the top right at STELLAR MODELS 123 Figure 4.5: Thickness of the ET layer, 𝑑 , with respect to the binary separation, 𝑎 , as a function of the radius of the primary star during the long lived contact phase. The gray line gives the physical size of the overflowing layers, and corresponds to the maximal width the ET layer can assume. 𝑞≈ 1 . 4 to near-equal mass ratio on the left. We see that the model not including ET follows closely a 𝑞2.2 relation, appropriate for single stars in the mass range of 10 – 30 𝑀⊙ . As the models including ET engage into deep contact however, their mass-luminosity ratio changes drastically from the 𝑞2.2 line in near-contact to the 𝑞1relation in full contact. We give (overplotted on Fig. 4.6) the measurements of several observed massive contact binaries of Abdul-Masih et al. (2021); Yang et al. (2019) and Lorenzo et al. (2014) (see also Fig. 2 of Langer,2022). Curiously, the luminosity ratios from Abdul-Masih et al. (2021), while they do follow an 𝐿∝𝑀 trend, are offset to lower 𝐿2/𝐿1 than predicted. Either the luminosity of the primary is overestimated or the uncertainties are underestimated. The systems included from Mahy et al. (2020) were categorized as ‘uncertain configurations’ since the measurement of the radius was consistent with being above as well as below the RL. In the context of the mass-luminosity relation however, we expect that the system from Mahy et al. (2020) at 𝑞≈ 1 . 3 (VFTS 563) is a true contact system, while the one at 𝑞≈ 1 . 2 (VFTS 217) is not (although within 1𝜎it could be either). 124 MODELING ENERGY TRANSFER IN CONTACT BINARIES Figure 4.6: Luminosity ratio versus mass ratio during the nuclear-timescale MT of the 25 𝑀⊙+ 20 𝑀⊙ systems of Sect. 4.4.1. Measurements from observed massive (near-)contact systems from Abdul-Masih et al. (2021), Yang et al. (2019), and Lorenzo et al. (2014) are overplotted. The systems from Mahy et al. (2020) were classified as ‘uncertain configurations’. 4.5 Conclusions In this chapter, we have taken a step forward in the detailed modeling of contact binaries, by implementing a model of ET in detailed stellar-structure and evolution models. From Fig. 4.4, we see that this relatively simple model (i.e., the inclusion of a heat source or sink in the stellar model as proxy for the ET) is capable of making the common layers of stellar components in a contact configuration shellular. Models without such ET do not exhibit their common layers to be shellular, which, from theoretical arguments, would drive strong horizontal flows equilibrating all gradients, in particular pressure. From Figs. 4.2 and 4.3, we see that the time spent in deep contact at mass ratios between ¯𝑞= 0 . 7 – 0 . 8 is extended when ET was included, versus when it was ignored. This is a promising result, in that if this trend persists across the parameter space of the total mass, initial mass ratio, and initial period, ET could provide an answer to the discrepancy between the observed mass-ratio distribution of massive contact systems and its predicted distribution. The computation and analysis of a full grid of models for a population synthesis study is the topic of the next chapter. Chapter 5 Energy transfer in a population of massive contact binaries “I know my apprehensions might never be allayed, and so I close, realizing that perhaps the ending has not yet been written.” — Atrus This chapter is mainly based on: Modeling contact binaries, III. Properties of a population of close, massive binaries. M. Fabry, P. Marchant, N. Langer and H. Sana in preparation for submission to ASTRONOMY & ASTROPHYSICS 5.1 Introduction In Chaps. 3and 4, we have developed the methodology to consistently take into account tidal deformation of close binary components and energy transfer in the common layers of contact binaries. In the introduction, we have also established that the state-of-the-art models have trouble reproducing the observed mass-ratio distribution of massive contact systems. While the models indicate that, once nuclear-timescale contact is engaged, the masses should equalize rapidly (Menon 125 126 ENERGY TRANSFER IN A POPULATION OF MASSIVE CONTACT BINARIES et al.,2021), long-term monitoring of the period derivative shows observed systems are stable on the main-sequence lifetime (Abdul-Masih et al.,2022). Still, in Chap. 4, we showed that for particular systems, energy transfer in radiative envelopes can alter the mass-ratio evolution of massive contact systems, and has the possibility to relieve the discrepancy between models and observations. In this chapter, we compute a grid of short-period, massive binary-evolution models and investigate the effect of energy transfer on the mass-ratio distribution of contact systems. In Sect. 5.2, we describe the modifications to the stellar evolution code with respect to the previous chapters, the initialization of the binary models, termination conditions, and we detail the calculations used for making population predictions. Next, Sect. 5.3 presents and discusses the results obtained from the population synthesis, while Sect. 5.4 gives concluding remarks. 5.2 Methodology We use the detailed binary evolution code MESA , version 22.11.1 (Paxton et al., 2011,2013,2015,2018,2019;Jermyn et al.,2023), to compute grids of binary models that undergo case A mass transfer, from which we assemble a synthetic population. The parameters to be varied in each grid are the initial primary mass, 𝑀1,init , initial mass ratio, 𝑞init =𝑀2,init/𝑀1,init , and the initial period, 𝑝init . For initial primary mass, we sample 𝑀1,init ={ 8 , 9 , 10 , ..., 19 }∪{ 20 , 22 , 24 , ..., 48 } ∪ { 50 , 52 . 5 , 55 , ..., 70 }𝑀⊙ . The initial mass ratio is 𝑞init = 0 . 6 – 0 . 975, spaced uniformly with 𝛥𝑞init = 0 . 025, while the initial period is 𝑝init = 0 . 5 – 8d, spaced logarithmically with 𝛥log(𝑝/ d )= 0 . 04. We compute two grids with these parameter variations: one that includes the effect of energy transfer (ET), while the other ignores it. This amounts to a total of 35712 detailed binary-evolution models. Throughout the discussion, we will refer to the models that include ET as the “ET models,” while those that do not the “no-ET models.” 5.2.1 Physical assumptions in MESA We use a similar setup of the MESA code as in Chap. 4,i.e., we use the same microphysics, convection and overshoot parameters (with two exceptions, see below), wind prescription, deformation geometry, and massand energy-transfer (MT, ET) calculations. In contrast to Chaps. 3and 4we make the following modifications to the modeling setup: • We do not assume rigid-body rotation throughout the evolution, so that differential rotation is allowed. However, we do include fast tidal synchronization by uniformly modifying the angular momentum of the star METHODOLOGY 127 on the timescale of the orbital period, 𝑝, 𝛥𝑗tides =1−exp −𝛥𝑡 𝑝(𝜔orb𝑖rot −𝑗).(5.1) Here 𝛥𝑡 is the timestep of the evolution, 𝑗 and 𝑖rot are the current angular momentum and the moment of inertia of the stellar layer, respectively, and 𝜔orb =2𝜋 𝑝 is the orbital angular velocity. We model spin-orbit coupling to conserve total angular momentum, so that the required angular momentum to synchronize both stars in this fashion is then subtracted from the orbital angular momentum. Given that 𝛥𝑡 is much longer than the orbital period in evolutionary calculations (years to thousands of years versus days), Eq. (5.1) specifies efficient tidal synchronization. We thus expect little departure from solid-body rotation in our very close binary simulations, since processes that cause departures from it operate on longer timescales than the orbital one. • We also allow for rotational mixing in our models. We include EddingtonSweet circulation, the Goldreich-Schubert-Fricke instability and both the dynamicaland secular-shear instability, all following the implementation of Heger et al. (2000). • We reduce the parameters of Mixing Length Theory to 𝛼MLT = 1 . 5, and the semiconvective efficiency to 𝛼sc =1, to be in line with Menon et al. (2021). 5.2.2 Binary initialization Once initial masses and a period are selected, the binary model is initialized as follows. Non-rotating, single-star models of mass 𝑀1,init and 𝑀2,init =𝑞𝑀1,init are interpolated from pre-computed zero-age main sequence (ZAMS) models of solar metallicity ( 𝑍⊙= 0 . 0142, with metal fractions following Asplund et al. 2009) and a helium content of 𝑌= 0 . 2703 and loaded into a binary model. We set the eccentricity to zero as we expect circular orbits for close binaries (Zahn,1975). The single star models are then spun up to the appropriate Keplerian angular velocity 𝜔orb . Since the models are brought out of thermal equilibrium due to the spin up, what follows is a thermal relaxation period onto the true ZAMS of the models, during which MT, ET or any other mass loss is turned off. We also keep the period constant to the initial value and enforce rigid body rotation at the orbital angular velocity. We define the ZAMS as the point where the luminosity of both stars is within 1% of their nuclear luminosity. When this condition is reached, we make a decision depending on the rate of overflow of the stars. First, if both stars are within their RL, the simulation proceeds as normal, without further intervention. Second, if at least one star overflows its RL, this means that MT must have started in the pre-main sequence (PMS) phase. To simulate 128 ENERGY TRANSFER IN A POPULATION OF MASSIVE CONTACT BINARIES this, we artificially boost the orbital period to put the components in a detached configuration, and then, on a thermal timescale of the most massive star ( 𝜏= 0 . 75 𝐺𝑀2 1/𝑅1𝐿1 ), drain the added angular momentum from the system to again reach the initial period. During this process, the components reach the RL, and MT and ET are allowed. Finally, if one of the stars overflows the second Lagrangian point (L2OF) at the ZAMS, we terminate the simulation immediately as we assume a merger would have happened on the PMS. 5.2.3 Outcomes and termination The outcome of our simulations broadly falls in three categories. • The binary reaches L 2 in a contact configuration. In this situation, we assume mass loss will start from L 2 , carrying away angular momentum, which we expect results in a merger on a dynamical timescale. We stop the simulation at this point and classify these systems as “mergers.” A variant of this outcome is when the MT calculation determines that ¤ 𝑀 is a very high value. We choose to take a limiting value of ¤ 𝑀= 10 −1𝑀⊙yr−1 across the whole grid. MT rates higher than this value are at least two orders of magnitude above the thermal timescale MT rates of massive main-sequence stars (except perhaps for the highest masses in our grid, where ¤ 𝑀thermal ≈ 5 × 10 −3𝑀⊙yr−1 ), and we expect MT to evolve toward dynamical timescale MT. • One of the stars (most likely the initially more massive one) reaches core hydrogen exhaustion, without triggering a merger. At this point, the star will start undergoing rapid evolution, and we stop the simulation calling these “(main-sequence) survivors.” • The simulation can encounter numerical difficulties, by failing to reach convergent models during the evolution. In this case, we halt the simulation and signal this system encountered an “error.” 5.2.4 Population synthesis computations When doing population synthesis calculations, each model in our grid needs to be assigned a relative weight to accurately predict the occurrence rate of a certain system in a population. Because the probability, P , of finding systems in a certain configuration 𝜗0 is proportional to the time the system spends in that configuration, we have P(𝜗0) ∝ ∫𝛿(𝜗−𝜗0)WNdN,(5.2) METHODOLOGY 129 where we sum the time each of the models, N , spends at the configuration 𝜗0 , using the Dirac delta function, 𝛿 . The integrand is then multiplied with its statistical weight, WN , coming from the likelihood of this model being created from the star and binary formation processes. Since we vary the initial primary mass, initial mass ratio and initial period, and take into account the star-formation rate (SFR) history, the total (relative) weight is going to be, WNdN=W𝑀1,init d𝑀1,initW𝑞init d𝑞initW𝑝init d𝑝initWSFR d𝑡. (5.3) The individual weights are taken from models of observed young populations. We take W𝑀1,init d𝑀1,init =𝑀−2.35 1,init d𝑀1,init,(5.4a) W𝑞init d𝑞init =𝑞0 init d𝑞init =d𝑞init,(5.4b) W𝑝init d𝑝init =(log 𝑝init)0dlog 𝑝init,(5.4c) WSFR d𝑡=d𝑡, (5.4d) coming from the Salpeter initial-mass function (Salpeter,1955;Kroupa,2001; Bastian et al.,2010) and the modeled mass-ratio distributions found from a sample of O-type stars from Sana et al. (2012), which are shown to be similar for Btype stars (Villaseñor et al.,2021;Banyard et al.,2022). No sample of observed binaries covers periods shorter than 1d, so we adopt Öpik’s law (Öpik,1924), W𝑝init =log 𝑝init , which is close to the findings of Almeida et al. (2017). Unless the index of the period-distribution power-law is positive and high, which is not supported by any observational sample of young, massive stars, short-period binaries are favored over longer ones. We use no twin binary excess in our setup, since high-mass stars have small twin fractions (Sana et al.,2012;Moe & Di Stefano, 2017). The age prior indicates we assume a constant star-formation history. Given we have a grid of a discrete number of models, each grid point represents a finite area of the birth distributions in Eqs. (5.4). In particular, we compute: ∫upper lower 𝑀−2.35 1,init d𝑀1,init ∝𝑀−1.35 1,init,upper −𝑀−1.35 1,init,lower,(5.5) for the contribution of the model of initial primary mass 𝑀1,init , and the upper and lower bounds are taken in the middle between our grid points. For points on the edge of the grid, we linearly extrapolate the spacing outward and then use Eq. (5.5) as usual. In this way, the model with initial primary mass 8 𝑀⊙ has 𝑀upper = 8 . 5 𝑀⊙ and 𝑀lower = 7 . 5 𝑀⊙ , while the model of highest initial mass has 𝑀upper = 71 . 125 𝑀⊙ and 𝑀lower = 68 . 875 𝑀⊙ . A similar computation for the bins of 𝑞init and 𝑝init is performed.