© ESO, 2014

1. Introduction

Binary population synthesis (BPS) codes enable the rapid calculation of the evolution of a large number of binary stars over the course of the binary lifetime. With such models, we can study the diverse properties of binary populations, e.g. the chemical enrichment of a region, or the frequency of an astrophysical event (for a review, see Han et al. 2001). We can learn about and study the formation and evolution of stellar systems that are important for a wide range of astronomical topics: novae, X-ray binaries, symbiotics, subdwarf B stars, gamma-ray bursts, R Coronae Borealis stars, AM CVn stars, Type Ia and Type Ib/c supernovae, runaway stars, binary pulsars, blue stragglers, etc.

To carefully study binary populations, in principle it is necessary to follow the evolution of every binary system in detail. However, it is not feasible to evolve a population of binary stars from the zero-age main-sequence (ZAMS) to remnant formation with a detailed stellar evolution code. Such a task is computationally expensive as there are many physical processes which must be taken into account over large physical and temporal scales, such as tidal evolution, Roche lobe overflow (RLOF), mass transfer. Moreover not all processes can be modelled with detailed codes or are quite uncertain (or both), e.g. common envelope evolution, contact phases. Therefore, simplifying assumptions are made about the binary evolution process and many of its facets are modelled by the use of parameters. This process is generally known as binary population synthesis. Examples of such parametrisation are straightforward prescriptions for the stability and rate of mass transfer. To some degree, the effects that are most important for the problem being studied will be more elaborately included in the corresponding BPS codes. For the evolution of an individual system, the above can of course be an oversimplification. However, for the treatment of the general characteristics of a large population of binaries this process works very well (e.g. Eggleton et al. 1989).

Recently, several BPS codes have been used to study the progenitors of Type Ia supernovae (e.g. Yungelson et al. 1994; Han et al. 1995; Jorgensen et al. 1997; Yungelson & Livio 2000; Nelemans et al. 2001a; Han & Podsiadlowski 2004; De Donder & Vanbeveren 2004; Yungelson 2005; Lipunov et al. 2009; Ruiter et al. 2009b, 2011; Mennekens et al. 2010; Wang et al. 2010; Meng et al. 2011; Bogomazov & Tutukov 2009, 2011; Ruiter et al. 2013; Toonen et al. 2012; Mennekens et al. 2012; Claeys et al. 2013). From these recent studies, it has become evident that the various codes show different results in terms of the SNe Ia rate (Nelemans et al. 2013), in particular for the single degenerate channel in which binary systems can produce a SNe Ia by accretion from a non-degenerate companion to a white dwarf (WD). The differences in the predicted SNIa rate are largely, but not completely, due to differences in the assumed retention efficiency of the accretion onto the WD (Bours et al. 2013). While it has long been expected by groups working on population synthesis that the differences in the BPS results were the result of different assumptions being made in these various studies rather than numerical in nature, it became ever more clear that a quantitative study of the nature and causes of these differences is necessary.

This paper aims to do this by clarifying, for four different BPS codes, the respective ingredients and assumptions included in the population codes and comparing models of several simulated populations for which all assumptions have been made the same as much as possible. We discuss the similarities and differences in the predicted populations and examine the causes for the differences that remain. The causes for differences are valuable information for interpreting BPS results, and as input for the astronomical community to increase our understanding of binary evolution. The project is known as PopCORN – Population synthesis of Compact Objects Research Network. It is not the purpose of this paper to discuss the advantages or shortcomings of the respective methods used in BPS codes, nor to judge which assumptions made for binary evolutionary aspects are the most desirable.

The paper focuses on low and intermediate mass close binaries, i.e. those with initial stellar masses below 10 M_⊙. The reason for this is twofold: firstly, as the project originates from differences in the predictions of SNe Ia rates, the systems that produce WDs are the main focus. Secondly, since the evolution of massive stars is even less straightforward, and its modelling includes even more uncertainties, comparing massive star population synthesis will be a whole new project.

In Sect. 2 we give an overview of the relevant processes for the evolution of low- and intermediate-mass binaries. Sect. 3 describes the codes involved in this project. The method we use to conduct the BPS comparison is described in Sect. 4. We compare the simulated populations of systems containing one WD in Sect. 5.1 and two WDs in Sect. 5.2. A more detailed comparison for the most important evolutionary paths is given in Appendix A. In Sect. 6 we summarise and discuss the causes for differences that were found in Sect. 5. Our conclusions are given in Sect. 7. An overview of the inherent and typical assumptions of each code can be found in Appendices B and C respectively.

2. Binary evolution

In this section we will give a rough outline of binary evolution and the most important processes that take place in low and intermediate mass binaries. The actual implementation in the four BPS codes under consideration in this study is described in Appendices B and C.

Low- and intermediate-mass systems with initial periods less than approximately 10 years and primary masses above approximately 0.8 M_⊙, will come into Roche lobe contact within a Hubble time. The stars in a binary system evolve effectively as single stars, slowly increasing in radius and luminosity, until one or both of the stars fills its Roche lobe. At this point mass from the outer layers of the star can flow through the first Lagrangian point leaving the donor star.

Depending on the reaction of the star upon mass loss and the reaction of the Roche lobe upon the rearrangement of mass and angular momentum in the system, mass transfer can be stable or unstable. When mass transfer becomes unstable, the loss of mass from the donor star will cause it to overfill its Roche lobe further. In turn this increases the mass loss rate leading to a runaway process. In comparison, when mass transfer is stable, the donor star will stay approximately within the Roche lobe. Mass transfer is maintained by the expansion of the donor star, or the contraction of the Roche lobe from the rearrangement of mass and angular momentum in the binary system.

RLOF influences the evolution of the donor star by the decrease in mass. The evolution of the companion star is affected too if some or all of the mass lost by the donor is accreted. This is particularly true if some of the accreted (hydrogen-rich) matter makes its way to the core through internal mixing, where it will thus lead to replenishment of hydrogen, a process known as rejuvenation (see e.g. Vanbeveren & De Loore 1994).

Orbits of close binaries are affected by angular momentum loss (AML) from gravitational wave emission (e.g. Peters 1964), possibly magnetic braking (Verbunt & Zwaan 1981; Knigge et al. 2011) and tidal interaction. Magnetic braking extracts angular momentum from a rotating star by a stellar wind that is magnetically coupled to the star. If the star is in corotation with the orbit, angular momentum is essentially also removed from the binary orbit. Tidal interaction plays a crucial role in circularising binaries and will strive to synchronise the rotational period of each star with the orbital period. While it is known that tidal effects will eventually achieve tidal locking of both components, the strength of tidal effects is still subject to debate (see e.g. Zahn 1977; Hut 1981).

2.1. Stable mass transfer

In the case of conservative RLOF the variation in the orbital separation a during the mass transfer phase is dictated solely by the masses. If the gainer star accretes mass non-conservatively, there is a loss of matter and angular momentum from the system. We define the accretion efficiency: (1)where M_d is the mass of the donor star and M_a is the mass of the accreting companion. If β < 1, it is also necessary to make an assumption about how much angular momentum J is carried away with it. We define this with a parameter η such that: (2)Several prescriptions for AML exist (Appendix B.5) and the amount of angular momentum that is lost from the system due to mass loss has a strong influence on the evolution of the binary.

Matter and angular momentum can also be lost through stellar winds. As these are usually assumed to be spherically symmetric, they will extract the specific orbital angular momentum of the donor star, and result in an increase in the orbital period. If, however, the wind is allowed to interact with the orbit of the binary, the result is entirely dependent on this interaction.

2.2. Unstable mass transfer

During unstable mass transfer, the envelope of the donor star engulfs the companion star. Therefore this phase is often called the common envelope (CE) phase (Paczynski 1976). A merger of the companion and the core of the donor star can be avoided, if the gaseous envelope surrounding them is expelled e.g. by viscous friction that heats the envelope. Because of the loss of significant amounts of mass and angular momentum the CE-phase can have a very strong effect on the binary orbit. In particular it plays an essential role in the formation of short period systems containing at least one compact object. Despite this, the phenomenon is not yet well understood, see Ivanova et al. (2013) for an overview.

There are several formalisms available to treat the orbital evolution during CE-evolution. The most popular ones are the α-formalism (Tutukov & Yungelson 1979) and the γ-formalism (Nelemans et al. 2000). The first considers the energy budget of the initial and final configuration, while the latter is based on the angular momentum balance. Both prescriptions include a parameter after which they are named, which determines the efficiency to remove the envelope. Because such an unstable mass transfer phase occurs on a short timescale, it is often assumed that the gainer does not have the time to gain an appreciable amount of mass during a CE-phase.

The α-parameter describes the efficiency of which orbital energy is consumed to unbind the CE according to: (3)where E_orb is the orbital energy, E_gr is the binding energy of the envelope and α_ce is the efficiency of the energy conversion. The subscript i and f represent the parameter before and after the CE-phase respectively. Several prescriptions for the quantities E_orb,i and E_gr have been proposed (Webbink 1984; Iben & Livio 1993; Hurley et al. 2002) resulting in de facto different α-formalisms. We assume E_orb,i and E_gr as given in the α-formalism of Webbink (1984), such that (4)and (5)where R is the radius of the donor star, M_d,env is the envelope mass of the donor and λ_ce depends on the structure of the donor (de Kool et al. 1987; Dewi & Tauris 2000; Xu & Li 2010; Loveridge et al. 2011).

In the case of mass transfer between two giants with loosely bound envelopes, both envelopes can be lost simultaneously. This process is considered by binary_c/nucsyn, SeBa and StarTrack. The envelopes are expelled according to (6)analogous to Eq. (3), where E_gr,d1 and E_gr,d2 represents the binding energy of the envelope of the two donor stars. This mechanism is termed a double CE-phase (Brown 1995).

3. Binary population synthesis codes

In this paper we compare the results of the simulations of four different BPS codes. These codes have been developed throughout the years with different scientific aims and philosophies, which has resulted in different numerical treatments and assumptions to describe binary evolution. An overview of the methods that are inherent to and the typical assumptions in the four BPS codes can be found in Appendices B and C. Below a short description is given of each code in alphabetical order.

3.1. binary_c/nucsyn

binary_c/nucsyn (binary_c for future reference) is a rapid single star and BPS code with binary evolution based on Hurley et al. (2000, 2002). Updates and relevant additions are continuously made (Izzard et al. 2004, 2006, 2009; Claeys et al. 2013) to improve the code and to compare the effects of different prescriptions for ill-constrained physical processes. The most recent updates (Claeys et al. 2013) that are relevant for this paper are a new formulation to determine the mass transfer rate, the accretion efficiency of WDs and the stability criteria for helium star donors and accreting WDs. The code uses analytical formulae based on detailed single star tracks at different metallicities (based on Pols et al. 1998; Karakas et al. 2002), with integration of different binary features (based on BSE, Hurley et al. 2002). In addition, the code includes nucleosynthesis to follow the chemical evolution of binary systems and their output to the environment (Izzard et al. 2004, 2006, 2009).

The code is used for different purposes, from the evolution of low-mass stars to high-mass stars. This includes the study of carbon- or nitrogen-enhanced metal-poor stars (CEMP/NEMP-stars, Izzard et al. 2009; Pols et al. 2012; Abate et al. 2013), the evolution of Barium stars (Bonačić Marinović et al. 2006; Izzard et al. 2010), progenitor studies of SNe Ia (Claeys et al. 2013), the study of rotation of massive stars (de Mink et al. 2013) and recently the evolution of triple systems (Hamers et al. 2013). Although the code has different purposes, the main strength of the code is the combination of a binary evolution code with nucleosynthesis which enables the study of not only the binary effects on populations, but also the chemical evolution of populations and its output to the environment.

3.2. The Brussels code

The Brussels binary evolution population number synthesis code has been under development for the better part of two decades, primarily to study the influence of binary star evolution on the chemical evolution of galaxies. A thorough review of the Brussels PNS code is given by De Donder & Vanbeveren (2004).

The population code uses actual binary evolution calculations (not analytical formulae) performed with the Paczyński-based Brussels binary evolution code, developed over more than three decades at the Astrophysical Institute of the Vrije Universiteit Brussel. An important feature is that the effects of accretion on the further evolution of the secondary star are taken into account. The population code interpolates between the results of several thousands of actual binary evolution models, calculated under the assumption of the “snowfall model” by Neo et al. (1977) in the case of direct impact, and assuming accretion induced full mixing (see Vanbeveren & De Loore 1994) if accretion occurs through a disk. The actual evolution models have been published by Vanbeveren et al. (1998). The research done with the Brussels code mainly focuses on the chemical enrichment of galaxies caused by intermediate mass and massive binaries. Therefore the interpolations contained in the population code do not allow for the detailed evolution of stars with initial masses below 3 M_⊙.

In recent years, the code was mainly used to study the progenitors of Type Ia supernovae (Mennekens et al. 2010, 2012), the contribution of binaries to the chemical evolution of globular clusters (Vanbeveren et al. 2012) and the influence of merging massive close binaries on Type II supernova progenitors (Vanbeveren et al. 2013).

3.3. SeBa

SeBa is a fast BPS code that is originally developed by Portegies Zwart & Verbunt (1996) with substantial updates from Nelemans et al. (2001a), Toonen et al. (2012) and Toonen & Nelemans (2013). Recent updates include the metallicity dependent single stellar evolution tracks of Hurley et al. (2000) for non-degenerate stars, updated wind mass loss prescriptions and improved prescriptions for hydrogen and helium accretion, and the stability of mass transfer.

The philosophy of SeBa is to not a priori define evolution of the binary, but rather to determine this at runtime depending on the parameters of the stellar system. When more sophisticated models become available of processes that influence stellar evolution, these can be included, and the effect can be studied without altering the formalism of binary interactions. An example of this is the stability criterion of mass transfer and the mass accretion efficiency.

SeBa has been used to study a large range of stellar populations: high mass binaries (Portegies Zwart & Verbunt 1996), double neutron stars (Portegies Zwart & Yungelson 1998), gravitational wave sources (Portegies Zwart & Spreeuw 1996; Nelemans et al. 2001c), double white dwarfs (Nelemans et al. 2001a), AM CVn systems (Nelemans et al. 2001b), sdB stars (Nelemans 2010), SNIa progenitors (Toonen et al. 2012; Bours et al. 2013), post-CE binaries (Toonen & Nelemans 2013) and ultracompact X-ray binaries (van Haaften et al. 2013).

As part of the software package Starlab, it has been used to simulate the evolution of dense stellar systems (Portegies Zwart et al. 2001, 2004). Recently, SeBa is incorporated in the Astrophysics Multipurpose Software Environment, or AMUSE. This is a component library with a homogeneous interface structure, and can be downloaded for free at amusecode.org (Portegies Zwart et al. 2009).

3.4. StarTrack

StarTrack is a Monte Carlo-based single and binary star rapid evolution code. Stars are evolved at a given metallicity (range: Z = 0.0001 − 0.03) by adopting analytical fitting formulae from evolutionary tracks of detailed single stellar models (Hurley et al. 2000), and modified over the years in order to incorporate the most important physics for binary evolution. The orbital parameters (separation, eccentricity and stellar spins) a,e,ω₁ and ω₂ are solved numerically as the system evolves, and re-distribution of angular momentum determines how the orbit behaves. As physical insights regarding various aspects of stellar and binary evolution become available in the literature, new input physics can be implemented into the code, and thus the code is continuously being updated and improved.

The StarTrack code was originally used to predict physical properties of compact objects such and single and double black holes and neutron stars, as well as gamma-ray bursts and compact object mergers in context of gravitational wave detection with LIGO (Belczynski et al. 2002a,b; Abbott et al. 2004). In more recent years, studies with the code have grown to include compact binaries in globular clusters (Ivanova et al. 2005), X-ray binary populations (Belczynski et al. 2004; Ruiter et al. 2006), sources of gravitational wave radiation for ground-based and space-based gravitational wave detectors (Ruiter et al. 2009a, 2010; Belczynski et al. 2010a,b), gamma-ray bursts (Belczynski et al. 2007; O’Shaughnessy et al. 2008; Belczynski et al. 2008b), Type Ia supernovae progenitors (Belczynski et al. 2005; Ruiter et al. 2009b, 2011, 2013) and core-collapse supernova explosion mechanisms (Belczynski et al. 2012). The most comprehensive description of the code to date can be found in Belczynski et al. (2008a), with some updates described in Ruiter et al. (2009b; SNe Ia), Belczynski et al. (2010c; stellar winds), and Dominik et al. (2012; wind mass-loss rates, CE).

4. Method

To examine the inherent differences between four BPS codes, we compare the results of a simulation made by these codes in which the assumptions are equalised as far as possible (Sect. 4.1). We consider two populations of binaries:

• single WDs with a non-degenerate companion (hydrogen-rich or helium-rich star) (SWDs);

• double WD systems (DWDs).

Of both populations we investigate the initial distributions and the distributions at the moment that the SWD or DWD system forms. We establish the similarities between the results of the different BPS codes. If we notice differences between the results, we analyse these in greater detail by comparing e.g. the evolutionary paths or individual systems. A more detailed comparison of the populations of the most important evolutionary paths is given in Appendix A.

In the simulation, we assume an initial primary mass M_1,zams between M_1,zams,min = 0.8 M_⊙ and M_1,zams,max = 10 M_⊙, an initial mass ratio q_zams = M_2,zams/M_1,zams between q_zams,min = 0.1 M_⊙/ M_1,zams and q_zams,max = 1 and an initial semi-major axis a_zams between a_zams,min = 5 R_⊙ and a_zams,max = 10⁴R_⊙ (Table 1). Furthermore we assume an initial eccentricity e_zams of zero. We consider SWDs and DWDs that are formed within a Hubble time, more specifically 13.7 Gyr. The initial distribution of the primary masses follows Kroupa et al. (1993), the initial mass ratio distribution is flat¹, and the initial distribution of the semi-major axis is flat in a logarithmic scale.

Not every BPS research group focuses on the full range of stellar masses. Consequently in their codes there are no (valid) prescriptions available for all stellar masses. The research group that uses the Brussels code, mainly focuses on the chemical enrichment of galaxies and therefore is not interested in the evolution of stars with a mass lower than 3 M_⊙ (Sect. 3.2). Consequently, in order to make the comparison with the results of the Brussels code we only compare with a subset of the SWD and DWD populations. We define this subset as the “intermediate mass range”, while the entire populations is considered as the “full mass range”. The “intermediate mass range” is defined in the two populations as follows:

• for the SWD population we only consider WDs originating from initial primary masses higher than 3 M_⊙;

• for the DWD population we only consider WDs originating from initial primary and secondary masses both higher than 3 M_⊙.

In addition, we refer to the “low mass range” or “low mass primaries” which encompasses the systems with an initial primary mass lower than 3 M_⊙.

BPS codes are ideal to investigate the effect of different assumptions on populations, since a different assumption can cause a shift in e.g. the mass or separation of the population under investigation. We do not have to agree on the exact evolution of individual systems. As long as the shift is small the characteristics of the population do not change. Keeping this in mind when comparing the results of the different BPS codes we define them to agree when simulated populations (of similar evolutionary paths) are recovered at the same regions in the mass and separation space.

4.1. Assumptions for this project

In order to compare the codes we make the most simple assumptions. These are not necessarily believed to be realistic, but are taken to make the comparison feasible. The assumptions for this project are discussed below and shown in Table 1. The typical assumptions taken by the authors in the corresponding BPS codes in their previous research projects are summarised in Table C.1 in Appendix C. For simplicity and brevity, we do not study the effect of these assumptions on the characteristics of SWD and DWD populations in this project.

• Mass transfer is assumed to be conservative (β = 1) during stable RLOF towards all types of objects. We emphasise that this is not a realistic assumption, especially in the case of a WD accretor. During the CE-phase no material is assumed to be accreted by the companion star (β = 0). In the Brussels code a constant accretion efficiency of a WD-accretor cannot be implemented and therefore for this study mass transfer to all compact objects is assumed to be unstable and evolve into a CE-phase in this code.

  	Fig. 1 Orbital separation versus WD mass for all SWDs in the full mass range at the time of SWD formation.
  Open with DEXTER

  	Fig. 2 Orbital separation versus WD mass for all SWDs in the intermediate mass range at the time of SWD formation.
  Open with DEXTER

  	Fig. 3 Secondary mass versus WD mass for all SWDs in the full mass range at the time of SWD formation.
  Open with DEXTER

• As no mass nor angular momentum is lost from RLOF, we do not require an assumption for the specific angular momentum loss of the material. During wind mass loss, we assume the wind matter leaves the system with the specific angular momentum of the donor star. However, this assumption is not possible in the Brussels code (for an overview of the assumptions see Sect. B.5).

• We use the α-prescription of Webbink (1984) to describe the CE-phase (Eqs. (3)–(5)). We assume that the parameters α_ce and λ_ce are equal to one, mainly for simplicity, but also because the prevalence of this choice in the literature allows for comparison between this and other studies.

• We assume that matter lost through winds cannot be accreted by the companion star.

• Due to the diversity of the prescriptions for magnetic braking and tides, we do not consider these effects and they are turned off for this paper. However, in StarTrack spin-orbit coupling is still taken into account, as it is firmly integrated with the binary evolution equations.

4.2. Normalisation

When calculating birthrates of evolutionary channels, the simulation has to be normalised to an entire stellar population (Table 1). For this work the initial distribution and ranges of M_1,zams, q_zams and a_zams are as discussed in Sect. 4 with the exception of the initial primary masses of a stellar population to vary between 0.1 and 100 M_⊙, and the semi-major axis between 5 and 10⁶R_⊙. We assume a binary fraction of 100%.

If the star formation rate S in M_⊙  yr^-1 is independent of time, the birthrate of a specific binary type X (e.g. systems evolving through a specific evolutionary channel) is given by: (7)with φ(X) the total number of systems of binary type X in the simulation, and M_tot the total mass of all stellar systems in the entire stellar population. More specifically, (8)with x = x(M_1,zams,q,a) equals 1 for binary systems of type X, and zero otherwise and Ψ is the initial distribution function of M_1,zams, q_zams and a_zams. Note that in this project we assume that the initial distribution for M_1,zams, q_zams and a_zams are independent (Table 1), such that Ψ is separable: (9)The total mass of all stellar systems assuming a 100% binary fraction is: (10)where M_t,zams = M_1,zams + M_2,zams.

For this project a constant star formation rate of 1 M_⊙  yr^-1 is assumed. This simple star formation rate is chosen to make the comparison with other codes easier.

5. Comparison

5.1. Single white dwarf systems

Systems containing a WD and a non-degenerate companion have typically undergone a one-directional mass transfer event i.e. one star has lost mass and possibly the other gained mass. The mass transfer event may consist of one or two episodes, either of which may have been stable or unstable. The characteristics of the population of SWD systems show the imprint of the mass transfer episodes. Figures 1 and 2 show the orbital separation a_swd as a function of primary mass M_1,swd at the moment of WD formation for the full and intermediate mass range respectively. Likewise Figs. 3 and 4 show the secondary mass M_2,swd as a function of primary mass at WD formation for the full and intermediate mass range. These figures show that in general the codes find very similar SWD systems.

In more detail, at large separations (a_swd ≳ 500 R_⊙ for the full mass range, and a_swd ≳ 2000 R_⊙ for the intermediate mass range) all codes find systems in which the stars do not interact. The population of SWDs with WD masses in the low mass range is very comparable in orbital separation, primary and secondary mass between the codes binary_c, SeBa and StarTrack. Intermediate mass systems can be divided in two groups, either in separation and/or in secondary mass. According to all codes, intermediate mass systems that undergo a CE-phase (for the first mass transfer episode) are compact with a_swd ≲ 200 R_⊙ and have secondary masses up to 10 M_⊙. Furthermore, the codes agree that in the intermediate mass range, systems for which the first phase of mass transfer is stable are in general more compact than non-interacting systems and less compact than the systems undergoing a CE-phase. The secondary mass is between 3 and 18 M_⊙ as it accretes conservatively during stable mass transfer.

The ZAMS configurations for progenitors of SWDs are shown in Figs. 5 and 6 with the separation a_zams versus primary mass M_1,zams. There is a general agreement between the codes about which progenitor systems lead to a SWD system and which systems do not. According to all codes, compact progenitor systems (a_zams ≲ 400 R_⊙ for the intermediate mass range, while a_zams ≲ 30 R_⊙ for the low mass range) undergo stable mass transfer for the initial mass transfer episode. Furthermore the codes agree that for most progenitor systems with orbital separations in the range a_zams ≈ (0.1−3)    ×    10³R_⊙ the first phase of mass transfer is unstable. Systems with orbital separations that lie between the ranges described above lead to a merging event, thus no SWD system is formed. Progenitor systems with a_zams ≳ 700 R_⊙ for the intermediate mass range (a_zams ≳ 250 R_⊙ for the low mass range) are too wide for the primary star to fill its Roche lobe.

Overall the simulations of the four codes show a good agreement on the characteristics of the population of SWDs in orbital parameters and birthrates (Table 2), however, differences can be noted. The most important causes are the relation between the initial and WD mass, the stability of mass transfer and the modelling of the stable mass transfer phase. The initial-final mass (MiMf)-relation of single stars (Fig. 7) is very similar between binary_c, SeBa and StarTrack, but different than the one from the Brussels code due to different single star prescriptions that are used in the latter code (see also Appendix A.1.1 for a discussion). The effect on the population of SWD progenitors can be seen in Fig. 6 in the maximum mass of the primary stars which is extended from about 8 M_⊙ in binary_c, SeBa and StarTrack to about 10 M_⊙ in the Brussels code. For binary stars the relation between WD mass and the initial mass is hereafter called the initial-WD mass (MiMwd)-relation (Appendix A.1.2). Differences in the MiMwd-relation lead to an increase of systems at small WD masses ≲0.64 M_⊙ in Fig. 2 in the Brussels code compared to the other codes. The gap in WD masses between 0.7−0.9 M_⊙ in the Brussels data in Fig. 2 is a result of a discontinuity in the MiMwd-relation between the WD masses of primaries that fill their Roche a second time, and those that do not. In the other codes, the primary WD masses of binaries that evolve through these two evolutionary channels are overlapping. Differences in the stability criteria of mass transfer can be seen in Figs. 2 and 4, where the StarTrack code shows a decrease of systems that underwent stable mass transfer (Appendix A.1.3). Mass transfer is modelled differently in the codes (Sect. B) leading to an extension to small separations in the Brussels data compared to the other codes (Fig. 6), and an increase in systems that underwent stable mass transfer at a_zams ≈ 10 R_⊙ for M_1,zams ≳ 4 M_⊙ in Fig. 4 (Appendix A.1.5).

For a more detailed comparison of the SWD population in the full and intermediate mass range, see Appendix A.1.

5.2. Double white dwarfs

In this section we compare and discuss the population of DWDs as predicted by binary_c, the Brussels code, SeBa and StarTrack. Prior to the formation of a second degenerate component, DWDs undergo the evolution as described in the previous section. Subsequently, they undergo a second intrusive (series of) event(s) at the time the secondary fills its Roche lobe. As a consequence the processes that influence the evolution of SWDs influence the DWD population as well. Here we will point out the evolutionary processes that are specifically important for DWDs.

The population of DWDs at DWD formation is shown in Figs. 8–11 where orbital separation and secondary mass respectively is shown as a function of primary mass for the full and intermediate mass range. The ZAMS progenitors of these systems are shown in Figs. 12 and 13 for the full and intermediate mass range respectively.

In the full mass range, the population of DWDs is comparable between binary_c, SeBa and StarTrack with WD masses of M_1,dwd ≈ 0.2−1.4 M_⊙ and M_2,dwd ≈ 0.1−1.4 M_⊙. At large separations (0.1−5)    ×    10⁴R_⊙ the codes find systems which are formed without any interaction, see Fig. 8. This figure also shows a population of interacting systems at lower separations, where the majority has separations of a ≈ 0.1 − 10 R_⊙. Furthermore there is a good agreement on which progenitors lead to a DWD system and which do not. Figure 12 shows several subpopulations of DWD progenitors with comparable binary parameters for binary_c, SeBa and StarTrack; a group of non-interacting systems (at a_dwd ≳ 5 × 10²R_⊙), a group of systems for which the first phase of mass transfer is stable (at a_dwd ≲ 25 R_⊙ for low mass primaries and a_dwd ≲ 2.5 × 10²R_⊙ for the full mass range), and a group of systems at intermediate separations that predominantly undergoes a CE-phase for the first phase of mass transfer.

Effects that play a role when comparing DWDs in the full mass range are the stability of mass transfer and differences in the interpretation of the double CE-phase in which both stars lose their envelopes (Eq. (5)). The most pronounced effect of the differences in the stability of mass transfer is the decrease of systems that underwent stable mass transfer in the StarTrack data compared to binary_c and SeBa. This can be seen in Fig. 12 in the lack of systems at M_1,zams > 3 M_⊙ and a_zams < 2.5 × 10²R_⊙ in the StarTrack data compared to binary_c and SeBa, and in Fig. 10 in the lack of systems with M_2,dwd > M_1,dwd. Furthermore differences in the stability of mass transfer lead to an increase in systems at a_dwd ≈ 10−50 R_⊙ and M_1,dwd ≈ 0.4 − 0.47 M_⊙ according to SeBa and StarTrack. Differences in the modelling of the double CE-phase result in larger separations at DWD formation and less mergers in StarTrack compared to binary_c and SeBa (Fig. 9 and Appendix A.2.2). At the same time, the initial separations of systems evolving through a double CE-phase can be smaller in the StarTrack data compared to binary_c and SeBa (a_dwd ≈ 25−100 R_⊙, see Fig. 12).

In the intermediate mass range at DWD formation, two groups of systems can be distinguished in all simulations (Fig. 9). Similar to the full mass range, there is one group of non-interacting systems at separations higher than 6    ×    10³R_⊙ and a group of interacting systems with separations ≲20 R_⊙. However, the distribution of systems in the latter range varies between the codes. Most DWD systems have primary and secondary WD masses above 0.6 M_⊙ in all the codes. The progenitors in the intermediate mass range show the same division in separation in three groups as the progenitors in the full mass range. DWD progenitors with separations a_zams ≲ 3    ×    10²R_⊙ undergo a stable first phase of mass transfer. The components of DWD progenitors with a_zams ≳ 1.5 × 10³R_⊙ do not interact. At intermediate separations the first phase of mass transfer is predominantly a CE-phase.

Comparing the Brussels code with binary_c and SeBa (differences with StarTrack have the same origin as discussed in previous paragraphs), the most important causes for differences in the DWD population in the intermediate mass range are the MiMf-relation, the MiMwd-relation, the modelling of the stable mass transfer phase and the survival of mass transfer. The effect of the first three causes on the DWD population is similar to the effect on the SWD population. Firstly, the differences in the MiMf-relation can be seen in the progenitor population of non-interacting binaries in Fig. 13 as an extension to higher primary masses in the Brussels data (8−10 M_⊙, see also Appendix A.2.1). Secondly, differences in the MiMwd-relation can be seen in Fig. 9 as an extension to lower primary WD masses M_1,dwd ≲ 0.64 M_⊙ and the discontinuity in primary WD masses around 0.7−0.9 M_⊙ (Appendices A.2.2 and A.2.3). The MiMwd-relation also effects the orbital separation distribution at DWD formation and results in a higher maximum separation in the Brussels code compared to binary_c and SeBa. Finally, due to the method of modelling of mass transfer there is a disagreement between the codes regarding which systems survive mass transfer, see Fig. 13 at a_dwd ≲ 20 R_⊙ (Appendix A.2.3). The survival of mass transfer is more important for the DWD population than for the SWD population, as the average orbital separation of DWDs is lower (Sect. 5.1 and also Appendices A.2.2 and A.2.3). As the formation of DWDs involves more phases of mass transfer than for SWDs, the differences in the SWD population carry through and are larger in the DWD population. The DWD population in the full and intermediate mass range are discussed in more detail in Appendix A.2.

6. Overview of critical assumptions in BPS studies

In the previous section we compared simulations from four different BPS codes and investigated the causes for the differences. The causes that we found are not numerical effects, but are inherent to the codes. In this section we list and discuss the underlying physical principles of the differences described in Sect. 5. The implementations of these principles in each code are described in Appendix B.

• Initial-WD mass-relation. For single stars or non-interacting stars, the initial-final mass relation for WDs (Fig. 7) is determined by the trade off between the growth of the core and how much mass is lost in stellar winds. The amount of mass a low or intermediate mass star loses in a stellar wind is small on the MS, but significant in later stages of its evolution. The amount of mass that is lost in the wind and in the planetary nebula influences the orbit directly, and indirectly through angular momentum loss (Sects. B.4 and B.5). The WD mass of primary stars is further affected by the mass transfer event, the moment and the timescale of the removal of the envelope mass. If the primary star becomes a hydrogen-poor helium burning star before turning into WD, the MiMwd-relation is influenced by helium star evolution. Of importance are the core mass growth versus the mass loss from helium stars and a possible second phase of mass transfer. A related issue, of particular importance for supernova Type Ia rates, concerns the composition of WDs; what is the range of initial masses for carbon-oxygen WDs or other types of WDs?

• The stability of mass transfer. For which systems does mass transfer occur in a stable manner and for which systems is it unstable? As BPS codes do not solve the stellar structure equations, and cannot model stars that are not in hydrostatic or thermal equilibrium, BPS codes rely on parametrisations or interpolations to determine the stability of mass transfer. Theoretical stability criteria for polytropes exist (Hjellming & Webbink 1987), but are lacking for most real stars (but see de Mink et al. 2007; Ge et al. 2010, 2013, for MS stars). The critical mass ratio for stable mass transfer with hydrogen shell-burning donors differs between the codes from q ≳ 0.2 in the Brussels code to q ≳ 0.6 in StarTrack. A difference for low mass stars between binary_c, SeBa and StarTrack arises from the uncertainty of the mass transfer stability of donors with shallow convective envelopes. In a recent paper, Woods et al. (2012) show that mass transfer between a hydrogen shell-burning donor (M_1,zams = 1 − 1.3 M_⊙) and a main-sequence star can be stable when non-conservative. The effect on the orbit is a modest widening.

• Survival of mass transfer. For which systems does mass transfer lead to a merger and which system survive the mass transfer phase, in particular when mass transfer is unstable? Different assumptions for the properties (e.g. radii) of stripped stars lead to differences in the results of the four BPS codes (e.g. channels II and III in Appendix A.2). For donor stars in which the removal of the envelope due to mass transfer leads to an end in nuclear (shell) burning and a WD is formed directly, it is unclear how much the core is bloated just after mass transfer ceases compared to a zero-temperature WD (Hurley 2000). For donor stars that are stripped of their hydrogen envelopes due to mass transfer, but helium burning continues, it is unclear how fast the transition takes place from an exposed core to an (evolved) helium star (channel 2b in Appendix A.1).

• Stable mass transfer. Modelling of the stable mass transfer phase in great detail is not possible in BPS codes, as for the stability of mass transfer. Therefore BPS codes rely on simplified methods to simulate stable mass transfer events. The evolution of the mass transfer rate during the mass transfer phase can have a strong effect on the resulting binary. However, in the current set-up of this project that assumes conservative mass transfer, the importance is greatly reduced. The mass transfer rates are only important when the timescale of other effects (e.g. wind mass loss or nuclear evolution) become comparable to the mass transfer timescale (channel 3b in Appendix A.1). A result of the approach is that mergers are less likely to happen in the Brussels code compared to the other codes (e.g. channel 5 in Appendix A.1). The approach of binary_c, SeBa and StarTrack is to follow the mass transfer phase in time, with approximations of the mass transfer rate. In the Brussels code, the mass transfer phases are not followed in detail. Instead it only considers the initial and final situation from interpolations of a grid of detailed calculations. Furthermore, it is important to better understand which contact systems lead to a merger and which do not. From observations, many Algol systems are found which have undergone and survived a phase of shallow contact.

• The evolution of helium stars. A large fraction of interacting systems go through a phase in which one of the stars is a helium star, for SWDs roughly 15% in the full mass range and roughly 50% in the intermediate mass range. These objects are not well studied and there remain several uncertainties, e.g. mass transfer stability. Also the mass transfer rate is important, in particular for evolved helium stars whose evolutionary and wind loss timescales can become comparable to the mass transfer timescales. Therefore small differences in the mass transfer rate can lead to large differences in the resulting WD. This is especially important for massive WDs, e.g. SNIa rates.

The influence of the parameters that were equalised in this project has not been studied here, neither qualitatively nor quantitatively. These parameters will lead to a larger diversity in the simulated populations as different groups make different assumptions in their codes (Appendix C) and these should be taken into account when interpreting BPS results. These assumptions are:

• the CE-prescription and efficiency;

• accretion efficiency;

• angular momentum loss during RLOF;

• tidal effects;

• magnetic braking;

• the initial distributions of primary mass, mass ratio and orbital separation.

Despite the significance of these phenomena to binary evolution and the efforts of the community to understand and quantify them, there remain questions about these phenomena. Several prescriptions exist for these phenomena and the effect on a binary population can be severe. For example regarding the CE-phase, it is unclear how efficient orbital energy can be used to expel the envelope and if other sources of energy can be used, such as recombination, rotational, tidal or magnetic energy (Iben & Livio 1993; Han et al. 1995; Politano & Weiler 2007; Webbink 2008; Zorotovic et al. 2010; De Marco et al. 2011; Zorotovic et al. 2011; Davis et al. 2012; Ivanova et al. 2013). Also, predictions for the efficiency of mass accretion onto WDs vary strongly and the SNIa rate is severely affected by this uncertainty (Bours et al. 2013). Furthermore, the adopted mode of angular momentum loss has a strong effect on the evolution of the orbit (Fig. C.1 and Appendix C.2). It is also not clear how the different prescriptions for tidal evolution affect the populations. However, in Appendix A.1.3 we find that spin-orbit coupling (assuming orbits are continuously synchronised), only has a small effect on the final separation of the SWD systems. The effect of different initial distributions (see e.g. Duquennoy & Mayor 1991; Kouwenhoven et al. 2007) of binary parameters can be severe with respect to the birthrate of a stellar population (see e.g. Eggleton et al. 1989; de Kool & Ritter 1993; Davis et al. 2010; Claeys et al. 2013). Furthermore, the importance of a certain channel is affected by the boundaries of the distribution through the normalisation of the simulation.

7. Conclusion

In this paper we studied and compared four BPS codes. The codes involved are the binary_c code (Izzard et al. 2004, 2006, 2009; Claeys et al. 2013), the Brussels code (De Donder & Vanbeveren 2004; Mennekens et al. 2010, 2012), SeBa (Portegies Zwart & Verbunt 1996; Nelemans et al. 2001a; Toonen et al. 2012; Toonen & Nelemans 2013) and StarTrack (Belczynski et al. 2002a, 2008a; Ruiter et al. 2009b; Belczynski et al. 2010c). We focused on low and intermediate mass binaries that evolve into SWD systems (containing a WD and a non-degenerate companion) and double white dwarf systems. These populations are interesting for e.g. post-CE binaries, cataclysmic variables, single degenerate as well as double degenerate supernova Type Ia progenitors. For this project input assumptions in the BPS codes were equalised as far as the codes permit. This was done to simplify the complex problem of comparing BPS codes that are based on many (often different) assumptions. In this manner inherent differences between and numerical effects within the codes were investigated.

Regarding the SWD population, there is a general agreement on what initial parameters of M_1,zams, M_2,zams and a_zams lead to SWD binaries and which parameters do not lead to SWDs. When the SWD system is formed, there is an agreement on the orbital separation range for those systems having undergone stable or unstable mass transfer. Furthermore there is a general agreement on the stellar masses after a phase of stable or unstable mass transfer and between the populations of the most common evolutionary channels.

Regarding the DWD population, there is an agreement on which primordial binaries lead to DWD systems through stable and unstable mass transfer respectively, and a rough agreement on the characteristics (M_1,dwd, M_2,dwd and a_dwd) of the DWD population itself. DWD systems go through more phases of evolution than SWD systems and therefore the uncertainty in their evolution builds up after each mass transfer phase. The WDs are formed with comparable masses, but at different separations. The most important evolutionary paths leading to DWDs are similar between the BPS codes.

We found that differences between the simulated populations are not due to numerical differences, but due to different inherent assumptions. The most important ones that lead to differences are the MiMf-relations (of single stars), the MiMwd-relation (of binary stars), the stability of mass transfer, the modelling of the mass transfer rate and the modelling of helium star evolution. Different assumptions between the codes are made for these topics as theory is poorly understood and sometimes poorly studied. Further research into these topics is necessary to eliminate the differences between BPS codes e.g. with a detailed (binary) stellar evolution code.

In addition some assumptions that affect the results of the codes were equalised for the comparison. These are the initial binary distributions, the CE-prescription and efficiency, the accretion efficiency, angular momentum loss during RLOF, tidal effects, magnetic braking and wind accretion. We leave the study of their effects on stellar populations for another paper.

In Sect. 3 a short description is given of each code. In Appendices B and C, a more detailed overview is given of the typical assumptions of each code outside the current project. These should be taken into account when interpreting results from the BPS codes. Furthermore, we recommend using these sections as a guideline when deciding which code or results to use for which project. Finally we would like to encourage other groups involved in BPS simulations, to do the same test as described in this paper and compare the results with the figures given in this paper. More detailed figures are available on request and on the website www.astro.ru.nl/~silviato/popcorn

Concluding, we found that when the input assumptions are equalised as far as possible within the codes, we find very similar populations and birthrates. Differences are caused by different assumptions for the physics of binary evolution, not by numerical effects. So although the four BPS codes use very different ways to simulate the evolution of these systems, the codes give similar and consistent results and are adequate for studying populations of low- and intermediate-mass stars.

---

Acknowledgments

The authors would like to thank G. Nelemans, M. van der Sluys, O. Pols, D. Vanbeveren, W. van Rensbergen, & R. Voss for very helpful discussions on this project. S.T. thanks S. Portegies Zwart, J.S.W.C. thanks R.G. Izzard and A.J.R. thanks I. Seitenzahl and K. Belczynski for thoughtful comments and general discussion. The authors are very grateful to those involved in the Lorentz Center workshops “Stellar mergers” (2009) and “Observational signatures of Type Ia supernova progenitors” (2010), who provided the inspiration to carry out this long-awaited work. Furthermore the authors are grateful for the contributions to an early version of this project in 2002−2003 from K. Belczynski, C. Fryer, J. Hurley, V. Kalogera, G. Nelemans, S. Portegies Zwart and C. Tout, and in 2009−2010 from A. Bogomazov, P. Kiel, B. Wang and L. Yungelson. The work of S.T. and J.S.W.C. was supported by the Netherlands Research School for Astronomy (NOVA).

References

Online material

Appendix A: Detailed comparison

Appendix A.1: Single white dwarf systems

In the next sections, we make a more detailed comparison between the simulated populations of SWDs of the four codes. We distinguish between the most commonly followed evolutionary paths with birthrates larger than 1.0    ×    10^-3   yr^-1 (Table A.1). We describe each evolutionary path, the similarities and differences, and investigate the origin of these differences. Specific examples are given and discussed for the most common paths. Abbreviations of stellar types are shown in Table A.2. Paragraphs explaining the evolutionary path, an example evolution and the comparison of the simulated populations for each evolutionary channel are indicated with Evolutionary path, Example and Population, respectively. For some channels, causes for differences between the populations are discussed separately in paragraphs that are indicated by Effects. Masses and orbital separations according to each code are given in vector form [c₁,c₂,c₃,c₄] where c₁ represents the value according to the binary_c code, c₂ according to the Brussels code, c₃ according to SeBa, and c₄ according to StarTrack. The examples are given to illustrate the evolutionary path and relevant physical processes. However, note that when comparing different BPS codes, achieving similar results for specific binary populations is more desirable and important than achieving a perfect match between specific, individual binary systems.

Appendix A.1.1: Channel 1: detached evolution

Evolutionary path. Most SWD binaries are non-interacting binaries where the stars essentially evolve as single stars. Most binary processes that are discussed in Sect. 2 do not play a role in channel 1.

Example. As an example of a system in channel 1, we discuss the evolution of a system that initially contains a 5 M_⊙ and 4 M_⊙ star in an orbit of 10⁴R_⊙ (and e_zams = 0 by assumption). When the primary star becomes a WD its mass is [1.0,0.94,1.0,1.0] M_⊙ in an orbit of [1.8,1.8,1.8,1.8]     ×    10⁴R_⊙. The differences in the resulting SWD system from different BPS codes are small and mainly due to different initial-final mass (MiMf)-relations (Fig. 7). The maximum progenitor mass to form a WD from a single star is [7.6,10,7.9,7.8] M_⊙ and corresponding maximum WD mass of [1.38,1.34,1.38,1.4] M_⊙ according to the four codes. The MiMf-relations of the binary_c code, SeBa and StarTrack are very similar. The similarities are not surprising as these codes are based on the same single stellar tracks and wind prescriptions of Hurley et al. (2000). However, small differences arise in the MiMf-relation as the prescriptions for the stellar wind are not exactly equal. The Brussels code is based on different models of single stars e.g. different stellar winds and a different overshooting prescription (Appendix B). The result is that the core mass of a specific single star is larger according to the Hurley tracks. In other words, the progenitor of a specific SWD is more massive in the Brussels code.

Population. Despite differences for individual systems, the population of non-interacting binaries at WD formation is very similar. The previously mentioned differences in the MiMf-relations are noticeable in the maximum initial primary mass in Figs. A.5 and A.6. The distribution of separations at WD formation (Figs. A.1 and A.2) are very similar between the codes. For the intermediate mass range, the separations at SWD formation are ≳ 4.5    ×    10³R_⊙ for the Brussels code and extend to slightly lower values of ≳ 2.0    ×    10³R_⊙ for binary_c, SeBa, and StarTrack. For the full mass range, the latter three codes agree that the separations can be as low as 5.0    ×    10²R_⊙. The progenitor systems of channel 1 have similar separations of ≳ 3.0    ×    10²R_⊙ for

low mass primaries. For intermediate mass stars binary_c, SeBa and StarTrack find that the initial separation is ≳ 0.7    ×    10³R_⊙ where the Brussels code finds a slightly higher value of ≳ 1.6    ×    10³R_⊙ (Figs. A.5 and A.6). The minimum separation (at ZAMS and WD formation) for a given primary mass depends on whether or not the primary fills its Roche lobe, which in turn depends on the maximum radius for that star according to the particular single star prescriptions that are used. Even though the progenitor populations are not 100% equal, the characteristics of the SWD population and the birthrates (Table A.1) in this channel are in excellent agreement.

Appendix A.1.2: Channel 2: unstable case C

Evolutionary path. One of the most common evolutionary paths of interacting binaries is channel 2, of which an example is shown in Fig. A.7. In this channel, the primary star fills its Roche Lobe when helium is exhausted in its core, so-called case C mass transfer (Lauterborn 1970). As the envelope of the donor star is deeply convective at this stage, generally mass transfer leads to an unstable situation and a CE-phase develops. While the orbital separation shrinks severely, the primary loses its hydrogen envelope. By assumption in this project, the secondary is not affected during the CE-phase. The primary can either directly become a WD or continue burning helium as an evolved helium star as shown in the example of Fig. A.7. If the primary becomes a WD directly, or indirectly but without further interaction, the evolutionary path is called channel 2a. Evolution according to channel 2b occurs if the primary fills its Roche lobe for a second time when it is a helium star. The second phase of mass transfer can be either stable or unstable.

Example. As an example of channel 2a, we discuss the evolution of the binary system in Fig. A.7 with initial parameters M_1,zams = 3.5  M_⊙, M_2,zams = 3  M_⊙ and a_zams = 350  R_⊙ in more detail. The primary star fills the Roche lobe early on the AGB before thermal pulses and superwinds occur. Wind mass loss prior to the CE-phase is small, [4.4,0,4.3,4.9]     ×    10^-2M_⊙. After the CE-phase the orbital separation is reduced to [14,9.1,14,14] R_⊙. In this example the primary continuous burning helium as an evolved helium star of mass [0.78,0.55,0.78,0.78] M_⊙. When the primary exhausts its fuel, it becomes a WD of [0.76,0.51,0.77,0.76] M_⊙ in an orbit of [14,9.1,14,14] R_⊙ with a 3 M_⊙ MS companion. The most important differences, to be seen between the Brussels code and the other codes, arises from the different single star prescriptions that are used. This affects the resulting mass of a WD from a specific primary, and the resulting orbital separation. Note that while the MiMf-relation for single stars depends on the single star prescriptions (i.e. core mass growth and winds), the MiMwd-relation is also affected by the companion mass and separation (which determine when and which kind of mass transfer event takes place), and the single star prescriptions for helium stars. In other words, the MiMwd-relation represents how fast the core grows on one hand, and the envelope is depleted by mass transfer and stellar winds on the other hand.

Population. Despite the differences between individual systems, the different BPS codes agree in which regions of phase space (M_1,swd, M_2,swd, a_swd) in Figs. A.8–A.11 the systems from channel 2 lie. The systems of channel 2 evolve towards small separations, with the majority in the range 0.2 − 150 R_⊙ at WD formation. In addition, the codes agree on the masses of both stars at formation of the SWD system. In the low mass range binary_c, SeBa and StarTrack find and agree that M_1,swd ≈ 0.5 − 0.7 M_⊙ and M_2,swd ≈ 0.1 − 2.7 M_⊙. In the intermediate mass range the different codes find that M_1,swd ≳ 0.64 M_⊙, however, the Brussels code finds primary WD masses down to 0.5 M_⊙ due to differences in MiMwd-relation. For secondary masses the codes find M_2,swd ≈ 0.1 − 7.0 M_⊙. The binary_c, SeBa, and StarTrack codes agree on the initial separation for low mass binaries, which is between (0.6 − 12)    ×    10²R_⊙ (Fig. A.12), M_1,zams ≈ 1.0 − 3.0 M_⊙ and M_2,zams ≈ 0.1 − 3.0 M_⊙. For intermediate mass binaries in channel 2, there is an agreement between all codes that the initial primary masses lie between M_1,zams ≈ 3 − 8.5 M_⊙ and M_2,zams ≈ 0.1 − 7.7 M_⊙. Due to the MiMwd-relation, the maximum initial primary mass extends to slightly higher values for the Brussels code in comparison with the other codes (Fig. A.13). However, for massive primary progenitors e.g. M_1,zams > 9 M_⊙ in the Brussels code, the envelope mass of the donor is large and therefore a merger is more likely to happen in the simulations of the Brussels code compared to those of the other three codes. The initial orbital separation lies between (0.1 − 2.4)    ×    10³R_⊙ (Fig. A.13) according to binary_c, SeBa and StarTrack, however, the range is extended to 3.2    ×    10³R_⊙ in the Brussels code due to the single star prescriptions of stellar radii.

Effects. Comparing channels 2a and 2b separately, the birthrates of SWDs (Table A.1) in the full mass range are close between the codes binary_c, SeBa, and StarTrack. In the intermediate mass range for channel 2a, the birthrates of binary_c, SeBa, and StarTrack are essentially identical, and within a factor of 2.5 lower compared to that of the Brussels code. The larger difference with the Brussels code are caused because this code assumes a priori that a WD is formed without a second interaction, thus there is no entry for the Brussels code in Table A.1 for channel 2b. The birthrates for channel 2b are very similar within a factor of about 1.2 between binary_c, SeBa and StarTrack. Comparing the total birthrate in channel 2 between all codes, the rate of binary_c, SeBa and StarTrack is only lower by about a factor 1.5 compared to the Brussels code, as some systems merge in the second interaction in the simulations of the former codes. Other differences in the simulated populations from this channel are due to the MiMwd-relation as seen in the example, but also due to differences in the criteria for the stability of mass transfer and the prescriptions for the wind mass loss (see below).

The effect of the stellar wind in the example above is negligible, but the effect of wind mass loss becomes more important for systems with more evolved donors. Mass loss from the primary either in the CE-phase or in foregoing wind mass loss episodes affects the maximum orbital separation of the SWD systems directly and through angular momentum loss. In the simulations of the Brussels code, the maximum orbital separations at WD formation are lower (a_swd ≲ 80 R_⊙ compared to ≲ 150 R_⊙ for the main group of systems in binary_c, SeBa and StarTrack), as winds are not taken into account and more mass is removed during the CE-phase in this code. More mass loss during a CE-phase leads to a greater shrinkage of the orbit, where as more wind mass loss with the assumption of specific angular momentum loss from the donor (Jeans-mode, see Eq. (C.1)), leads to an orbital increase.

Another effect arises from the stellar wind in combination with the stability criterion of mass transfer. For systems with high wind mass losses in which the mass ratio has reversed, the first phase of mass transfer can become stable according to binary_c, SeBa, and StarTrack. Systems in which this happen are not included in channel 2, however, the birthrates are low ([1.3, − ,6.5,4.7]     ×    10^-4   yr^-1 in the full mass range and [5.4, − ,10,9.1]     ×    10^-5   yr^-1 in the intermediate mass range). In general, when a AGB star initiates mass transfer, stable mass transfer is more readily realised in SeBa and StarTrack than in binary_c. Therefore the maximum separation of SWDs in channel 2 is highest in the binary_c data (up to 650 R_⊙). However, only about 1% of systems in channel 2 in the binary_c code lie in the region with a separation larger than 70 R_⊙ and a WD mass higher than 0.6 M_⊙.

The stability of mass transfer is another important effect for the population of systems in channel 2b during the second phase of mass transfer. We only compare the binary_c code, SeBa, and StarTrack, as the Brussels code does not consider this evolutionary path. Whether or not the second phase of mass transfer is stable affects the resulting distribution of orbital separations. This effect is shown in Fig. A.8 as an extension to lower separations a_swd ≲ 10 R_⊙ for M_1,swd ≳ 0.8 M_⊙ in the binary_c data due to unstable mass transfer.

There is a difference between StarTrack on one hand, and binary_c and SeBa on the other hand regarding the survival of systems in channel 2b during the first phase of mass transfer. Due to a lack of understanding of the CE-phase, generally BPS codes assume for simplicity that when the stars fit in their consecutive Roche lobes after the CE is removed, the system survives the CE-phase. However, this depends crucially on the evolutionary state of the stars after the CE. For channel 2b in which the primary continues helium burning in a shell as a non-degenerate helium star, the response of the primary to a sudden mass loss in the CE-phase is not well known. The StarTrack code assumes the stripped star immediately becomes an evolved helium star and corresponding radius, while binary_c and SeBa assume the stripped star is in transition from an exposed core to an evolved helium star with a radius that can be a factor of about 1−15 smaller. The uncertainty in the radii of the stripped star mostly affect systems with M_1,zams ≳ 5 M_⊙ at separations ≳450 R_⊙ that merge according to StarTrack, and survive according to binary_c and SeBa.

Included in channel 2 are systems that evolve through a double CE-phase² in which both stars lose their envelope described in Sect. 2 and in Eq. (6). The double CE-mechanism is taken into account by the binary_c, SeBa and StarTrack code. However, there is a difference between StarTrack on one hand, and binary_c and SeBa on the other hand regarding the binding energy of the envelope of the secondary star. In StarTrack the binding energy is calculated according to Eq. (5)) with R₂ = R_RL,2, where as in binary_c and SeBa the instantaneous radius at the start of the double CE-phase is taken for the secondary radius. This can have a significant effect on the orbit of the post-double CE-system, leading to an increase of systems at low separations (approximately 1 R_⊙) in the binary_c and SeBa data compared to the StarTrack data.

Appendix A.1.3: Channel 3: stable case B

Evolutionary path: for channel 3, mass transfer starts when a hydrogen shell burning star fills its Roche lobe in a stable way before core helium-burning starts (Kippenhahn & Weigert 1967, case Br). This can occur when the envelope is radiative or when the convective zone in the upper layers of the envelope is shallow. In this project we assume that stable mass transfer proceeds conservatively and so the secondary significantly grows in mass. Because mass transfer is conservative, the orbit first shrinks and when the mass ratio has reversed the orbit widens. Mass transfer continues until the primary has lost (most of) its hydrogen envelope. At this stage the primary can become a helium WD or, if it is massive enough, ignite helium in its core. In the latter scenario the primary is a He-MS star. Like for channel 2, there are two sub-channels depending on whether the primary star fills the Roche lobe for a second time as a helium star. If the primary does not go through a helium-star phase or does not fill its Roche lobe as a helium star, the system evolves according to channel 3a. In channel 3b there is a second phase of mass transfer.

Example of channel 3a: Fig. A.14a shows an example of the evolution of a binary system of channel 3a with initial parameters M_1,zams = 4.8  M_⊙, M_1,zams = 3  M_⊙ and a_zams = 70  R_⊙. The masses of He-MS and secondary star are very similar in the BPS codes [0.82,0.83,0.82,0.82] M_⊙ and [6.9,7.0,7.0,7.0]M_⊙ respectively. The separations at the moment the helium star forms are [4.2,4.3,4.3,4.7]     ×    10²R_⊙ and are similar as well. In the subsequent evolution, the primary star effectively evolves as a single helium star before becoming a carbon-oxygen WD and loses [0.038,0.14,0.043,0.038] M_⊙ during that time and the orbit does not change significantly. changes by [2.1„ 2.5,4.2] R_⊙.

Population from channel 3a: regarding channel 3a, not all codes agree on the ranges of separation and masses (Figs. A.19 and A.20). However, there is an agreement between binary_c, the Brussels code and SeBa that majority of intermediate mass systems originate from systems with M_1,zams between 3 and 5 M_⊙ and a_zams between 10 and 100 R_⊙. The SWD population at WD formation is centred around systems with M_1,swd ≈ 0.6 M_⊙ for the binary_c, Brussels and SeBa codes, and with the majority of separations between about 20 − 1000 R_⊙. The SWD systems and their progenitors that are just described are not SWD progenitors according to StarTrack. According to this code, mass transfer is unstable and the system merges. The birthrates of binary_c, the Brussels code and SeBa differ within a factor of about 4 (Table A.1). In addition binary_c, SeBa and StarTrack show a good agreement on the different sub-populations for the full mass range. At WD formation these codes show a subpopulation between 15 to about 200 R_⊙ with WD masses of between 0.17 and 0.35 M_⊙. There is a second subpopulation at about 1 R_⊙ with most systems having a WD between 0.4 and 0.8 M_⊙. A third population shows mainly WD masses of more than 0.8 M_⊙ at separations of more than 300 R_⊙, where the population is extended to higher separations and WD masses in the results of SeBa and StarTrack. The third population is also visible in the progenitor population in Fig. A.19 with primary masses of more than 5 M_⊙ and separations of more than about 70 R_⊙. Again this population is more extended to high masses and separations according to SeBa and StarTrack. The low mass range of the progenitor population shows predominantly systems in orbits of 5−15 R_⊙. SeBa and StarTrack agree that there is an extra group at high orbital separations a_zams ≈ (1.3 − 4.6)    ×    10²R_⊙.

Example of channel 3b: an example of the evolution in channel 3b is shown in Fig. A.14b. Initially the system has M₁ = 7  M_⊙, M₂ = 5  M_⊙ and a = 65  R_⊙. After the first phase of mass transfer the primary masses M₁ =  [1.4,1.5,1.4,1.4] M_⊙, the secondary masses M₂ =  [11,11,11,11] M_⊙ and separations a =  [3.8,3.3,3.8,4.1]     ×    10²R_⊙. When the primary fills its Roche lobe again, it has lost [5.8, − ,6.8,7.3]     ×    10^-2M_⊙ in the wind. The mass transfer phase is stable and the secondary increases in mass to [11,11,11,11] M_⊙. The primary becomes a WD of [1.1,1.0,0.99,1.0] M_⊙ in an orbit of [4.5,6.5,5.9,6.2]     ×    10²R_⊙.

Population from channel 3b: the binary_c, Brussels and SeBa codes agree well on the initial systems leading to SWDs through channel 3b. This holds for both the initial mass, namely between about 5 and 10 M_⊙ and the initial separation between 0.1 − 3.0    ×    10²R_⊙. The population of progenitors of channel 3b according to the StarTrack code lies inside the previously mentioned ranges, however, the parameter space is smaller. In addition the four codes agree that at WD formation the majority of companions that are formed through channel 3b are massive, about 6 to 18 M_⊙ (for StarTrack 8−18 M_⊙.) The orbits of these systems are wide around 10³R_⊙, however, the ranges in separation and WD mass differs between the codes and will be discussed in the next paragraphs. Binary_c, SeBa and StarTrack also show a group of lower mass companions. For binary_c and SeBa these lie in the range 0.8−4.5 M_⊙ with separations of 0.5−30 R_⊙ and M_1,swd mainly between 0.6 and 1.0 M_⊙. The population of StarTrack agrees with these ranges, however, the parameter space for this population is smaller.

Effects: the population of SWDs from channels 3a and 3b are influenced by the MiMwd-relation. An important contribution to the MiMwd-relation comes from the assumed mass losses for helium stars and its mechanism, i.e. in a fast spherically symmetric wind or in planetary nebula (Appendix B). There is not much known about the mass loss from helium stars either observationally or theoretically. The differences in the MiMwd-relation affect for example the distribution of separations in Fig. A.16. For channel 3b the separation is ≲1400 R_⊙ for binary_c, SeBa, and StarTrack, but is extended to 6600 R_⊙ in the Brussels code. Binaries become wider in the Brussels code, as the WD masses in channel 3 are in general smaller compared to the other three codes.

There is also a difference in the MiMwd-relation between StarTrack on one hand, and binary_c and SeBa on the other hand regarding primaries that after losing their hydrogen envelopes become helium stars. For massive helium stars, binary_c and SeBa find that these stars will collapse to neutron stars, where as in StarTrack these stars form WDs. For channel 3a the difference occurs for the range of helium star masses of 1.6−2.25 M_⊙. As a result, systems containing massive helium stars are not considered to become SWD systems in binary_c and SeBa. These systems are present in the SWD data of StarTrack at M_1,swd ≳ 1.38 M_⊙ in Fig. A.3 for channels 3a and 3b. The progenitors lie at M_1,zams ≳ 8 M_⊙ with mostly a_zams ≈ 65 − 220 R_⊙ for channels 3a and 3b.

Another effect on the SWD population is the modelling of the mass transfer phases which is inherent to the BPS codes. The value of the mass transfer rate or the length of the mass transfer phase, however, do not have a large effect on the population or the evolution of individual systems from channel 3b in the set-up of the current study. This is because a priori conservative mass transfer is assumed, and therefore the accretion efficiency is not affected by the mass transfer rate. The mass transfer timescale only affects the binary evolution when other evolutionary timescales (such as the wind mass loss timescale or nuclear evolution timescale) are comparable. For example, while for M₁ ≪ M₂ the orbit increases strongly during RLOF, the orbit increases only moderately during wind mass loss assuming Jeans mode angular momentum loss. The range of separations in Fig. A.16 is therefore, besides the MiMwd-relation, also affected by the amount of wind mass and wind angular momentum leaving the system during RLOF. The binary_c, SeBa, and StarTrack codes assume wind mass takes with it the specific angular momentum of the donor star (Jeans mode), where as the Brussels code does not take wind mass loss into account during stable mass transfer.

Generally, no significant evolution of the donor star takes place during the mass transfer phase. Therefore with the current set-up, the post-mass transfer masses are determined by their initial mass and for binary_c, SeBa and StarTrack also the evolutionary moment the donor star fills its Roche lobe. However, an exception to this occurs for channel 3b during the second phase of mass transfer. Here the length of the mass transfer phase is important, as the evolutionary time scale of an evolved helium star is very short (of the order of few Myr) and the core grows significantly during this period. As a result the duration of the mass transfer phase becomes important for the resulting WD mass and separation in binary_c, SeBa and StarTrack (e.g. the example of channel 3b).

A crucially important assumption for the evolutionary outcome of channel 3 are the adopted stability criteria. Despite the importance of the stability criteria, the various implementations have not been compared until this study. We find a clear disagreement between the codes; stable mass transfer is possible in systems with mass ratios q_zams ≳ 0.6 according to StarTrack, in SeBa q_zams > 0.35, in binary_c q_zams > 0.25 and q_zams > 0.2 in the Brussels code. The effect of the relative large critical mass ratio for StarTrack results in a low birthrate in particular in the intermediate mass range (Table A.1 and Fig. A.20), which results in a lack of SWD systems lower than 300 R_⊙ (Fig. A.15) and fewer SWD systems with M_2,swd ≲ 1.0 M_⊙ (Fig. A.17). The effect of the relative low critical mass ratio for the Brussels code can be seen in Fig. A.18 as an extension in the Brussels code to lower separations a_swd ≲ 50. Systems with lower mass companions initially, go through mass ratio reversal and subsequent expansion of the orbit later in the mass transfer phase.

In the low mass range we find that the stability criteria vary most strongly for donor stars that are early on the first giant branch when they have shallow convective zones in the upper layers of the envelope. In general, stable mass transfer from this type of donor for the same conditions is more readily realised in StarTrack than in binary_c, and it is even more readily realised in SeBa. Systems with this kind of donor show in Fig. A.19 at M_1,zams < 3 M_⊙ a larger range in initial separations for SeBa (5−25 R_⊙) than for binary_c and StarTrack (5−18 R_⊙). There is also an extra population of SWD systems in the SeBa and StarTrack data with high initial separation a_zams ≈ (1.3−4.6)    ×    10²R_⊙ and high initial mass ratio q_zams ≈ M_2,zams/M_1,zams > 0.8. In these systems the primary fills its Roche lobe stably on the giant branch after the mass ratio has reversed due to wind mass losses. When donors with shallow convective zones are excluded, the birthrate in the full mass range in channel 3a decreases to 1.4  ×  10^-3   yr^-1 for SeBa and 7.3    ×    10^-4   yr^-1 for StarTrack, which is comparable to the birthrate predicted by binary_c (Table A.1).

The long-term behaviour of the orbit can be effected by tides. If energy is dissipated, angular momentum can be exchanged between the orbit and the spin of the stars. For this project the binary_c, Brussels and SeBa code assume that the spin angular momentum of the stars can be neglected compared to the orbital angular momentum³. As such in their simulations orbital angular momentum is conserved. In the StarTrack code, the orbital and spin angular momentum combined are conserved, under the assumption that the stars are and remain in a synchronized orbit. As a consequence after the first phase of mass transfer in channel 3, the orbits are slightly larger in StarTrack compared to those of the other codes (see the example system of channels 3a and 3b).

Whether or not a primary fills its Roche lobe for a second time is modelled different in the Brussels code than in the other three codes. In the Brussels code stars with an initial mass less than 5 M_⊙ are assumed to evolve through channel 3a, and stars with a higher mass evolve through channel 3b. The binary_c and SeBa simulations roughly agree with this, see Fig. A.19. However, the boundary between channels 3a and 3b is determined at run time in binary_c, SeBa, and StarTrack. It is dependent on the evolution of the radii and wind mass loss of helium stars, the stability criterion and the separation after the first phase of mass transfer. Therefore differences exist between these codes in the upper limits for ZAMS masses and separations in channel 3a in Fig. A.19 and at WD formation in Fig. A.17. Binary_c, SeBa, and StarTrack also find that systems that evolve through channel 3a or 3b overlap in WD mass at WD formation (Figs. A.16 and A.18). In the data from the Brussels code, the boundary at WD formation is discontinuous in primary mass causing a gap between 0.7 and 0.9 M_⊙ (Figs. A.16 and A.18). The gap in WD mass in the data from the Brussels code originates as a considerable amount of mass is lost during the planetary nebula phase of a star that does not initiate a second mass transfer phase. In the other three codes, the mass loss in winds from helium stars is less strong compared to the mass loss in the planetary nebula phase of helium stars in the Brussels code.

The evolution of helium stars (their radii, core masses, wind mass losses, and if they fill their Roche lobes also the stability and mass transfer rates) are important in channel 3. A difference arises between the Brussels code and the others, because of the way helium stars are simulated. In binary_c, SeBa, and StarTrack it is possible that after the first phase of mass transfer, the secondary fills its Roche lobe before the primary moves off the He-MS and becomes a white dwarf. Subsequently the primary becomes a WD before the secondary evolves significantly⁴. This reversal can occur because the evolutionary timescale of a low-mass helium star is very long (about 10⁸yr), while that of the secondary that gained much mass is reduced. As a result, when the first WD is formed, the mass of the secondary and the orbital separation has decreased substantially. These systems lie according to binary_c, SeBa and StarTrack at separations ≲20 R_⊙, primary WD masses of ≲1.0 M_⊙ and secondary masses of ≲4.5 M_⊙ in Figs. A.15 and A.17. The birthrates of the systems in binary_c, SeBa and StarTrack, are low ([1.1, − ,8.6,0.4]     ×    10^-4   yr^-1 in the full mass range). In the Brussels code, the evolution of the stars is not followed in time, and this evolutionary track is not considered. As a result in the range of 0.45−0.7 M_⊙ for the WD mass, the range in secondary masses is broader in the Brussels code.

Appendix A.1.4: Channel 4: unstable case B

Evolutionary path: in this path, a hydrogen shell burning star fills its Roche lobe (Kippenhahn & Weigert 1967, case Bc), but the mass transfer is unstable. After the CE-phase the primary becomes a helium WD or a He-MS. Again, we differentiate two evolutionary paths within a channel. In channel 4a, the primary becomes a WD directly or the primary becomes a helium star that will evolve into a WD without any further interaction with the secondary. If the primary star fills its Roche lobe for a second time, the system evolves through subchannel 4b. An example of the evolutionary path of channel 4b in shown in Fig. A.21.

Example: Fig. A.21 shows the evolution of a system of channel 4b, that starts its evolution with M_1,zams = 6  M_⊙, M_1,zams = 3  M_⊙ and a_zams = 320  R_⊙. The primary fills its Roche lobe as it ascends the first giant branch. After mass transfer ceases the primary has become a He-MS of mass [1.1,1.1,1.1,1.1] M_⊙ in an orbit with a separation of [7.0,7.1,7.0,7.0] R_⊙. As the helium star evolves and increases in radius, it initiates the second phase of mass transfer. Soon after mass transfer ceases, the primary becomes a WD with M_1,swd =  [0.81,0.91,0.77,0.79] M_⊙. The secondary is still on the MS with M_2,swd =  [3.2,3.2,3.3,3.3] M_⊙ and the orbital separation is a_swd =  [10,9.4,11,11] R_⊙. The differences in this example are caused by effects discussed before; the MiMwd-relation including the mass transfer rates from helium rich donors.

Population: the codes agree well on the location of the SWDs at WD formation from channels 4a and 4b in Figs. A.1–A.4, their progenitor systems in Figs. A.5 and A.6 and the birthrates of the channels (Table A.1). For channel 4a, which predominantly contains low mass binaries, there is an excellent agreement between binary_c, SeBa and StarTrack in the previously mentioned figures as well as in the birthrates (Table A.1). The low mass SWDs at WD formation have WDs of 0.25 − 0.48 M_⊙, companions of <1.8 M_⊙, in an orbit of 0.5 − 100 R_⊙, and progenitor systems with a_zams ≈ (0.3 − 4.0)    ×    10²R_⊙ for M_1,zams ≈ 1 − 2 M_⊙.

The population of systems that evolve through channel 4b are primarily intermediate mass binaries of mass M_1,zams ≈ 4.5 − 10 M_⊙ that become WDs of M_1,swd ≈ 0.7 − 1.3 M_⊙. The majority of systems have initial separations of a_zams ≈ (0.2 − 1.0)    ×    10³R_⊙. At WD formation the range of separations according to binary_c, SeBa and StarTrack is 4 − 1.0    ×    10²R_⊙, however, for the Brussels code it is extended to 0.9 − 1.4    ×    10²R_⊙.

Effects: there is a difference between the Brussels code on one hand and the other three codes on the other hand, regarding the survival of systems with low initial secondary masses M_2,zams < 3 M_⊙ in channel 4b. This is predominantly due to the difference in the single star prescriptions for the radii of stars. The radius of low-mass secondary-stars are in general larger in binary_c, SeBa and StarTrack than in the Brussels code. Therefore in the former three codes, the stars are more likely to fill their Roche lobe at the end of the CE-phase resulting in a merger. In the Brussels data, these systems survive at small separations (≲10 R_⊙ at 1 M_⊙, see Fig. A.2). Note that the Brussels code was written for intermediate mass stars (Sect. 3), and in principal the code does not allow for the detailed evolution of stars with initial masses below 3 M_⊙.

In addition, the stability of the second phase of mass transfer affects the SWD population of channel 4. If this phase is unstable, the system will evolve to a merger. In the Brussels code, it is assumed that the second phase of mass transfer is always stable, however, this is not the case in the three other codes. Differences in the stability criteria affect the orbital separation of SWD formation for all codes.

Appendix A.1.5: Channel 5: case A

Evolutionary path: in channel 5 mass transfer starts during the core hydrogen burning phase of the donor (Case A, Kippenhahn & Weigert 1967).

Population: the birthrates in the full mass range differ within a factor 2.5 between binary_c, SeBa and StarTrack (Table A.1). According to these codes, the progenitors of the primaries in channel 5 are stars of low mass (1−4 M_⊙) in small orbits (5−13 R_⊙), see Fig. A.5. There is a good agreement that the majority of SWDs from channel 5 at WD formation consists of a primary of mass 0.2−0.35 M_⊙, a secondary of mass 1.8−5.5 M_⊙ in an orbit with a separation of 30−240 R_⊙ (Figs. A.1 and A.3). Binary_c, SeBa and StarTrack further agree on a subchannel (a_swd ≈ 0.4 R_⊙ and M_1,swd ≈ 0.3 M_⊙ in Fig. A.1) in which the secondary is a hydrogen-poor helium-star at WD formation (see also channel 3). The birthrates of this subchannel are low ([4.5, − ,4.8,18]     ×    10^-6   yr^-1).

The birthrate of channel 5 in the Brussels code is higher by over a factor 20 compared to binary_c, SeBa and StarTrack (Table A.1). For the Brussels code the intermediate mass primaries have an initial mass M_1,zams ≈ 3 − 10 M_⊙ and WD mass M_1,swd ≈ 0.45 − 1.3 M_⊙, while the other codes show smaller ranges: for the main group of progenitors M_1,zams ≈ 3 − 4 M_⊙ and WD mass M_1,swd ≲ 0.35 M_⊙ (Figs. A.1 and A.5). The initial separation in the Brussels code a_zams ≈ 5 − 22 R_⊙, while in binary_c, SeBa and StarTrack a_zams ≈ 8 − 13 R_⊙. The separation at SWD formation a_swd in the Brussels code is between 20−350 R_⊙, while in the other codes the separation is mainly between 100−250 R_⊙. The range of secondary masses is M_2,swd ≈ 3 − 18 M_⊙ in the Brussels code, but only M_2,swd ≈ 4 − 6 M_⊙ in the other codes. Note that the region indicated by the dash-dotted contours in Fig. A.1, contains systems from channel 5 as well as from channel 3, however, this does not change our conclusion regarding the extended range and birthrates in the Brussels code compared to the other codes.

Effects: the differences between the Brussels code and the other codes is caused by the fact that the Brussels population code does not follow the mass transfer event and its mass transfer rate in detail. It considers only the initial and final moment of the mass transfer phase, therefore any intermediate steps in which the system can be closer are disregarded. For example, during conservative mass transfer to an initially less massive companion, the orbital separation first decreases and then increases again after mass ratio reversal. As the orbital separation decreases, the secondary can fill its Roche lobe leading to a contact system, especially as it grows in mass and radius due to the accretion. In the binary_c, SeBa and StarTrack code, it is assumed that the contact phase will lead to a merger or CE-phase for evolved secondaries. The Brussels code assumes that for shallow contact, the merger can be avoided. In other words, the codes have different assumptions for the stability of mass transfer.

Appendix A.2: Double white dwarfs

In the next sections, we differentiate four different evolutionary paths of DWDs. This is based on whether or not mass transfer occurs and if so, if the mass transfer initiated by the primary and secondary is stable or unstable. For clarity we do not distinguish the evolutionary path further e.g. by separating channels 3 and 5, nor channels 3a and 3b. Channels I, II and III represent the most commonly followed evolutionary paths with birthrates larger than 1.0    ×    10^-3   yr^-1. Channel IV is included because it stands out in Figs. A.22 and A.24, even though the birthrates in this channel are low (Table A.1). In each section we describe a specific evolutionary path (marked as Evolutionary path), we compare the simulated populations from each code (marked as Population) and investigate where differences between the populations come from (marked as Population and Effects).

Appendix A.2.1: Channel I: detached evolution

Evolutionary path: channel I involves non-interacting binaries. An example of a system was given for channel 1 in Sect. A.1.1: a 5 M_⊙ and 4 M_⊙ star in a circular orbit of 10⁴R_⊙. When the first WD is born, the orbit has increased to [1.8,1.8,1.8,1.8]     ×    10⁴R_⊙. When the second WD is born, the orbit has increased even more to [4.9,5.0,4.9,4.9]     ×    10⁴R_⊙ with primary and secondary masses of [1.0,0.94,1.0,1.0] M_⊙ and [0.87,0.86,0.87,0.87] M_⊙ respectively.

Population: there is a good agreement between the codes on the separations and masses of non-interacting DWDs, initially and at DWD formation. In the full mass range, initial separations are a_zams ≈ (0.5 − 10)    ×    10³R_⊙. The codes binary_c, SeBa, and StarTrack find non-interacting DWDs with WD masses between 0.5−1.4 M_⊙ in wide orbits of a_dwd ≈ (0.1 − 5.4)    ×    10⁴R_⊙. In the intermediate mass range, the initial separations are a_zams ≳ 1.5    ×    10³R_⊙ for binary_c, SeBa and StarTrack. Both WD masses are ≳0.75 M_⊙ and orbits are wide with separations a_dwd ≈ (0.6 − 5.4)    ×    10⁴R_⊙ for binary_c, SeBa and StarTrack. For the Brussels code, the separations are slightly higher at a_zams ≳ 2.8    ×    10³R_⊙ and a_dwd ≈ (1.3 − 7.2)    ×    10⁴R_⊙, and both WD masses extend to slightly lower values of ≳0.65 M_⊙. Small differences between the populations are due to different MiMf-relations and different prescriptions for single stars (e.g. stellar radii), as for SWDs from channel 1. The birthrates in channel I are very similar in the full mass range as well as in the intermediate mass range (Table A.1).

Appendix A.2.2: Channel II: CE + CE

Evolutionary path: the classical formation channel for close DWDs involves two CE-phases. First the primary star evolves into a WD via a phase of unstable mass transfer, i.e. via the evolutionary path described in Sects. A.1.2 and A.1.4 as channel 2 or 4 respectively. Subsequently the secondary initiates a CE-phase. It should be noted that for binary_c, SeBa and StarTrack this channel includes, systems that evolve through one CE-phase in which both stars lose their (hydrogen) envelope, the so-called double CE-phase described in Sect. 2 and in Eq. (6). Note that in the Brussels code, the double CE-phase is not considered.

Population: in the full mass range there is a good agreement between the progenitors according to the binary_c and SeBa code and a fair agreement with the StarTrack code. These three codes find that primaries of M_1,zams ≈ 1 − 8 M_⊙ contribute to this channel. For the majority the primaries have initial separations of a_zams ≈ (0.1 − 2.5)    ×    10³R_⊙. The DWD populations as predicted by binary_c and SeBa are similar, and comparable with the population of StarTrack. WD masses range from 0.35−1.4 M_⊙ for primaries and 0.19−0.9 M_⊙ for secondaries for binary_c and SeBa, and slightly larger ranges for StarTrack of 0.2−1.4 M_⊙ for primaries and 0.1−0.8 M_⊙ for secondaries. The orbital separation of DWDs from channel II is between a few tenths of solar radius to a few solar radii, however, the specific ranges of the three codes differ. The birthrates in channel II are similar between the three codes (Table A.1). In the intermediate mass range the codes agree that primaries and secondaries with initial mass between about 3 to 8 M_⊙ can contribute to channel II. In the Brussels code the mass range is slightly extended to higher masses of 10 M_⊙ for primaries due to the MiMwd-relation. There is an agreement on the initial separation of the majority of system, although the range of separations differs between the codes. For binary_c and SeBa a_zams ≈ (0.7 − 2.5)    ×    10³R_⊙, however, the range for StarTrack is extended to lower values of a_zams ≈ (0.4 − 20)    ×    10²R_⊙ as noted above. Comparing with the Brussels code, the range is extended to lower as well as higher values ((0.3 − 3.2)    ×    10³R_⊙). The higher maximum initial separation depends on the maximum radius in the single star prescriptions as discussed in channel 1. The difference in the lower minimum initial separation for the Brussels code has been noted for the SWDs in channel 2 as well. The Brussels code assumes that the primaries in these systems become WDs without a second interaction, where as in binary_c, SeBa and StarTrack these systems merge in the second interaction of the primary star. The separations of DWDs are centred around 0.5 R_⊙, however, the distribution of separations is different between the codes: 0.17−10 R_⊙ for binary_c, 0.06−1.18 R_⊙ for the Brussels code, 0.14−3.6 R_⊙ for SeBa and 0.05−11 R_⊙ for StarTrack. Primary WD masses are comparable between the codes, [0.8 − 1.4,0.5 − 1.3,0.7 − 1.4,0.7 − 1.4] M_⊙ where the ranges are the largest for the Brussels code. The maximum WD mass in the Brussels code is lower compared to the other codes due to the MiMwd-relation, see channel 1. The secondary WD masses at a given primary WD mass are lower in binary_c, SeBa and StarTrack (≲0.9 M_⊙) compared to the Brussels code (≲1.3 M_⊙).

Effects: several effects influence the distribution of separations in Fig. A.23. Even though the codes agree that the majority of DWDs from channel II have separation around 0.5 R_⊙, the spread around this value varies between the codes. In the full mass range the maximum separation is 8 R_⊙ in the SeBa data, 22 R_⊙ in the StarTrack data and 31 R_⊙ for binary_c. In the intermediate mass range it is 1 R_⊙ for the Brussels results, 4 R_⊙ in the SeBa data, 10 R_⊙ for binary_c and 11 R_⊙ in the StarTrack data. The maximum separation is affected by the MiMwd-relation and winds. As seen in channel 2, the maximum orbital separation in the Brussels code is lower as winds are not taken into account and more mass is removed during the CE. The distribution of orbital separation in the Brussels data is also affected in a different way than in the others codes as this code assumes that AGB donors become WDs directly without a second phase of interaction (see also channel 2). In binary_c, SeBa and StarTrack AGB donors can become helium stars, that fill their Roche lobes for a second time, resulting in lower average masses. This effect can be seen in Fig. A.25 where the secondary mass in binary_c, SeBa and StarTrack is ≲0.9 M_⊙ where as it is extended to ≲1.3 M_⊙ in the Brussels data. Mass loss in combination with the stability criteria, as also discussed for channel 2 causes high separations in the binary_c data. However, the relatively high maximum separations found by the StarTrack code is not affected much by the difference in the MiMwd-relation and winds, but are affected by differences in the double CE-formalism (see below).

All codes find that initially many DWD systems have high mass ratios, that in binary_c, SeBa and StarTrack lead to a double CE-phase. As discussed for channel 2, there is a difference in the formalism of the double CE-phase between StarTrack on one hand, and binary_c and SeBa on the other hand. As a result the separation after the double CE-phase is smaller according to the latter two codes, and a merger is more likely to happen. The birthrates of systems in the full (intermediate) mass range that evolve through a double CE-phase is 7.2    ×    10^-4   yr^-1  (7.9    ×    10^-5   yr^-1) according to StarTrack, while it is 4.6    ×    10^-5   yr^-1  (2.5    ×    10^-5   yr^-1) and 1.1    ×    10^-4   yr^-1  (3.2    ×    10^-5   yr^-1) for binary_c and SeBa respectively. An example of systems that merge according to binary_c and SeBa, but form a DWD according to StarTrack are the systems at a_zams ≲ 120 R_⊙ in Fig. A.26 which lie at a_dwd ≈ 0.07 − 1.2 R_⊙ for M_1,dwd ≲ 0.35 M_⊙ in Fig. A.22. An example of systems that survive according to all codes, however, at smaller separations for binary_c and SeBa compared to StarTrack, are systems with M_1,dwd ≳ 0.7 M_⊙ and a_dwd ≈ 4 − 10 R_⊙ for StarTrack and a_dwd ≲ 2 R_⊙ for binary_c and SeBa.

An effect that plays a role in channel II concerns the survival of a system during the mass-transfer event. As explained for channel 2, BPS codes compare the radius of the stripped star (i.e. exposed cores) to the corresponding Roche lobe to determine whether or not a merger takes place during the CE-event. For donor stars that become WDs directly after mass transfer ceases, i.e. without a hydrogen-poor helium burning phase, the Brussels and StarTrack code assume a zero-temperature WD where as binary_c and SeBa assume the exposed core is expanded due to previous nuclear shell burning. The effect of this is that the radius of the stripped star is a factor of about 5 smaller in the Brussels and StarTrack code than in binary_c and SeBa. Therefore a merger is less likely to take place. Therefore the minimum separation in the intermediate mass range is 0.06 and 0.05 R_⊙ in the Brussels code and the StarTrack code, respectively. While the minimum separation is about 0.15 R_⊙ in binary_c and SeBa.

Appendix A.2.3: Channel III: stable + CE

Evolutionary path: in channel III, the primary initiates stable mass transfer (alike channels 3 or 5 which are described in Sect. A.1.3 and A.1.5). When the secondary fills its Roche lobe mass transfer is unstable⁵.

Population: there is an agreement between the codes about the main parameter space occupied by the DWDs from channel III and their progenitors, however, the codes do not agree completely. The causes for differences in channel III have been discussed previously in the context of SWDs (see the discussion on channels 3 and 5), but they lead to more pronounced differences in the DWD population than in the SWD population.

In the intermediate mass range, the binary_c, Brussels and SeBa code agree on the orbital characteristics of the main progenitors, M_1,zams ≈ 4 − 9 M_⊙ and a_zams ≈ (0.2 − 2)    ×    10²R_⊙. There is also a rough agreement between these codes on the range of masses of both WDs. For primaries binary_c and SeBa find M_1,dwd ≳ 0.65 M_⊙ and the Brussels code M_1,dwd ≳ 0.45 M_⊙ due to differences in the MiMwd-relation. For secondaries these three codes find M_2,dwd ≳ 0.7 M_⊙. The maximum mass of the primary and secondary WDs varies between 1.2 and 1.4 M_⊙. The birthrates are high (a few times 10^-4   yr^-1) in this channel according to binary_c, the Brussels code and SeBa, however, the birthrate is a factor 1000 lower according to StarTrack. In the StarTrack simulation there are only two systems in channel III in the intermediate mass range, and therefore we refrain from showing contours for this channel for the StarTrack data in Figs. A.23, A.25 and A.27. Figure A.27 shows an increase of progenitor systems at separations a_zams ≲ 20 R_⊙ and primary masses M_1,zams ≈ 3 − 5.5 M_⊙ in the Brussels simulation compared to those from the other codes. The effect carries through into the DWD population as the increase of systems in the data from the Brussels code with WD primary masses between 0.45−0.7 M_⊙. The orbital separation of DWDs in channel III is very similar between binary_c and SeBa, a_dwd ≈ 0.1 − 1.1 R_⊙, however, for the Brussels code a_dwd ≈ 0.3 − 20 R_⊙. The existence of wide systems in the Brussels code is not surprising, as this code also finds the widest SWDs from channel 3 in comparison with binary_c and SeBa. As discussed previously in Sect. A.1.3, this is related to differences in the MiMwd-relation and angular momentum loss from winds. The gap at M_1,dwd ≈ 0.7 − 0.9 M_⊙ in Fig. A.25 in the data from the Brussels code, is caused by the boundary between channels 3a and 3b, as in Fig. A.18 (Sect. A.1.3).

Regarding the populations of progenitors for low mass primaries, binary_c, SeBa and StarTrack agree reasonably well. They both show that most DWDs in channel III have initial separations of 5−20 R_⊙. However, the range of initial separations is extended to 25 R_⊙ in the population simulated by SeBa. SeBa and StarTrack also show an extra population compared to binary_c (a_zams ≈ 140 − 270 R_⊙ and M_1,zams ≲ 1.2 M_⊙). These two differences are due to differences in the stability of mass transfer for donors with shallow convective envelopes, as discussed for channel 3. Comparing the population of DWDs itself for low mass primaries, binary_c, SeBa and StarTrack agree well. The codes show a population of DWDs with primary mass M_1,dwd ≈ 0.2 − 0.44 M_⊙ at a separation of a_dwd ≈ 0.1 − 1.5 R_⊙, with secondary masses M_2,dwd around 0.6 M_⊙. The extra population in the SeBa and StarTrack data lies at a_dwd ≈ 10 − 50 R_⊙ and M_1,dwd ≈ 0.4 − 0.47 M_⊙. The three codes show systems at M_2,dwd about 0.3 M_⊙, where in the binary_c data this group is extended to higher primary WD masses of M_1,dwd ≈ 0.2 − 0.7 M_⊙ in Fig. A.24. These systems in binary_c mainly evolve through a specific evolutionary path in which there is a phase of stable mass transfer between a He-MS and a WD, a so-called AM CVn-system. The birthrate of these systems is 5.0    ×    10^-4   yr^-1 according to binary_c and negligible according to the other codes.

Effects: the extremely low birthrate of StarTrack in the intermediate range is caused by a combination of effects discussed previously. Firstly, stable mass transfer is more readily realised in the other codes compared to StarTrack (see channel 3). Only systems with q_zams ≳ 0.6 undergo stable mass transfer and become SWD systems according to StarTrack. For about 60% of these systems, the secondary becomes massive enough to collapse to a neutron star after nuclear burning ceases (in accordance with the other codes). Secondly, the remaining systems merge when the secondary star fills its Roche lobe. For AGB donors this is more likely to happen in the StarTrack data, because of the difference in the radii of stripped stars compared to binary_c and SeBa (see channel 2b).

The different methods of calculating mass transfer between the Brussels code and the other codes, cause an increase in systems in the data from the Brussels code, similar to channel 5. In particular for DWDs, it is important how the secondary responds to mass gain. The systems that survive in the Brussels code have q_zams > 0.85, such that the orbit widens severely due to the mass transfer. However, according to binary_c and SeBa, when the secondary accretes a significant amount of mass and is rejuvenated, its evolutionary timescale is reduced. As the secondary evolves, the system comes in contact and merges. The Brussels code assumes that the merger can be avoided for phases of shallow contact.

The evolution of and mass transfer rates from evolved helium stars donors (see channel 3) are important for channel III. It affects the DWD systems with high masses of the primary progenitor and primary WD, see Fig. A.25 and A.27. The range of primary WD masses is extended to 1.2 M_⊙ according to SeBa, and 1.3 M_⊙ according to the Brussels code and 1.4 M_⊙ according to binary_c. Contrary to stable mass transfer from hydrogen rich donors, the core of evolved helium stars can grow significantly during stable mass transfer phases as the timescale for mass transfer can become comparable to the timescale of wind mass loss or nuclear evolution. If the mass transfer phase is relatively short, the core of the donor star does not have time to grow significantly and little mass is lost in the wind. With the assumption of conservative mass transfer, most of the envelope is then transferred to the secondary star which then is more likely to become a neutron star instead of a WD.

Differences in the radii of stripped stars causes a relative lack of close systems for the Brussels code compared to the other codes. For channel II this was discussed in the context of donor stars that become a WD directly. However, in channel III in the intermediate mass range many donor stars continue burning helium after the mass transfer event ceases. The radius of the stripped donor star depends on its mass, and for binary_c, SeBa and StarTrack also on the evolutionary state of the donor stars (see also channel 2). When the donor star is stripped of its envelope before the AGB-phase, the core radius is a factor of about 4−5 larger in the Brussels code compared to binary_c, SeBa and StarTrack. Therefore a merger is more likely to take place in the Brussels code.

Appendix A.2.4: Channel IV: CE + stable

Evolutionary path: in the final evolutionary channel for DWDs, when the primary fills its Roche lobe, mass transfer proceeds in an unstable manner (according to channels 2 or 4 which are described in Sects. A.1.2 and A.1.4). However, when the secondary fills its Roche lobe mass transfer mass to the primary is stable. As a result the primary accretes mass.

Population: the systems of channel IV lie in a small and specific region of DWD parameter space (Figs. A.22 and A.24). The birthrates are low, 1.6    ×    10^-4   yr^-1 and 5.5   ×    10^-5   yr^-1 for binary_c and SeBa respectively in the full mass range. We do not compare the population of this channel with the Brussels code as the progenitors according to binary_c and SeBa are low mass binaries and the birthrate in the Brussels code is zero per definition (Sect. 4.1). We cannot compare the characteristics of the population of binary_c and SeBa with that of the StarTrack code as in the simulations of the latter code there are no systems evolving through channel IV indicating a birthrate of <4 × 10^-7   yr^-1. The birthrate is low according to StarTrack as unstable mass transfer is more readily realised in this code compared to binary_c and SeBa (see also channel 3). The binary_c and SeBa code agree well on the binary parameters of the population of DWDs at DWD formation from this channel: separations of 10−30 R_⊙ and primary WD masses of 1.1−1.4 M_⊙, and secondary WD masses of 0.15−0.20 M_⊙. The progenitors systems in this channel are similar, M_1,zams ≈ 1 − 3 M_⊙ and a_zams ≈ 50 − 400 R_⊙. Differences in the population of DWD systems from this channel, their progenitors and the birthrates occur due to the uncertainty in the stability of mass transfer and the mass transfer rate (see also the discussion for channel 3). Note that in the current study we have assumed conservative mass transfer to all accretors, including WDs. This is not a physical picture, so a warning of caution needs to be given to trust the parameters of this population, nonetheless the similarities between the binary_c and SeBa codes are striking.

Appendix B: Backbones of the BPS codes

The structure of BPS codes can vary strongly, which complicates the process of comparing BPS codes. Some aspects of the code are relatively simple to adapt in order to let assumptions of different groups converge, where as other aspects are inherent to the code and are not straightforward to change. For example, where some codes use results from detailed single star evolution codes, written down in analytical formulae (e.g. Eggleton et al. 1989; Hurley et al. 2000) to compute stellar parameters, others use the results of detailed binary evolution codes – a grid over which one can interpolate – and those results are integrated into the population code (e.g. De Donder & Vanbeveren 2004). The inherent differences will create differences between the results of the different groups. The main differences are summarised in Table B.1 and a more complete overview is given below. For most of the points the influence on a population it not immediately clear, therefore their effects are discussed in Sect. 6.

Appendix B.1: Single star prescriptions

The single star prescriptions, either given by analytical formulae or included in a grid of binary systems over which can be interpolated, determine which mass the WD will have when the star loses its envelope. Furthermore they determine the radii during the evolution of the star and therefore the moment at which the star fills its Roche lobe.

• binary_c, SeBa, StarTrack: the codes use analytical fitting formulae (Hurley 2000) from detailed single star evolution tracks, with an overshooting constant δ_ov = 0.12 (based on Pols et al. 1998). In binary_c different AGB models can be used, based on detailed models of Karakas et al. (2002) for thermally pulsating AGB stars (TP-AGB). However, these are not used for this work. Prior to the work of Toonen et al. (2012), the single star prescriptions in SeBa were based on Eggleton et al. (1989).

  Table B.1

  Numerical treatments in the different codes which are inherent to them.

• Brussels code: intermediate mass single star prescriptions are taken from Schaller et al. (1992). These tracks include convective overshooting by means of the following prescription: the overshooting distance d_over is directly proportional to the pressure scale height H_p according to d_over = 0.2H_p. This corresponds to a slightly lower degree of overshooting than in the codes that use the overshooting constant δ_ov = 0.12 in the stability criteria, the latter corresponding to a d_over/H_p between 0.22 and 0.4 depending on mass (see Hurley 2000). Stellar parameters which do not depend on whether the star is part of an interacting binary system are taken directly from this reference. Other stellar parameters, such as the remnant mass after RLOF, are taken from the detailed binary evolution code.

Appendix B.2: Stability of mass transfer

At the moment that one of the stars fills its Roche lobe mass transfer can proceed in a stable manner or the system can evolve into a CE-phase (Sect. 2). In the simulation of the evolution of a binary system the entire stellar structure is not explicitly followed in detail, and consequently, “stability checks” must be built-in to BPS codes to determine if RLOF will lead to a CE-phase.

• binary_c: for every type of donor star and type of accretor star acritical mass ratio (q_crit) is given. The mass ratio of the system during mass transfer is compared with the critical mass ratio for stable mass transfer and determines if mass transfer will proceed in a stable manner or not. An overview can be found in Claeys et al. (2013). Note that in that paper two possibilities are described for the stability of Roche lobe overflowing helium stars to non-degenerate accretors. For this project the criterion as described in Hurley et al. (2002) is used.

• Brussels code: the boundary between stable and unstable RLOF is determined by whether the outer layers of the donor star are radiative or deeply convective respectively. For each stellar mass, the minimum stellar radius R_conv is given for which the envelope is convective. If the orbital period of the system under investigation is smaller than the theoretical orbital period at the time when R_RL = R_conv, mass transfer will proceed in a stable way. If the mass ratio between the two stars is extreme (q = M_a/M_d < 0.2 = q_crit) at the onset of mass transfer, this will result in an instability (Darwin 1879). The donor star will be unable to extract sufficient angular momentum from the orbit to remain in synchronized rotation, resulting in the mass transfer episode quickly becoming dynamically unstable. Tidal interaction will cause the secondary to spiral into the donor’s outer layers, a process that is treated identically to the CE-evolution (hence with β = 0).

• SeBa: the stability and rate of mass transfer are dependent on the reaction to mass change of the stellar radii and the corresponding Roche lobes. The change in the Roche radius R_RL due to loss and transfer of mass M is given by (B.1)\newpagethe adiabatic (i.e. immediate) response of the donor star’s radius R is given by (B.2)For every Roche lobe filling system, ζ_RL and ζ_ad are compared at every timestep. If ζ_RL < ζ_ad we assume mass transfer proceeds in a stable manner (e.g. Webbink 1985; Pols & Marinus 1994). When ζ_RL > ζ_ad, mass transfer is dynamically unstable leading to a CE-phase. The value of ζ_RL is calculated numerically by transferring a test mass of 10^-5M_⊙. The advantage of this is that, because ζ_RL = ζ_RL(M_d,M_a,a) and so ζ_RL is dependent on the mass accretion efficiency of the secondary, the (de)stabilising effect (see Soberman et al. 1997) of non-conservative stable mass transfer is taken into account automatically. Appropriate recipes of ζ_ad are implemented in the code for every type of donor star. An overview can be found in Toonen et al. (2012, Appendix A3 therein). Furthermore, the orbital angular momentum is compared with the stellar spin angular momenta, to check whether a Darwin instability is encountered (Darwin 1879).

• StarTrack: when a non-degenerate star fills its Roche lobe, ζ_ad and ζ_RL are calculated, similar to the case of SeBa. The value of ζ_ad is determined by removing mass from the star over a 1-year timestep (Belczynski et al. 2008a). The value of ζ_RL is determined by transferring a small amount (1%) of the star’s mass toward the companion. In cases where the mass loss is so rapid such that the star loses thermal equilibrium, a “diagnostic diagram” is used to predict the stability of mass transfer (see description in Belczynski et al. 2008a, sect. 5.2). The diagnostic diagram is a numerical tool that was first calibrated using detailed stellar evolution calculations of massive stars, and is currently being updated to include a range of stellar models for low- and intermediate-mass stars. In addition, there is also a check for a possible delayed dynamical instability. This occurs for stars with M_d/M_a > q_ddi, with q_ddi based on Hjellming & Webbink (1987), or when a Darwin instability is encountered, or when the trapping radius of the accretion stream (King & Begelman 1999) exceeds the Roche radius of the accreting star (see Ivanova et al. 2003; Belczynski et al. 2008a, sect. 5.4). This latter point, however, is not considered for this work.

Appendix B.3: Stable mass transfer

To take into account various driving mechanism of stable RLOF, such as thermal readjustment or nuclear evolution of the donor, approximate prescriptions are used to determine the mass transfer rate. Note that mass transfer rate refers to the mass lost by the donor, which will always be equal to or greater than the mass accretion rate, which refers to the mass gained by the companion.

• binary_c: the mass transfer rate is calculated as a function of the ratio of the stellar radius and the Roche radius (based on Whyte & Eggleton 1980). A function is generated which follows the radius more closely during mass transfer on a thermal timescale and more loosely when the star is in thermally equilibrium. A smooth transition is built-in between the two. The formulation can be found in Claeys et al. (2013). That paper also shows that the resulting mass transfer phases are comparable to that of the detailed binary stellar evolution code STARS (based on Eggleton 1971) in the duration of the mass transfer phases and the mass transfer rates for a set of models. This method indirectly considers mass transfer on the nuclear and thermal timescale, but also on the timescale of gravitational wave radiation or magnetic braking are considered.

• Brussels code: the mass transfer rates are not explicitly calculated in the population code. It considers merely the initial and final masses. These are interpolated from the results of the detailed binary evolution code. The latter calculates the mass transfer rate during stable RLOF iteratively, by investigating how much mass needs to be lost during the current timestep for the donor star to remain confined by its Roche lobe (within the order of a few percent).

• SeBa: ζ_RL is compared with appropriate values of ζ_eq to determine if mass transfer is driven by the thermal readjustment or the nuclear evolution of the donor star. ζ_eq represents the response of the donor star’s radius R as is adjusts to the new thermal equilibrium: (B.3)Appropriate recipes of ζ_eq are implemented for every type of donor star. If ζ_RL < min(ζ_ad,ζ_eq), mass transfer is driven by the nuclear evolution of the donor star and we assume mass transfer proceeds on the nuclear timescale of the donor star (e.g. Webbink 1985; Pols & Marinus 1994). If ζ_eq < ζ_RL < ζ_ad, RLOF is dynamically stable and driven by thermal readjustment of the donor, so that mass transfer proceeds on the thermal timescale of the donor star. In addition, stable mass transfer can be driven by angular momentum loss from magnetic braking or gravitational wave emission. When the timescale of angular momentum loss is shorter than the mass loss timescale determined above, we assume mass transfer is driven by angular momentum loss. For more detail see Appendix A.3 of Toonen et al. (2012).

• StarTrack: for non-degenerate donors ζ_RL and ζ_ad are calculated, along with the thermal timescale τ_KH (based on Kalogera & Webbink 1996). Additionally, the equilibrium mass transfer timescale τ_eq is calculated as a combination of RLOF both driven by angular momentum loss and the nuclear evolution of the star and/or the changes due to magnetic braking and gravitational wave radiation (see Belczynski et al. 2008a). If τ_eq > τ_KH the mass losing star is in thermal equilibrium and mass transfer proceeds on Ṁ_eq = M/τ_eq. If τ_eq ≤ τ_KH mass transfer proceeds on a thermal timescale, given by Ṁ_KH = M/τ_KH. If Ṁ_eq becomes positive the star falls out of equilibrium and the stability of mass transfer is determined by the diagnostic diagram (see Belczynski et al. 2008a). In the case of WD donors, the mass transfer rate is always driven by gravitational radiation.

Appendix B.4: Wind mass loss

The driving mechanisms of the wind and the explicit rate at which this material is lost are not yet completely understood. This results in different prescriptions to describe the rate of wind mass loss and the amount that can be lost (e.g. Wachter et al. 2002). We only discuss the wind-prescriptions that are relevant for low and intermediate mass stars.

• binary_c, SeBa, StarTrack: for stars up to the early AGB theprescription of Reimers (1975) is adopted (withη = 0.5). To describe the wind mass loss of stars on the TP-AGB a prescription based on Vassiliadis & Wood (1993) is implemented. Both prescriptions are defined in Hurley et al. (2000). In binary_c and StarTrack different prescriptions for the wind mass loss are available that used by different users of the respective codes.

• Brussels code: for intermediate mass interacting binaries, the initial-final mass relation of WDs is determined by assuming the wind prescription of van den Hoek & Groenewegen (1997). However, it should be noted that in the BPS code a star in an interacting binary does not have wind mass loss. For the most massive stars, wind mass loss is as is described in Vanbeveren et al. (1998).

Appendix B.5: Angular momentum loss from winds

Section 2.1 and Appendix C.2 describe the importance of angular momentum loss (AML) and the effect on the orbit. Not only mass lost during RLOF, but also wind carries angular momentum, which is lost when it leaves the system. The same prescriptions as described in Appendix. C.2 can be applied to AML when material is lost through a wind and different prescriptions are used in the BPS codes.

• binary_c: different prescriptions of angular momentum loss through a stellar wind are available in binary_c. In this study, as in Claeys et al. (2013), wind angular momentum loss is as described in Hurley et al. (2002). When no material is accreted by the companion star, the wind takes specific angular momentum of the donor.

• Brussels code: mass lost by a stellar wind in non-interacting systems is lost through the Jeans mode. Interacting systems do not have wind mass loss prior to interaction.

• SeBa, StarTrack: the material lost by a wind that is not accreted by the companion is lost from the system with specific angular momentum from the donor.

Appendix B.6: Evolution of helium stars

A helium star is formed after a hydrogen-rich star with a helium core loses its hydrogen-rich envelope. When the core is not degenerate at that time, the evolution of the star continuous as a helium-burning star. Uncertainties in the evolution of helium stars encompasses the growth of this star, the wind mass loss and mass transfer phase, such as the stability and rate.

• binary_c, SeBa, StarTrack: the evolutionary tracks and windprescription are based on Hurleyet al. (2000). The stability of masstransfer and the rate are described in previous sections.

• Brussels code: helium star evolution is not explicitly included in the code. It is assumed that the donor star always loses its entire H-rich envelope in one episode and becomes a WD afterwards, except in the case where a donor fills its Roche lobe for a second time as a helium star. In this case mass transfer is followed as described in Sect. B.3, however, time-dependent evolutionary aspects of the helium star are not followed. This simplification is made because the intermediate step is not believed to have a large influence on the eventual masses and separation. However, this implicitly means that the most massive star will always become a WD first, which is not necessarily the case when helium star evolution is explicitly followed. For stars that lose mass during the planetary nebula phase, no resulting angular momentum loss is taken into account.

Appendix B.7: Generating the initial stellar population

The initial population can be chosen by a Monte Carlo method, or the choice can be made grid-based. Nevertheless, if the method is well performed both methods should give the same results for a high enough resolution.

• binary_c: N_M1,zams × N_M2,zams × N_azams binaries are simulated, with M_1,zams, M_2,zams, a_zams chosen in logarithmic space. A probability is calculated for every system determined by the defined initial distributions.

• Brussels code: the code works with a three-dimensional grid of initial parameters: primary mass M_1,zams, mass ratio q_zams and orbital period P_zams. According to the initial mass function, initial mass-ratio distribution and initial orbital period (or separation) distribution, each grid point is assigned a certain weight. Every system corresponding to such a grid point is then taken through its evolution.

• SeBa, StarTrack: initial parameters M_1,zams, M_2,zams, a_zams and the initial eccentricity e_zams are chosen randomly on a Monte-Carlo based-approach where the probability functions are given by the initial distributions. With this method, the resolution is highest in those regions of parameter space where most systems lie.

Some aspects of the codes that are not straightforward to change have been discussed in the previous section. However, other aspects of the codes are relatively simple to adapt. These aspects are often contained in relatively isolated and parametrised functions. For this project we equalised these aspects in the codes as far as possible. However, we do not believe that all the assumptions made for this project are realistic. Previous publications of results from these BPS codes are based on different assumptions. Although we do not compare the effect of the different assumptions on stellar populations in this work, it is good to realise which assumptions are generally used. Therefore the usual assumptions made by the authors in their corresponding BPS code are summarised in Table C.1 and are discussed in more detail below. Typical assumptions may vary between different users of the BPS codes.

Appendix C: Typical variable assumptions in BPS codes

Appendix C.1: Accretion efficiency

In this project mass transfer is assumed to be conservative to all types of stars. However, in general, the accretion efficiency depends on the type of accreting star and the mass transfer rate.

• binary_c, SeBa: in the case of non-degenerate accretors with radiative envelopes, the accretion efficiency mainly depends on the mass transfer rate and the thermal timescale of the accreting star. In the case of non-degenerate objects with convective envelopes, mass is transferred conservatively. In the case of a degenerate accretors, the accretion efficiency depends on the mass of the degenerate object and the mass transfer rate.

• Brussels code: the accretion efficiency onto a non-degenerate object is taken to be constant. If the mass ratio is below 0.2, mass transfer is unstable and the accretion efficiency is assumed to be zero (Sect. B.2). To ensure continuity, between mass ratios 0.2 and 0.4 a linear interpolation is used for the accretion efficiency, between 0 and β (usually 1). Note that for popcorn this transition was not implemented and the accretion efficiency is one between 0.2 and 0.4. In case of a degenerate accreting object, the regions in the (companion mass, orbital period)-parameter space from Hachisu et al. (2008) are used to determine in which cases the WD can stably accrete up to 1.4 M_⊙. In all other cases, mass transfer towards WDs is assumed to become unstable, and is treated as a CE-phase.

• StarTrack: the accretion efficiency onto a non-degenerate object is taken to be constant. In the case of a degenerate accreting object, the accretion efficiency depends on the mass of the accreting object and the mass transfer rate (see Belczynski et al. 2008a, sect. 5 therein).

Appendix C.2: Angular momentum loss during RLOF

In BPS codes a wide range of prescriptions are used to describe angular momentum loss when material is lost in a phase of stable RLOF. They can be divided in four modes of AML or combinations of these modes (see e.g. Soberman et al. 1997, for an overview of the effect of the different prescriptions on the stability of the system).

• Orbital angular momentum loss mode. In this mode the mass is assumed to leave the binary system, with (a multiple of) the specific orbital angular momentum of the binary, i.e η = const.

• Jeans mode. Mass is assumed to leave the system from the vicinity of the donor star in a fast spherically symmetric wind. In this mode, the wind matter does not interact with the system. It takes with it the specific orbital angular momentum of the donor in its relative orbit around the centre of mass. Making the assumption that the donor star can be approximated by a point mass, the specific angular momentum loss is as in Eq. (2) with: (C.1)

• Isotropic re-emission. In this case mass is assumed to leave the system from the position of the gainer in a spherically symmetric way (or at least symmetric with respect to the equatorial plane of the star). Possible scenarios are an enhanced stellar wind or bipolar jets. Further assumptions are as in the previous case, resulting in: (C.2)

• Circumbinary ring. Finally, it is possible to assume that the matter will leave the binary through the formation of a non-corotating circumbinary ring, after passing through the second Lagrangian point L₂. The amount of angular momentum lost then depends on the radius of this ring a_ring compared to the orbital separation a: (C.3)While an absolute minimum for a_ring is the distance from the center of mass to L₂ (which can be shown to vary only very slightly during a mass transfer episode), it was shown by Soberman et al. (1997) that a more realistic value is 2.3 times the orbital separation.

Figure C.1 shows, for the four different AML modes, the angular momentum loss as a function of mass ratio q = M_a/M_d. It is clear that the assumption of AML from a circumbinary ring always leads to the largest AML. The Jeans mode causes the least AML for systems with low mass ratios q < 1, because the donor is then close to the center of mass of the system. As the mass ratio increases during mass transfer, the AML increases as well since the donor recedes from the center of mass. Conversely, the isotropic re-emission mode causes a large AML for low mass ratio systems, as the gainer is far from the center of mass. As the mass ratio rises, the gainer closes in on the center of mass and AML decreases. The orbital AML assumption results in an intermediate case between the two.

The choice of AML mode is critical for both the stability and the orbital evolution of mass transfer. To illustrate, in the case of the circumbinary ring mode (extracting the most angular momentum), a given amount of mass loss will lead to much more AML than in the case of Jeans mode AML. The former mode will thus result in a far greater number of systems that merge than the latter.

The typical assumptions in the BPS codes are:

• binary_c: in this work and the standard model in Claeyset al. (2013), the material notaccreted during the stable RLOF phase is lost as isotropicre-emission.

• Brussels code: the material is lost through the second Lagrangian point such that angular momentum is lost from a circumbinary ring with a_ring = 2.3.

• SeBa: when the accretor is a non-degenerate star, the material lost carries 2.5 times the specific orbital angular momentum of the binary (Portegies Zwart 1995; Nelemans et al. 2001a). In the case of a degenerate accretor, the material lost carries specific orbital angular momentum of the accreting star.

• StarTrack: when the accretor is a non-degenerate star, the material lost carries one time the specific orbital angular momentum. In the case of a degenerate accretor, the material lost carries specific orbital angular momentum of the accreting star.

Appendix C.3: Common envelope evolution

There remain several uncertainties in the evolution of a CE-phase. For this reason, various BPS codes employ different CE-prescriptions (Sect. 2.2) and CE-efficiencies (such as α_ce) and both aspects are often varied within a BPS study for comparison. Here, we briefly describe the CE-parametrisations that are implemented most often by the authors in the four different codes.

• binary_c: to describe CE-evolution the prescription based onHurley et al. (2002) is used. In thestandard model of Claeyset al. (2013),α_ce is one, while λ_ce depends on the type of star, its mass and luminosity (see Izzard 2004; Claeys et al. 2013). However, in the BPS code also the γ-prescription can be used.

• Brussels code, StarTrack: for standard calculations, the prescription based on Webbink (1984) is used, where α_ce and λ_ce are both one. In both codes different values for α_ce and λ_ce can be implemented, as well as the γ-prescription (for further information about the version of the γ-prescription implemented in StarTrack see Belczynski et al. 2008a; Ruiter et al. 2011).

• SeBa: the standard model for simulating CE-evolution in SeBa is the γ-prescription, unless the binary contains a compact object or the CE is triggered by a Darwin instability (Darwin 1879) for which the α-formalism based on Webbink (1984) is used. The γ-formalism is introduced by Nelemans et al. (2000) in order to better reproduce the mass ratio distribution of observed DWDs. The mass loss reduces the angular momentum of the system according to: (C.4)where J_i and J_f are the angular momenta of the pre- and post-mass transfer binary respectively. The motivation for this formalism is the large amount of angular momentum available in binaries with similar mass objects that possibly can be used to expel the envelope. In SeBa γ is taken to be equal to 1.5, and α_ce × λ_ce is equal to two.

Appendix C.4: Wind accretion

Material that is lost in the form of a stellar wind can be partly accreted by the companion star. The amount depends on properties of the wind (e.g. the velocity), the accreting star and the binary system (e.g. the separation). However, the exact amount accreted is ill-constrained.

• binary_c, SeBa: the accretion efficiency of wind ma-terial is determined by the Bondi-Hoyle prescription(Bondi & Hoyle 1944). In binary_c theaccretion efficiency based on the wind Roche-lobeoverflow model can be used (Mohamed & Podsiad-lowski 2007, 2012; Abateet al. 2013), however, is not used forthis work.

• Brussels code: no material lost in the form of a stellar wind is accreted by the companion star.

• StarTrack: material lost through a wind is in general not accreted by the companion star, except when the companion star is a neutron star or a black hole.

Appendix C.5: Tides

The general picture of tidal effects is clear, however, uncertainties remain due to missing knowledge about for example some dissipative processes.

• binary_c: tidal evolution is implemented as described by Hur-ley et al. (2002), which is based onHut (1981); Zahn (1977).

• Brussels code: tidal evolution is a described by Zahn (1977).

• SeBa: tidal evolution is implemented as described by Portegies Zwart & Verbunt (1996).

• StarTrack: tidal evolution is implemented as described by Claret (2007), as well as Hurley et al. (2002), which is based on Hut (1981); Zahn (1977).

Appendix C.6: Magnetic braking

Magnetic braking is important for low mass stars with convective envelopes. Nevertheless, this process is not fully understood and different prescriptions co-exist.

• binary_c, SeBa: both codes use the prescription of Rappaportet al. (1983).

• Brussels code: the code is not used for the evolution of stellar objects with a mass lower than 3 M_⊙, therefore magnetic braking is not considered.

• StarTrack: the prescription of Ivanova & Taam (2003) is used in standard calculations.

Appendix C.7: Initial population

The choice for an initial distribution and the respective boundaries can severely affect the importance of a certain evolutionary channel through the normalisation of the simulation. For example, changing the upper boundary of the distribution of orbital separations from 10⁶R_⊙ to 10⁴R_⊙, increases the birthrates of interacting binaries by about 70%. On the other hand, changing the upper boundary for the primary mass distribution within 80 M_⊙ to 150 M_⊙ (see Table C.2), does not affect the normalisation.

The assumptions made by the authors with their respective codes are summarised in Table C.2. Different aspects which need extra clarification are discussed below. Note that other users of the BPS codes under study here, other than the authors, may use different distribution functions and/or ranges.

• binary_c: the initial eccentricity is zero, based on the work of Hurley et al. (2002). The minimum initial separation is varied between 5 R_⊙ or the minimum separation at which a binary system with a certain mass is initially detached. The minimum and maximum separations and secondary masses are based on the work of Kouwenhoven et al. (2007).

• Brussels code: the initial eccentricity is zero. No minimum and maximum separation is assumed for binaries, but a minimum and maximum initial orbital period, more specifically one day and 3650 days. In order to compare with the other codes, a conversion of orbital period to separations is given in Table C.2.

• SeBa, StarTrack: a distribution for the initial eccentricities is assumed (Table C.2). The initial semi-major axis is chosen between 10⁶R_⊙ and the minimum initial separation is the minimum separation at which a binary system with a certain mass is initially detached.

All Tables

All Figures