A&A 601, A29 (2017)
DOI: 10.1051/0004-6361/201629685
c
ESO 2017
Astronomy
&
Astrophysics
Delay-time distribution of core-collapse supernovae with late
events resulting from binary interaction
E. Zapartas
1
, S. E. de Mink
1
, R. G. Izzard
2
, S.-C. Yoon
3
, C. Badenes
4
, Y. Götberg
1
, A. de Koter
1, 5
, C. J. Neijssel
1
,
M. Renzo
1
, A. Schootemeijer
6
, and T. S. Shrotriya
6
1
Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands
e-mail: [E.Zapartas;S.E.deMink]@uva.nl
2
Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
3
Astronomy Program, Department of Physics and Astronomy, Seoul National University, 151–747 Seoul, Korea
4
Department of Physics and Astronomy & Pittsburgh Particle Physics, Astrophysics, and Cosmology Center (PITT-PACC),
University of Pittsburgh, Pittsburgh, PA 15260, USA
5
Institute of Astronomy, KU Leuven, Celestijnenlaan 200 D, 3001 Leuven, Belgium
6
Argelander-Institut für Astronomie, Universität Bonn, Auf dem Hügel 71, 53121 Bonn, Germany
Received 11 September 2016 / Accepted 1 January 2017
ABSTRACT
Most massive stars, the progenitors of core-collapse supernovae, are in close binary systems and may interact with their companion
through mass transfer or merging. We undertake a population synthesis study to compute the delay-time distribution of core-collapse
supernovae, that is, the supernova rate versus time following a starburst, taking into account binary interactions. We test the systematic
robustness of our results by running various simulations to account for the uncertainties in our standard assumptions. We find that
a significant fraction, 15
+9
−8
%, of core-collapse supernovae are “late”, that is, they occur 50–200 Myr after birth, when all massive
single stars have already exploded. These late events originate predominantly from binary systems with at least one, or, in most cases,
with both stars initially being of intermediate mass (4–8 M
⊙
). The main evolutionary channels that contribute often involve either
the merging of the initially more massive primary star with its companion or the engulfment of the remaining core of the primary by
the expanding secondary that has accreted mass at an earlier evolutionary stage. Also, the total number of core-collapse supernovae
increases by 14
+15
−14
% because of binarity for the same initial stellar mass. The high rate implies that we should have already observed
such late core-collapse supernovae, but have not recognized them as such. We argue that φ Persei is a likely progenitor and that
eccentric neutron star – white dwarf systems are likely descendants. Late events can help explain the discrepancy in the delay-time
distributions derived from supernova remnants in the Magellanic Clouds and extragalactic type Ia events, lowering the contribution of
prompt Ia events. We discuss ways to test these predictions and speculate on the implications for supernova feedback in simulations
of galaxy evolution.
Key words. supernovae: general – binaries: close – stars: massive – stars: evolution
1. Introduction
Core-collapse supernovae (ccSNe) are bright explosions that
mark the end of the lives of massive stars (e.g.,
Heger et al.
2003
; Smartt 2009) and the birth of neutron stars or black holes
(e.g., Ertl et al. 2016). They play a crucial role as sources of
chemical enrichment (e.g., Arnett 1973; Woosley et al. 2002)
and feedback, driving the evolution of their host galaxies (e.g.,
Hopkins et al. 2014). Their extreme brightness also allows us to
use them as probes of star-forming galaxies out to appreciable
redshifts (e.g., Strolger et al. 2015). These explosions are usually
attributed to stars with birth masses larger than approximately
8 M
⊙
(Heger et al. 2003), although the exact value depends
on model assumptions concerning core overshooting, stellar-
wind mass-loss, and metallicity (e.g., Poelarends et al. 2008;
Jones et al. 2013; Takahashi et al. 2013; Doherty et al. 2015).
Observing campaigns of young massive stars in our Galaxy
and the Magellanic Clouds show that a very large fraction
have one or more companions, forming a close binary system
where severe interaction between the stars during their lives
is unavoidable (e.g., Kobulnicky & Fryer 2007; Mason et al.
2009; Sana et al. 2012; Chini et al. 2012). Such interaction can
be the exchange of mass and angular momentum through
Roche-lobe overflow, common envelope evolution, and merg-
ing of the two stars (
Wellstein & Langer 1999; de Mink et al.
2013; De Marco & Izzard 2017). This interaction can drasti-
cally affect the further evolution of both stars and thus the
properties of their possible supernovae. Pioneers in modeling
the effects on (samples of) ccSNe include
Podsiadlowski et al.
(1992), De Donder & Vanbeveren (2003b), Yoon et al. (2010),
Eldridge et al. (2008, 2013).
Our understanding of the endpoints of massive stars is radi-
cally being transformed by the rise of (automated) transient sur-
veys, which enable the efficient detection of ccSNe and other
transients in large numbers. Examples are the Lick Observatory
Supernova Search (LOSS,
Filippenko et al. 2001), the Palomar
Transient Factory (PTF,
Rau et al. 2009; Law et al. 2009), and its
near-future upgrade, the Zwicky Transient Facility, the All-Sky
Automated Survey for SuperNovae (ASAS-SN, Shappee et al.
2014
), Pan-STARRS (Kaiser et al. 2002), and eventually the
Large Synoptic Survey Telescope (LSST, Ivezic et al. 2008). The
datasets provided by these facilities will be large, but may not
necessarily provide very detailed information about individual
Article published by EDP Sciences A29, page 1 of 22
A&A 601, A29 (2017)
events, since this typically requires more intensive follow up, to
obtain spectra, for example. The large potential of these datasets
will be the statistical constraints that they can provide, allowing
for new constraints on theoretical models for both common and
rare events. Fully harvesting these datasets will require adapta-
tions from the theory side and thus predictions of the statistical
properties for large samples will be needed.
Motivated by the technological and observational develop-
ments, as well as the insight into the large importance of binarity,
we have started a systematic theoretical investigation aiming to
quantify the impact of binarity on the statistical properties ex-
pected for large samples of ccSNe. This paper is the first in a
series in which we describe the motivation and setup of our sim-
ulations (Sect.
2). In two papers that were completed ahead of
this one, the lead authors of this team demonstrated the early ap-
plication of these new simulations against observations of two
individual events.
In
Van Dyk et al. (2016), we compared these simulations
with new deep Hubble Space Telescope observations of the site
of the now faded stripped-envelope type Ic supernova SN1994I
in search of a surviving companion star. While no companion
was detected, the data provided new strong upper brightness lim-
its, constraining the companion mass to less than 10 M
⊙
. This
result is consistent with the theoretical predictions of our simu-
lations and allowed a subset of formation scenarios to be ruled
out. In
Margutti et al. (2017), we used these simulations to inter-
pret the multi-wavelength observations of supernova SN2014C
which over the timescale of a year underwent a complete meta-
morphosis from an ordinary H-poor type Ib supernova into a
strongly interacting, H-rich supernova of type IIn. These simu-
lations helped us to estimate the possibility that the surround-
ing hydrogen shell originated from a prior binary interaction (as
opposed to ejection resulting from instabilities during very late
burning phases, e.g.,
Quataert et al. 2016).
In this paper, we focus on the distribution of the expected
delay time between formation of the progenitor star and its fi-
nal explosion, extending the work of
De Donder & Vanbeveren
(2003b). We investigate how binary interaction affects the delay-
time distribution of ccSNe. A significant fraction of ccSNe are
expected to be “late”, that is, they occur with delay times longer
than approximately 50 Myr, which is the maximum delay time
expected for single stars. We show that these late events origi-
nate from progenitors in binary systems with most of them be-
ing of intermediate mass. We discuss these late ccSNe in Sect. 3
and describe the various evolutionary channels that produce
them.
We further describe the outcome of an extensive study of the
robustness of our results against variations in the model assump-
tions and we compare with earlier work in Sect.
4. In Sect. 5,
we discuss (possible) observational evidence. We argue that the
well known binary φ Persei provides a direct progenitor sys-
tem that is expected to result in a late ccSN and we discuss
how the observed eccentric neutron star – white dwarf systems
may well provide the direct remnants. We then compare our re-
sults directly with the inferred delay time measured from su-
pernova remnants in the Magellanic Clouds, showing that they
are consistent with our predictions. We finish with a brief dis-
cussion on possible implications for feedback in star-forming
regions in galaxies by showing the differences with the widely
used single star predictions by the STARBURST99 simulations
(
Leitherer et al. 1999). We end with a summary of our findings
in Sect. 6.
2. Method
We use a binary population synthesis code, binary_c (ver-
sion 2.0, SVN revision 4105), developed by
Izzard et al. (2004,
2006, 2009) with updates described in de Mink et al. (2013) and
Schneider et al. (2015). The code employs rapid algorithms by
Tout et al. (1997) and Hurley et al. (2000, 2002) based on ana-
lytical fits to the detailed non-rotating single stellar models com-
puted by
Pols et al. (1998).
The code enables us to efficiently simulate the evolution of
single stars and binary systems from the zero-age main sequence
until they leave behind compact remnants. This allows us to
make predictions for an entire population of massive stars by
spanning the extensive parameter space of initial properties that
determine their evolution. It also allows us to explore the robust-
ness of our results against variations in our assumptions.
In Sect.
2.1, we discuss the initial conditions and in Sect. 2.2,
we discuss the physical assumptions. From now on when we
mention “standard models” or “standard simulations”, we refer
to the simulations where we followed our main assumptions in
all the key parameters that we discuss below. There are two stan-
dard models with one simulating only single stars and the other
including binaries (discussed also below in the paragraph for bi-
nary fraction). A summary of the key parameters, their values
for our standard assumptions and the model variations that we
consider is provided in
Table 1.
2.1. Initial conditions
Initial distributions. We assume that the distribution of the ini-
tial mass, M
1
, of primary stars (the initially most massive star
in a binary system) and of single stars follows a Kroupa (2001)
initial mass function (IMF),
dN
dM
1
∝ M
α
′
1
, (1)
where,
α
′
=
−0.3
−1.3
−2.3
α
0.01 < M
1
/ M
⊙
< 0.08,
0.08 < M
1
/ M
⊙
< 0.5,
0.5 < M
1
/ M
⊙
< 1,
1 < M
1
/ M
⊙
< 100.
(2)
In our standard models, we adopt α = −2.3. When assessing
the uncertainties, we consider variations in which α = −1.6 and
−3.0 following the uncertainty given in
Kroupa (2001) as well
as one model in which α = −2.7 (e.g., Kroupa et al. 1993).
For the initial mass ratio q ≡ M
2
/M
1
, where M
2
is the initial
mass of the secondary star, we take
dN
dq
∝ q
κ
. (3)
We adopt a uniform distribution in our standard simulation,
for example, κ = 0 for q ∈ [0.1, 1], consistent with
Kiminki & Kobulnicky (2012) and Sana et al. (2012). We also
consider the variations κ = −1 and 1.
For the initial orbital period distribution, we assume
dN
d log
10
P
∝
log
10
P
π
. (4)
We adopt π = 0, also known as
Öpik’s law (1924), for sys-
tems with primary masses up to 15 M
⊙
(Kobulnicky et al. 2014;
Moe & Di Stefano 2015). To account for the strong preference of
A29, page 2 of 22
E. Zapartas et al.: Delay-time distribution of core-collapse supernovae
Table 1. Summary of the key parameters adopted in our standard simulations and the variations that we consider.
Symbol Standard models
a
Model variations
Physical assumptions
- Mass transfer efficiency β β
th
0, 0.2, 1
- Angular momentum loss γ γ
orb,acc
0, γ
disk
- Mass loss during merger of two MS stars µ
loss
0.1 0, 0.25
- Mixing during merger of two MS stars µ
mix
0.1 0, 1
- Natal kick compact remnant (km s
−1
) σ 265 σ
0
, ∞
- Common envelope efficiency α
CE
1 0.1, 0.2, 0.5, 2, 5, 10
- Envelope binding energy λ
CE
λ
Dewi+00
0.5
- Critical mass ratio for contact for MS donor q
crit,MS
0.65 0.25, 0.8
- Critical mass ratio for unstable mass transfer for HG donor q
crit,HG
0.4 0, 0.25, 0.8, 1
- Stellar-wind mass-loss efficiency parameter η 1 0.33, 3
- Maximum single-star equivalent birth mass for ccSN (M
⊙
) M
max,cc
100 20, 35
- Minimum metal core for ccSN (M
⊙
) M
min,metal
1.37 1.3, 1.4
Initial conditions
- Slope initial mass function α −2.3 –1.6, –2.7, –3.0
- Slope initial mass ratio distribution κ 0 –1, 1
- Slope of initial period distr. π π
Opik24, Sana+12
–1, 1
- Metallicity Z 0.014 0.0002, 0.004, 0.008, 0.02, 0.03
- Binary fraction
a
f
bin
0.7, 0.0
a
0.3, 1, f
bin
(M
1
)
- Normalization parameter (M
⊙
) M
low
2 1, 3
Notes. See Sect. 2 for a description of the symbols and further assumptions.
(a)
The difference between our two standard models is that in one we
simulate only single stars and in the other we assume a binary fraction of 0.7.
more massive stars to reside in short period systems, we adopt
π = −0.55 when M
1
> 15 M
⊙
as found by
Sana et al. (2012). The
range of initial periods we consider is log
10
(P/day) ∈ [0.15, 3.5]
as given by
Sana et al. (2012). When assessing the uncertainties,
we consider π = −1 and 1 over the full mass range.
For the initial spin period of the stars, we follow
Hurley et al. (2000). Although this does not account for the
full distribution (e.g., Huang et al. 2010; Dufton et al. 2013;
Ramírez-Agudelo et al. 2013, 2015), this is sufficient for inves-
tigating the role of binarity as the impact of the adopted birth
spin is negligible compared to the angular momentum the star
later receives as a result of interaction by tides and mass transfer
(de Mink et al. 2013).
Although we account for the effects of eccentricity, we chose
to adopt circular orbits at birth to limit the number of dimensions
that our grid of models spans. This is justified as most systems
circularize shortly before the onset of mass transfer by Roche-
lobe overflow as a result of tides (Portegies Zwart & Verbunt
1996
; Hurley et al. 2002). However, with this approach, we do
not account for systems that are too wide to strongly interact
when circular, but where eccentricity implies periastron sepa-
rations small enough to trigger Roche-lobe overflow. We may
therefore slightly underestimate the impact of binary interaction.
See, for example, the interacting systems arising from binaries
with initial orbital periods well in excess of 10
3.5
days depicted
in Fig. 2 of de Mink & Belczynski (2015). We explore the un-
certainties arising from this assumption indirectly when we vary
the initial orbital period distribution and the total binary fraction.
Binary fraction. In our two standard models we either simulate
only single stars or we adopt a binary fraction of f
bin
= 0.7.
Here, we define a binary as a system with initial mass ratio q ∈
[0.1, 1] and initial period log
10
(P/day) ∈ [0.15, 3.5] based on
Sana et al. (2012) and consistent with the ranges adopted above.
We consider variations of f
bin
= 0.3 and 1.0.
The binary fraction for intermediate-mass stars is less well
constrained. We therefore consider a model variation where
we adopt a binary fraction that decreases with mass based on
Moe & Di Stefano (2013). These authors provide the inferred
fraction of systems with a companion in very close orbit, P = 2–
10 days, and q > 0.1. They find that this fraction drops from
0.22 for early B-type to 0.16 for late B-type stars. Information
on systems with orbits longer than 10 days are not available from
this study. These results may either indicate that the binary frac-
tion decreases towards late spectral types or that there is simply
a preference for systems with orbital periods longer than 10 days
in these stars. We use these results to construct a mass-dependent
binary fraction assuming that the binary companions of B-type
stars still follow an (
Öpik 1924) law over the full period range
0.15 < log
10
P < 3.5. We construct a mass-dependent binary
fraction f
bin
(M) referring to the binary fraction for the full period
range, such that the binary fraction for periods P = 2–10 days
are as in Moe & Di Stefano (2013). This results in
f
bin
(M
1
) =
0.44
0.61
0.7
M
1
/M
⊙
< 6
6 ≤ M
1
/M
⊙
< 15
15 ≤ M
1
/M
⊙
∼late B,
∼early B,
∼O,
(5)
which we adopt as one of the model variations. We want to stress
the importance of further observational campaigns aimed to con-
strain the initial binary distributions and the binary fraction for
the full M
1
, q and period range (e.g.,
Moe & Di Stefano 2017).
Normalization. When quoting absolute rates, we express our
results normalized by the total mass formed in stars in units of
10
6
M
⊙
. For this, we integrate over the full range of the IMF
A29, page 3 of
22
A&A 601, A29 (2017)
as specified in Eq. (
2). For stars with masses above M
low
, we
account for the mass contained in the companion star as specified
in Eq. (3). Effectively, we assume that low-mass stars, with M
1
<
M
low
, do not have companions massive enough to significantly
contribute to the mass of the stellar population. In our standard
assumptions we adopt M
low
= 2 M
⊙
and we vary this parameter
to 1 and 3 M
⊙
to check that this choice does not have a large
influence on our results.
Metallicity. We assume solar metallicity in our standard pop-
ulation, adopting a mass fraction of elements heavier than he-
lium of Z = 0.014 (Asplund et al. 2009) because present day
transient surveys focus on larger galaxies with metallicities that
are comparable to solar. We also consider low metallicities of
Z = 2 × 10
−4
relevant for metal-poor progenitors of globular
clusters and populations formed at high redshift. We addition-
ally test Z = 0.004 and Z = 0.008, relevant to nearby dwarf
galaxies similar to the Small and Large Magellanic Clouds. We
further provide results for the former canonical solar abundance
of Z = 0.02 (e.g.,
Grevesse et al. 1996) for comparison with ear-
lier studies and one super solar metallicity, Z = 0.03, relevant to
the central regions of large galaxies.
2.2. Physical assumptions
For a full description of the code, we refer to the references cited
at the start of this section. Here we discuss the main assumptions
that are of direct relevance to this study.
Stellar lifetimes. The evolutionary tracks and stellar lifetimes
(until reaching the white dwarf phase or core collapse) of single
stars in our simulations originate from the grid of detailed non-
rotating stellar evolutionary models of Pols et al. (1998) com-
puted with an updated version of the STARS code (Eggleton
1971, 1972; Pols et al. 1995). For stars up to 20 M
⊙
, we use
the fitting formulae of Hurley et al. (2000) for these models. At
higher mass, we switch to a logarithmic tabular interpolation of
the lifetimes by
Pols et al. (1998), and above 50 M
⊙
, we extrap-
olate as described in Schneider et al. (2015).
Our resulting mass-lifetime relation is shown in Fig. 1. We
find good agreement with simulations with the evolutionary code
MESA, version 7184 (
Paxton et al. 2011, 2013, 2015) for non-
rotating stars when using our standard metallicity, Z = 0.014,
and the Schwarzschild criterion for convection with a step-
overshooting parameter α
ov
= 0.335H
p
, where H
p
is the pressure
scale height, as calibrated by
Brott et al. (2011).
The widely-used Geneva models of
Schaller et al. (1992)
predict lifetimes that are 10–15% shorter, as can be seen in
Fig. 1. These are models of non-rotating stars with a metallicity
Z = 0.02. Our MESA simulations give similar lifetimes when
assuming no overshooting, α
ov
= 0. Given the evidence for extra
mixing processes beyond the convective core based on calibra-
tions of the overshooting parameter α
ov
(e.g., Ribas et al. 2000;
Claret 2007; Brott et al. 2011), we consider the lifetime predic-
tions by our models with overshooting to be more realistic.
Stellar winds. We include updated mass-loss prescriptions
as described in de Mink et al. (2013), which include the
recipes of
Vink et al. (2000). At luminosities in excess of
4000 L
⊙
, we switch to the empirical mass-loss rates of
Nieuwenhuijzen & de Jager (1990) when these rates exceed
those by Vink et al. (2000). To account for the empirical
10
0
10
1
10
2
10
3
10
4
Lifetime (Myr)
This work
MESA (α
ov
=0.335)
MESA (no overshooting)
Schaller+92
2 4 6 10 20 40 60 100
0.8
1.0
1.2
Initial mass (M
¤
)
Δ τ / τ
core-collapsewhite dwarfs
τ
max,cc
≃ 48 Myr
τ
min,cc
≃ 3 Myr
M
min,cc
≃ 7.5 M
¤
Fig. 1. Lifetime τ (until the white dwarf phase or core collapse) as a
function of initial mass for single stars adopted in this work (black line)
is compared with predictions that we obtained using the MESA stel-
lar evolutionary code (
Paxton et al. 2011) with and without overshoot-
ing (dark and lightblue dots) and with Geneva models of
Schaller et al.
(1992, black dots). Single stars with masses less than M
min,cc
end their
lives as white dwarfs instead of ccSNe (hashed region). ccSNe are ex-
pected between τ
min,cc
and τ
max,cc
(yellow shaded region), which refer to
the lifetimes of the most and least massive star to undergo core collapse.
The bottom panel shows the relative difference in lifetimes with respect
to the lifetimes used in our work (Sect.
2.2).
boundary of stars in the upper part of the Hertzsprung-Russell
diagram as described in Humphreys & Davidson (1994), we add
a factor in the mass loss as described in Hurley et al. (2000) to
simulate the enhanced mass loss of Luminous Blue Variables
(LBV) that are thought to reside near this boundary. For stars that
are stripped from their hydrogen envelopes, we adopt the Wolf-
Rayet (WR) mass-loss prescription by
Hamann et al. (1995) and
Hamann & Koesterke (1998) reduced by a factor of 10 to ac-
count for the effect of wind clumping (Yoon 2015). For post-
main-sequence stars, Asymptotic Giant Branch (AGB) stars,
and thermally pulsating AGB stars, we use Kudritzki & Reimers
(1978), Vassiliadis & Wood (1993), and Karakas et al. (2002) re-
spectively, as described in Izzard et al. (2009).
Our mass-loss prescriptions scale with metallicity as
˙
M ∝
(Z/Z
⊙
)
m
where m = 0.69 in main-sequence stars (
Vink et al.
2001; Mokiem et al. 2007). In post-main-sequence phases, we
adopt m = 0.5 (Kudritzki et al. 1989). In the WR phase, mass
loss scales with metallicity assuming a power-law index of 0.86
(
Vink & de Koter 2005). In the LBV phase, mass loss is assumed
to be invariant for metallicity.
The mass-loss rate by stellar winds as well as eruptive events
is uncertain, in particular for the late phases and the most mas-
sive stars (
Smith 2014). In the mass range we are most interested
in, mass loss during the late phases only affects the stellar enve-
lope. It does not have a large impact on the core of the stars
and thus on the remaining lifetimes, which are the main focus
of this work. Nevertheless, we explore the impact of changes in
the mass-loss rate by multiplying the mass loss by an efficiency
factor, η, which we set to unity in our standard simulations. We
consider the variations η = 3 and η = 0.33 for all mass-loss rates
simultaneously.
A29, page 4 of
22
E. Zapartas et al.: Delay-time distribution of core-collapse supernovae
Tides. We account for the effect of tides on the stellar spins
and the stellar orbits of stars in binary systems (
Zahn 1977;
Hurley et al. 2002) and the transfer of angular momentum dur-
ing mass transfer via an accretion disk or the direct impact of the
accretion stream onto the surface as described in de Mink et al.
(2013) following Ulrich & Burger (1976) and Packet (1981). We
assume that the stellar spins are aligned with the orbit (Hut
1981
).
Mass transfer. When a star fills its Roche lobe and mass trans-
fer is stable, we compute the mass-loss rate from the donor star
by removing as much mass as needed for the star to remain in-
side its Roche lobe. The resulting mass transfer rates are capped
by the thermal timescale of the donor. We define the mass trans-
fer efficiency, β, as the fraction of the mass lost by the donor that
is accreted by the companion,
β ≡
˙
M
acc
˙
M
don
· (6)
If mass is transferred on a timescale that is much shorter than
the thermal timescale of the accreting star τ
th,acc
, the star will be
driven out of thermal equilibrium and expand (
Neo et al. 1977).
Although our simulations do not follow this phase in detail, it is
expected that the companion can only accrete a fraction of the
transferred material when |
˙
M
don
| ≫
˙
M
acc,th
≡ M
acc
/τ
th,acc
, where
M
acc
is the mass of the accreting star. In line with this physical
picture, we limit the mass accretion rate to
|
˙
M
acc
| = min
|
˙
M
don
|, f
M
acc
τ
th,acc
!
, (7)
where f is an efficiency parameter for which we adopt 10 in our
standard simulation (to reproduce the mass transfer efficiency of
Schneider et al. 2015). The mass transfer efficiency in this case,
which we refer to as β
th
, varies between 0 and 1 depending on the
physical properties of the donor and the accretor. The efficiency
of mass transfer is poorly constrained (see, e.g., the discussion
in de Mink et al. 2007). We therefore also explore the extreme
case where none of the transferred mass is accreted, β = 0, the
case of very inefficient mass accretion, β = 0.2, as well as fully
conservative mass transfer, β = 1.
Angular momentum loss. Mass that is lost from the system
also takes away angular momentum. The specific angular mo-
mentum h, carried away from the system during mass loss, is
parametrized by,
h = γ
J
orb
M
acc
+ M
don
, (8)
where J
orb
is the total orbital angular momentum, M
acc
and M
don
are the masses of the accretor and donor star respectively and γ
is a free parameter. In our standard simulation, we assume that
mass lost from the system is emitted in a spherical wind or bipo-
lar outflow originating from the accreting star (
van den Heuvel
1994). Thus, the specific angular momentum, h, that the lost
mass carries is equal to the specific orbital angular momentum
of the accreting star, which yields γ = γ
orb,acc
≡ M
don
/M
acc
.
We also consider the extreme limiting case of negligi-
ble angular momentum transported by the mass lost from the
system during mass transfer (γ = 0). We further consider
the case where mass is lost through the outer Lagrangian
point, forming a circumbinary disk. Based on simulations by
Artymowicz & Lubow (1994), we consider that the binary sys-
tem will clear out the inner portion of the disk by resonance
torques. We explore the case that an inner region of size r
min
=
2a is cleared, where r
min
is the inner radius of the circumbinary
disk and a the separation of the binary system. This is consistent
with typical values found by Artymowicz & Lubow (1994). We
thus consider the case of γ = γ
disk
,
γ
disk
≡
(M
acc
+ M
don
)
2
M
acc
M
don
r
r
min
a
· (9)
Contact and common envelope evolution. To decide which
systems come into contact or experience common envelope (CE)
evolution, we consider a critical mass ratio, q
crit
. In binary sys-
tems with a post-main-sequence donor, we assume that systems
with M
acc
/M
don
< q
crit
enter a common envelope phase. We fol-
low the prescriptions of
Hurley et al. (2002) for q
crit
, except for
the case when the donor fills its Roche lobe while experiencing
hydrogen shell burning and crossing the Hertzsprung gap (HG),
where we use q
crit,HG
= 0.4 (de Mink et al. 2013). We use the
same value for the naked helium star donors that experience he-
lium shell burning (HeHG).
Common envelope evolution may either lead to the removal
of the envelope or, if the ejection is not successful, a merger. In
our treatment of common envelope evolution we use the formal-
ism described in
Tout et al. (1997) based on Webbink (1984),
Livio & Soker (1988) and de Kool (1990). In this formalism,
two parameters are introduced, the efficiency parameter of ejec-
tion, α
CE
, and λ
CE
which parametrizes the binding energy of the
envelope (see Eqs. (73) and (69) in Hurley et al. 2002, respec-
tively). In our standard model, we assume that α
CE
is unity (e.g.,
Webbink 1984; Iben & Tutukov 1984; Hurley et al. 2002), but
we also run models with a range of values (0.1, 0.2, 0.5, 2, 5,
10) to probe the large uncertainties associated with this phase
of evolution. By varying the efficiency parameter, we also im-
plicitly consider the effect of uncertainties in the binding en-
ergy λ
CE
as α
CE
λ
CE
appears as a product in the expression.
Values of α
CE
> 1 account for possible extra energy sources
used to unbind the envelope apart from the orbital energy (e.g.,
De Marco et al. 2011; Ivanova & Nandez 2016). To compute the
envelope binding energy parameter, λ
CE
, we use fits to detailed
models (
Dewi & Tauris 2000, 2001; Tauris & Dewi 2001). We
also consider a model variation where we adopt a constant value
λ
CE
= 0.5 (e.g., de Kool 1990). For a discussion of the limita-
tions of this formalism we refer to Ivanova et al. (2013b).
In systems with a main-sequence donor, we adopt q
crit,MS
=
0.65 to account for systems that come into contact during the
rapid thermal timescale mass transfer phase, which is consis-
tent with the detailed models by de Mink et al. (2007). Alterna-
tively, binary systems may come into contact because of their
own nuclear timescale evolution. In main-sequence stars, we as-
sume that contact leads to a merger.
Mergers and rejuvenation. In case a merger occurs, we fol-
low Table 2 of
Hurley et al. (2002) to determine the outcome.
When two main-sequence stars (MS+MS) merge, we follow
the updated algorithm by de Mink et al. (2013, 2014), based on
Glebbeek et al. (2013), to account for mass loss and internal
mixing. This algorithm uses two parameters, µ
loss
= 0.1, which
is the fraction of the total mass lost from the system during the
merger, and µ
mix
= 0.1, which is the fraction of the remaining
envelope mass that is mixed into the convective core. The values
A29, page 5 of
22