# The role of core-collapse physics in the observability of black-hole
neutron-star mergers as multi-messenger sources

Jaime Román-Garza Observatory of Geneva, University of Geneva, Chemin Pegasi
51, 1290, Versoix, Switzerland.  Simone S. Bavera Observatory of Geneva,
University of Geneva, Chemin Pegasi 51, 1290, Versoix, Switzerland.  Tassos
Fragos Observatory of Geneva, University of Geneva, Chemin Pegasi 51, 1290,
Versoix, Switzerland.  Emmanouil Zapartas Observatory of Geneva, University of
Geneva, Chemin Pegasi 51, 1290, Versoix, Switzerland.  Devina Misra
Observatory of Geneva, University of Geneva, Chemin Pegasi 51, 1290, Versoix,
Switzerland.  Jeff Andrews Center for Interdisciplinary Exploration and
Research in Astrophysics (CIERA) and Department of Physics and Astronomy,
Northwestern University, 1800 Sherman Avenue, Evanston, IL 60201, USA  Scotty
Coughlin Center for Interdisciplinary Exploration and Research in Astrophysics
(CIERA) and Department of Physics and Astronomy, Northwestern University, 1800
Sherman Avenue, Evanston, IL 60201, USA  Aaron Dotter Center for
Interdisciplinary Exploration and Research in Astrophysics (CIERA) and
Department of Physics and Astronomy, Northwestern University, 1800 Sherman
Avenue, Evanston, IL 60201, USA  Konstantinos Kovlakas Physics Department,
University of Crete, GR 71003, Heraklion, Greece  Institute of Astrophysics,
Foundation for Research and Technology-Hellas, GR 71110 Heraklion, Greece
Juan Gabriel Serra  Center for Interdisciplinary Exploration and Research in
Astrophysics (CIERA) and Department of Physics and Astronomy, Northwestern
University, 1800 Sherman Avenue, Evanston, IL 60201, USA  Ying Qin Center for
Interdisciplinary Exploration and Research in Astrophysics (CIERA) and
Department of Physics and Astronomy, Northwestern University, 1800 Sherman
Avenue, Evanston, IL 60201, USA  Department of Physics, Anhui Normal
University, Wuhu, Anhui 241000, China  Kyle A. Rocha Center for
Interdisciplinary Exploration and Research in Astrophysics (CIERA) and
Department of Physics and Astronomy, Northwestern University, 1800 Sherman
Avenue, Evanston, IL 60201, USA  Nam Hai Tran  DARK, Niels Bohr Institute,
University of Copenhagen, Jagtvej 128, 2200 Copenhagen, Denmark

(Received …; Revised …; Accepted …)

###### Abstract

Recent detailed 1D core-collapse simulations have brought new insights on the
final fate of massive stars, which are in contrast to commonly used parametric
prescriptions. In this work, we explore the implications of these results to
the formation of coalescing black-hole (BH) – neutron-star (NS) binaries, such
as the candidate event GW190426_152155 reported in GWTC-2. Furthermore, we
investigate the effects of natal kicks and the NS’s radius on the synthesis of
such systems and potential electromagnetic counterparts linked to them.
Synthetic models based on detailed core-collapse simulations result in an
increased merger detection rate of BH-NS systems ($\sim
2.3\,\mathrm{yr}^{-1}$), 5 to 10 times larger than the predictions of
“standard” parametric prescriptions. This is primarily due to the formation of
low-mass BH via direct collapse, and hence no natal kicks, favored by the
detailed simulations. The fraction of observed systems that will produce an
electromagnetic counterpart, with the detailed supernova engine, ranges from
$2$–$25$%, depending on uncertainties in the NS equation of state. Notably, in
most merging systems with electromagnetic counterparts, the NS is the first-
born compact object, as long as the NS’s radius is $\lesssim 12\,\mathrm{km}$.
Furthermore, core-collapse models that predict the formation of low-mass BHs
with negligible natal kicks, such as the ones from detailed core-collapse
studies, increase the detection rate of GW190426_152155-like events to $\sim
0.6\,\mathrm{yr}^{-1}$; with an associated probability of electromagnetic
counterpart $\leq 10$% for all supernova engines. However, increasing the
production of direct-collapse low-mass BHs also increases the synthesis of
binary BHs, over-predicting their measured local merger density rate. In all
cases, models based on detailed core-collapse simulation predict a ratio of
BH-NSs to binary BHs merger rate density that is at least twice as high as
other prescriptions.

Gravitational waves – Stars: black holes – Stars: neutron – Supernovae:
general – Stars: massive – Binaries

††journal: AJ††thanks<EMAIL_ADDRESS>
##  1 Introduction

The recently released catalogue of the LIGO Scientific and Virgo Collaboration
(LVC), GWTC-2, includes for the first time an event, GW190426_152155,
classified as a black hole (BH) – neutron star (NS) merger (Abbott et al.,
2020a). In addition, GW190814, an extreme mass-ratio merger event, has an
estimated mass of 2.59${}^{+0.08}_{-0.09}$ for the lower-mass compact object,
making it unclear whether it is a binary BH (BBH) or a BH–NS (BHNS) merger
(Abbott et al., 2020b). Due to the relatively low significance of
GW190426_152155 and the unclear nature of GW190814, no BHNS merger rate
density was estimated based on GWTC-2 (Abbott et al., 2020c), with the older
estimates from GWTC-1 setting only an upper limit of $<610\,\rm
Gpc^{3}\,yr^{-1}$ (Abbott et al., 2019). Nevertheless, taking into account
that the first half of the LVC’s third observing run (O3a) included 177.3 days
of data suitable for coincident analysis and assuming 1 or 2 detections, one
can estimate a detection rate of $\sim$2–4$\,\mathrm{yr}^{-1}$.

The first detection of a binary NS merger was accompanied by an
electromagnetic counterpart (EMC), which was observed in the whole
electromagnetic spectrum as a kilonova and a short Gamma-ray burst (Abbott et
al., 2017b, a). The merger of a BH with a NS is also expected to be
accompanied by a similar EMC, if the tidal disruption radius of the NS is
outside the innermost stable circular orbit (ISCO) of the BH. The maximum mass
of a non-spinning BH for this to happen, assuming a 1.4 M⊙ NS, is $\sim$3.5
M⊙, with the exact value depending on the adopted NS equation of state (e.g.
Capano et al., 2020). If, however, the BH is spinning, then the ISCO moves
closer to the BH and the corresponding BH mass limit becomes as high as $\sim
20\,\rm M_{\odot}$ for a maximally spinning BH (see Foucart et al., 2018),
vastly increasing the probability of an electromagnetic counterpart.

Under the assumption of efficient angular momentum transport in the interior
of stars, the first-born compact object in a binary is expected to be formed
with negligible spin (Fragos & McClintock, 2015; Qin et al., 2018). This is
because as the progenitor star expands to become a supergiant, most of its
angular momentum is transported to its outer layers, which are then removed
via winds and Roche-lobe overflow. On the other hand, the immediate progenitor
of the second-born compact object, in the isolated binary formation channels,
is a stripped helium (He) star in a close orbit with its first-born compact-
object companion. There, the He-star has a chance to be spun up via tides, and
thus give rise to a compact object with a significant spin (e.g., Van den
Heuvel & Yoon, 2007; Qin et al., 2018; Bavera et al., 2020a, b). Hence, the
only way to have a highly spinning BH in a BHNS system is to have the NS be
the first-born compact object, which will then tidally spin up the BH’s
progenitor star.

Coalescing BHNSs formed via isolated binary evolution are thought to be
sufficiently abundant, with theoretical estimates of their merger rate density
covering the whole range of 0.1–1000 Gpc3 yr-1 (Giacobbo & Mapelli, 2018;
Belczynski et al., 2020; Drozda et al., 2020). Among them, however, BHNSs in
which the NS is the first-born system should be rare for two main reasons: (i)
the initially more massive primary star is typically the progenitor of the BH
which tends to form first, and (ii) even if binary interactions reverse the
pre-core-collapse mass of the primary and secondary, kicks imparted on newly-
born NSs (Hobbs et al., 2005), which are typically higher than those on BHs,
tend to disrupt the binary when the system is in a wide orbit, as is typically
the case when the first compact object is formed.

The details of the general statements above depend crucially on the physics of
core-collapse and compact object formation. The vast majority of binary
population synthesis studies use the “rapid” and “delayed” mechanisms (Fryer
et al., 2012) to prescribe the fate of massive stars. Both of them are
parametric descriptions of the convection-enhanced supernova (SN) engine
driven by neutrino losses (see Herant et al., 1994). In contrast to the
“dealayed”, the “rapid” prescription predicts a mass gap between BHs and NSs
due to stronger convection which allows instabilities to grow rapidly and
produces more energetic SN explosions. Furthermore, both mechanisms predict a
unique boundary in the core mass of the pre-core-collapse star that leads to
the formation of a BH or a NS, above or below that boundary respectively.

Detailed 1D core-collapse simulations, on the other hand, show that there is
no unique boundary in the core mass of the pre-core-collapse star that
transitions between the formation of NSs and BHs. Instead, these studies find
successive islands of successful and failed explosions leading to the
formation of NSs and BHs via direct collapse (e.g., Pejcha & Thompson, 2015;
Sukhbold et al., 2016; Ertl et al., 2020; Patton & Sukhbold, 2020; Schneider
et al., 2020; Ugliano et al., 2012; O’Connor & Ott, 2011). This is the result
of the non-monotonic behaviour between the central carbon-burning phase and
the final core properties, linked to the convective episodes developed during
the burning phase (e.g. Sukhbold et al., 2018). Another significant difference
between these more recent calculations and the “rapid” and “delayed”, is in
the formation of BHs via successful explosions and significant fall-back.
Whereas in the “rapid” and “delayed” prescriptions there is a wide range of
pre-core-collapse core masses that lead to the formation of BHs via accretion
of fall-back mass from a successful explosion. Detailed 1D simulations find
that these cases are very rare and virtually all BHs are formed via direct
collapse. The latter becomes important in the context of population synthesis
studies, where the natal kicks imparted on BHs are most often normalized to
the fraction of fall-back mass, while BHs formed via direct collapse receive
no kick (e.g., Belczynski et al., 2008).

In this work we explore the effects of the core-collapse mechanism on the
formation of coalescing BHNSs, such as GW190426_152155, and their potential
observability as electromagnetic transients. For the first time we consider a
core-collapse prescription based on detailed 1D core-collapse simulations and
study the effect of non-monotonic stellar explodability, with respect to the
pre-collapse mass of the core. Finally, we discuss what we can learn in the
future in terms of formation pathways, core-collapse physics and NS equation
of state, once observations put a firmer constraints on the BHNS merger rate
density and the fraction of them accompanied by an EMC.

##  2 Methods

We use the software framework POSYDON (Fragos et al. 2021, in prep.) to evolve
populations of binaries for this study. POSYDON allows, among other
functionalities, to couple parametric population synthesis codes with such
models for different phases of a binary’s evolution. In this work, we coupled
the parametric code COSMIC (Breivik et al., 2019) to evolve binaries from the
zero-age main sequence (ZAMS) until the first compact object strips its
companion; and a grid of $\sim$170,000 detailed binary evolution models ran
with the MESA code (Paxton et al., 2011, 2013, 2015, 2018, 2019) to follow the
final evolutionary phase of a BHNS progenitor, i.e. that of a binary
consisting of stripped He-star and a compact object in a close orbit. This
allows us to accurately predict the spin of the second born compact object
(see Bavera et al., 2020b, for a detailed description of the simulation
setup).

For all populations evolved, we consider the following initial binary
properties: the mass of the most massive star $m_{1}$ is distributed by the
Kroupa (2001) initial mass function in a mass range of
$[5,150]\,\mathrm{M}_{\odot}$. We assume the mass ratio at birth is
distributed uniformly as $q\in[0,1]$ (Sana et al., 2012). The initial orbital
periods are distributed in the range $[0.4,10^{5.5}]$ days as in Sana et al.
(2012), extending the distribution for for low values with a flat distribution
as in Bavera et al. (2020b). All binaries have zero birth eccentricity, and we
assume an overall binary fraction $f_{b}=0.7$ (see Bavera et al., 2020b; Sana
et al., 2012). For each model, we evolve $5\times 10^{6}$ binaries per
metallicity, for 10 different metallicity values, $Z\in\ $[0.0001, 0.0002,
0.0003, 0.0006, 0.001 , 0.0018, 0.0031, 0.0055, 0.0098, 0.0174]. This
corresponds to a total stellar mass of $\sim 5\times 10^{8}\mathrm{M}_{\odot}$
for the underlying stellar population, per metallicity bin; and corresponds to
a fraction of the initial mass function of $f_{\mathrm{corr}}=0.212$.

We use the same set of physical model parameters as in Bavera et al. (2020b).
Specifically, we adopt mass-transfer (MT) stability according to the values of
$q_{\mathrm{crit}}$ as described there. Unstable MT is modelled with the
classical $\alpha_{\mathrm{CE}}-\lambda$ common-envelope (CE) formalism (van
den Heuvel, 1976; Webbink, 1984), where assume an $\alpha_{\mathrm{CE}}=1$ and
$\lambda_{\mathrm{CE}}$ fits from Claeys et al. (2014) without taking into
account the ionization energy of the envelope. Moreover, we assume the
pessimistic common-envelope scenario, namely all systems that start common
envelope evolution with a star in the Hertzsprung’s gap are considered to
merge due to the unsuccessful envelope ejection (Ivanova & Taam, 2004;
Belczynski et al., 2007). The electron-capture SN (ECSN) prescription
described in Podsiadlowski et al. (2004) is used, which maps Helium-core
masses in the range $[1.4,2.5]\,\mathrm{M}_{\odot}$ to remnant baryonic mass
$1.38\,\mathrm{M}_{\odot}$ as in Giacobbo & Mapelli (2020). For the pair-
instability and pulsational pair-instability SNe we consider the prescription
by Marchant et al. (2019) which limits the maximum BH mass at $\sim
44\,\mathrm{M}_{\odot}$.

To model the core-collapse, in addition to the “rapid” and “delayed”
mechanisms by Fryer et al. (2012), we implement a new prescription based on
detailed 1D core-collapse simulations. We use the publicly available data on
pre-SN models and the remnant properties produced by the N20 engine of
Sukhbold et al. (2016). We consider the He-core mass of the star at the pre-SN
phase, $M_{\mathrm{He,\,core}}$, to predict the remnant’s baryonic mass,
$M_{\mathrm{rem,\,bar}}$, taking into account whether the SN explosion is
predicted to be successful or not. We only consider the successful explosions
that will produce NSs, as the ones that produce BHs are rare. A complete
description of the implementation of the engine in our population synthesis
study is in the Appendix A.

For each mechanism three populations were produced with the assumptions
described in Table 1. Our fiducial model, called STANDARD, considers a
Maxwellian distribution with $\sigma_{\mathrm{ECSN}}=20$ km/s for kicks
imparted on NS formed from ECSN and $\sigma_{\mathrm{FeCCSN}}=265$ km/s for Fe
core-collapse SN (FeCCSN); the kicks in this model are fallback-weighted as in
Fryer et al. (2012). The FULL-ECSN-KICK model considers that NS formed from
ECSN and FeCCSN receive the same kicks
($\sigma_{\mathrm{ECSN}}=\sigma_{\mathrm{FeCCSN}}=265$ km/s). Finally, the NO-
BH-KICK model considers BHs receive no natal kicks, even if they are not
produced by direct collapse.

While the evolution of binaries from ZAMS to the formation of a compact object
plus a He-star was computed with the parametric code COSMIC, the last phase of
the evolution of a close compact object plus He-star was performed by
interpolating a grid constituted by 172,570 detailed MESA binary evolution
models. This last step, allows us to derive accurate predictions of the
second-born compact object’s spin (Qin et al., 2018; Bavera et al., 2020a).
The details on the grid and its interpolation are discussed in Appendix B.

We extract the BHNS mergers, systems that merge due to gravitational wave
radiation emission in less time than the current age of the Universe, and
compute the number of BHNS mergers per unit mass. The merging timescale, by
gravitational wave radiation, is computed according to Peters (1964). Adopting
the $\Lambda$CDM cosmology, we distribute the BHNS across the metallicity-
dependent cosmic star-formation history; assuming metallicities follow a
truncated log-normal distribution with standard deviation $0.5~{}\mathrm{dex}$
around the empirical mean metallicity function derived by Madau & Fragos
(2017). We compute the merger rate densities and detection rates as in Bavera
et al. (2020b), assuming the simulated “mid high/late low” LVC O3 detector
sensitivity (Abbott et al., 2018), by considering a single detector signal-to-
noise ratio threshold $>8$ that simulates a 2 network detector (see Barrett et
al., 2018). Finally, we compute the fraction of events that produce EMCs. We
assume that a BHNS merger will produce an electromagnetic counterpart if the
mass outside the BH ISCO after the merger, $M^{\mathrm{ejecta}}$, is greater
than zero. The detailed description on the computation of
$M^{\mathrm{ejecta}}$ is described in Appendix C. Taking into account
uncertainties in NS equation of state (e.g. Capano et al., 2020;
Chatziioannou, 2020, and references there in), we consider three different
constant values for the NS readius, $R_{\mathrm{NS}}\in[11,12,13]$ km, to
compute $M^{\mathrm{ejecta}}$ and, hence, predict the occurrence EMCs.

Table 1: Definition of the models used to evolve the binary stellar populations.  Model | $\sigma_{\mathrm{ECSN}}$ | $\sigma_{\mathrm{FeCCSN}}$ | NS | BH
---|---|---|---|---
name | [km/s] | [km/s] | kicks | kicks
STANDARD | 20 | 265 | fallback | fallback
|  |  | weighted | weighted
FULL-ECSN | 265 | 265 | fallback | fallback
-KICK |  |  | weighted | weighted
NO-BH | 20 | 265 | fallback | No
-KICK |  |  | weighted | kicks

Figure 1: BHNS mergers per unit mass in terms of relative metallicity to the
solar value (left) and its translation to the merger density rate history as a
function of redshift $z$ (right) for the STANDARD model. The dashed lines
represent the whole population of BHNSs mergers, while the contours enclose
the number of electromagnetic counterparts by considering the three values of
NS radius as described in Sec. 2; where upper and lower boundaries of the
contours are the predictions taking into account $R_{\mathrm{NS}}=13$km and
$R_{\mathrm{NS}}=11$km respectively. In the right panel, the gray vertical bar
is the $90\%$ confidence interval of local BHNS merger density rate reported
in GWTC-1.

##  3 Results

###  3.1 The theoretical merger rates of BHNS

The left panel of Fig. 1 shows the number of BHNS mergers per unit mass,
namely the total evolved mass corrected by the initial mass function and the
binary fraction, in terms of metallicity for the STANDARD model. In that panel
the dashed lines represent the whole population of BHNSs for each core-
collapse mechanism. While the shaded regions represent the subset of systems
that will reproduce an EMC. The upper boundary of a shaded region is defined
for $R_{\mathrm{NS}}=13$ km, while the lower boundary is delimited by
considering $R_{\mathrm{NS}}=11$ km. By distributing the mergers on the cosmic
star-formation history we find the merger density history in terms of
redshift, shown in the right panel of Fig. 1. There, we see that the local
merger density for all core-collapse mechanisms is consistent with the upper
limit reported on GWTC-1 (see Abbott et al., 2019) and plotted as the gray
shaded region. Furthermore, these results are consistent with other recent
population synthesis studies (e.g. Giacobbo & Mapelli, 2018; Neijssel et al.,
2019; Drozda et al., 2020; Belczynski et al., 2020).

We show the detection rate of BHNSs for simulated O3 LIGO/Virgo sensitivity
for all the populations in the upper panels of Fig. 2. The estimated BHNSs
detection rates are plotted with unfilled diamonds, while the filled symbols
represent the subset of systems that will produce an electromagnetic
counterpart for different values of $R_{\mathrm{NS}}$. All the rates shown on
the upper panels of Fig. 2 are also summarized in Table 2. The most striking
feature of Fig. 2 is that detection rates from the N20 engine (the purple
diamonds) are higher by a factor of $\sim 2-10$ with respect to the
predictions for the “rapid” and “delayed” mechanisms in all models. This is
because the N20 engine predicts the formation of low mass BHs by direct
collapse. These BHs do not receive natal kicks, that would otherwise disrupt
the binary, and they are produced by less massive stars whose number is
favored by the initial mass function.

To explore further the role of BH natal kicks on the detection rates of BHNS
mergers, we focus on the NO-BH-KICK model (upper right panel in the Fig. 2).
The assumption of no BH kicks increases the detection rate for the “delayed”
mechanism by a factor of ${\sim}4$ with respect to the STANDARD model. In
contrast the rate for the “rapid” prescription does not have a significant
increase between models. This difference is due to the fact that the pre-core-
collapse He-core mass, $M_{\mathrm{He,\,core}}$, range where the “rapid”
mechanism predicts BH formation via partial fall-back and non-zero kick
velocities is smaller ($M_{\mathrm{He,\,core}}\in[9,13]\mathrm{M}_{\odot}$)
compared to the “delayed” mechanism
($M_{\mathrm{He,\,core}}\in[6,13]\mathrm{M}_{\odot}$), as shown in Fig. 3.
Furthermore, since this region is located at higher $M_{\mathrm{He,\,core}}$
values for the “rapid” mechanism, the steepness of the initial mass function
results, in any case, in less BHs with non-zero kicks. Despite all of this,
the N20 engine produces more than twice as many BHNS mergers, even when only
considering the NO-BH-KICK models. This is explained by the fact that low mass
BHs formed by the “delayed” mechanism in the NO-BH-KICK model are not produced
by direct collapse. Meaning that only a fraction (as low as $\sim 30\%$) of
the pre-SN mass of their progenitors will collapse and the orbit of the binary
needs to be readjusted even if the BH does not receive a natal kick. The
readjustment of the orbit in such cases reduces the number of BHNS mergers
with respect the amount expected if the low mass BHS were produced by direct
collapse.

###  3.2 Rates of electromagnetic counterparts linked to BHNS mergers

The predicted rates of EMCs considering three different values for NS radii,
11, 12 and 13 km, are shown in the upper panels of Fig. 2 with the triangle,
square and circle markers, respectively. The predicted EMC rates for the N20
engine in the STANDARD population remains higher than the ones predicted using
the other core-collapse prescriptions and, for all NS radii. Note, however,
that it is the “delayed” mechanism that has the largest fraction of BHNS
mergers with EMCs. This can be understood from the fact that among the three
mechanisms considered, it is only the “dalayed” that can produce low-mass BHs,
in the range of $2.5-3.5\rm\,M_{\odot}$. These low-mass BHs, even when they
are non-spinning, are able to disrupt a relatively compact ($\lesssim 12\,\rm
km$) NS outside their ISCO.

For the “rapid” and N20 engines and assuming NS radii $\lesssim 12\,\rm km$,
BHNS mergers with EMCs are only produced when the NS is the first-born
compact-object which subsequently tidally spins up the BH progenitor star and
produces a higly spinning BH. The bottom row of Fig. 2 shows that although the
fraction of first-born NS in BHNS merges, for the “rapid” and N20 engines, is
$\sim 10\%$, the fraction of BHNS mergers with EMC that had a first-born NS is
close to 100%.

The formation of the first compact object occurs in a wide orbit. Since the
probability of the binary to remain bound after a SN kick scales with the
ratio of the orbital velocity over the kick velocity (Kalogera, 1996), a NS
formed via FeCCSN, that typically receives a natal kick of hundreds of km s-1,
will disrupt the binary. With the N20 engine, however, as BHs may be produced
by $M_{\mathrm{He,\,core}}$ as low as $\sim 4.5\rm\,M_{\odot}$, a non-
negligible fraction of first-born NSs may come from progenitors with low
enough pre-core-collapse $M_{\mathrm{He,\,core}}$ that results to an ECSN. The
low kicks associated with ECSN increase significantly the survivability of
these relatively rare binary configurations.

The FULL-ECSN-KICK models allow us to quantitatively explore the role of ECSN
on the formation of BHNS with EMCs, which as explained above affects models
with the N20 engine. Increasing the ECSN kick velocities has primarily an
impact on the systems with first-born NSs rather than the whole population of
BHNSs. From Table 2, we see that the BHNS merger detection rate slightly
decrease for N20 from STANDARD to FULL-ECSN-KICK. Contrarily, such rates for
the other core-collapse mechanisms remain similar, as the fraction of first-
born NSs from ECSN is negligible in those cases111In fact, we do see a very
small increase in the detection rate for those models, which is however above
the Poisson error of our simulations. This very small increase may be
explained by the fact that in those models some of the second-born NSs are
produced by ECSN. Imparting a larger natal kick on those NSs will not disrupt
the binary as its orbit is very close at that point in time, but will impart
some eccentricity in the post-SN orbit which may shorten the time to merge due
to GW emission. This in turn can have a small effect on the overall rate..
Most importantly though, high ECSN kicks decrease the rate of BHNS with EMCs,
in the N20 engine by a factor of $\sim 4$, for NS with radii $\lesssim 12\,\rm
km$, see Table 2.

Figure 2: Predicted merger rate of BHNS populations for O3a (upper panels) and
the fraction of systems with first-born NSs $f_{\mathrm{first\ born\ NS}}$
(lower panels). This predicted rates are plotted each core-collapse mechanism
mechanisms and each model respectively, see labels. In the upper panels, the
unfilled diamonds indicate the predicted detection rate for the whole
population of BHNSs, while the filled markers indicate the EMCs considering
different NS radii, see legend above. In the lower panels, the unfilled
hexagons denote the detected fraction of BHNSs with first-born NSs with
respect to the whole population of detected mergers, while the filled markers
are the fractions of EMC from systems with firstborn NSs with respect to all
of them, see legend above for the shape description.

Current constraints on the NS equation of state indicate that NS radii can be
as high as $13$ km (see Capano et al., 2020; Chatziioannou, 2020, and
references therein). This parameter has a crucial impact on the population of
EMCs from BHNS mergers where the BH is the first-born compact object. In such
cases the BH spin is negligible; therefore a $R_{\mathrm{NS}}$ on the high end
of the allowed parameter space is crucial to achieve a NS tidal disruption
(Foucart et al., 2018). In fact, for $R_{\mathrm{NS}}=13\,\rm km$ the maximum
BH mass that can tidally disrupt a NS is below the minimum BH mass predicted
by the N20 engine. This allows for BHNS mergers with EMCs originating from
binaries with non-spinning, first-born BHs and translates to an increase by a
factor of 15 and 2.5 of EMCs, compared to assumed $R_{\mathrm{NS}}$ of 11 km
and 12 km, respectively.

Table 2: The predicted local detection rates, $\mathcal{R}^{\mathrm{det}}$ (detected mergers per year), for the LVC’s third observing run. Results here are for the whole population of BHNSs (All) the electromagnetic counterparts (EMC), and for systems like GW190426_152155 (GW190426-like), and the subset of GW190425-like systems those systems that will produce an electromagnetic counterpart (GW190426 EMC) for each physical model.  |  |  | STANDARD |  | FULL-ECSN-KICK |  | NO-BH-KICK
---|---|---|---|---|---|---|---
Population | $R_{\mathrm{NS}}$ |  | $\mathcal{R}^{\mathrm{det}}_{\mathrm{N20}}$ | $\mathcal{R}^{\mathrm{det}}_{\mathrm{Rapid}}$ | $\mathcal{R}^{\mathrm{det}}_{\mathrm{Delayed}}$ |  | $\mathcal{R}^{\mathrm{det}}_{\mathrm{N20}}$ | $\mathcal{R}^{\mathrm{det}}_{\mathrm{Rapid}}$ | $\mathcal{R}^{\mathrm{det}}_{\mathrm{Delayed}}$ |  | $\mathcal{R}^{\mathrm{det}}_{\mathrm{N20}}$ | $\mathcal{R}^{\mathrm{det}}_{\mathrm{Rapid}}$ | $\mathcal{R}^{\mathrm{det}}_{\mathrm{Delayed}}$
All |  |  | 2.37 | 0.56 | 0.24 |  | 2.25 | 0.57 | 0.25 |  | 2.37 | 0.85 | 1.10
EMC |  $11$km |  | 0.04 | 0.01 | 0.01 |  | 0.01 | 0.01 | 0.01 |  | 0.04 | 0.01 | 0.05
$12$km |  | 0.23 | 0.01 | 0.03 |  | 0.06 | 0.01 | 0.03 |  | 0.23 | 0.01 | 0.13
$13$km |  | 0.6 | 0.01 | 0.04 |  | 0.42 | 0.01 | 0.05 |  | 0.6 | 0.01 | 0.21
GW190426-like |  |  | 0.56 | 0.03 | 0.06 |  | 0.55 | 0.05 | 0.07 |  | 0.56 | 0.12 | 0.42
GW190426 EMC |  $11$km |  | 0.00 | 0.00 | 0.00 |  | 0.00 | 0.00 | 0.00 |  | 0.00 | 0.00 | 0.00
$12$km |  | 0.01 | 0.00 | 0.00 |  | 0.00 | 0.00 | 0.00 |  | 0.01 | 0.00 | 0.00
$13$km |  | 0.05 | 0.00 | 0.02 |  | 0.03 | 0.00 | 0.01 |  | 0.05 | 0.00 | 0.03

###  3.3 The role of the SN engine on the synthesis of GW190426_152155-like
events

The estimation for the BHNS merger rate form the LVC O3a run, $2$ – $4$ yr-1,
is closer to the results from the N20 engine on all populations, as well as
for the “delayed” for the NO-BH-KICK population, favoring the engines that
predict low-mass BHs with negligible SN kicks. We compute the merger rate of
systems like GW190426_152155, the only observed system labeled as a BHNS
merger (Abbott et al., 2020a), by considering the $90\%$ confidence intervals
of the event’s measured total mass, chirp mass and the effective spin. Those
results are shown in Table 2 as the population “GW190426-like”. The detection
of a system similar to GW190426_152155 is favored by the N20 engine in all
populations and by the “delayed” mechanism in the NO-BH-KICK populations;
while the “rapid” prescription predicts a rate of only $0.03$ – $0.12$ yr-1,
indicating that the detection of an event with GW190426_152155-like properties
is rare in such case.

We also computed the rate of GW190426_152155-like systems that will produce an
EMC (population “GW190426 EMC” in Table 2). In all cases the rates of EMCs
linked to events like GW190426_152155 are less than $0.1$ yr-1, translating to
a probability of $0$ – $25$ % of EMCs per observed system. This low rate
agrees with the lack of an electromagnetic signal linked to the observed
event.

###  3.4 Comparison to BBH merger density rates

In order to compare our models with measurements with better constraints, we
compute the local BBH and BHNS merger density rates as shown in Table 3 of
Appendix C. The N20 STANDARD model, which increases the likelihood of events
like GW190426_152155, over-predicts the measured BBH merger density, as
reported by Abbott et al. (2020a), by $\sim 5$ times compared to the other
collapse mechanisms which over-predict the rate by $<2$ times. Again, this is
because the N20 engine assumes BH are born without kicks. Similarly, “rapid”
and “delayed” NO-BH-KICK over-predict the local BBH rate by the same amount.
When looking carefully at the BBH’s progenitor evolutionary pathway we find
that it is the CE channel that produces the majority of these merging BBHs
contributing to the observed local rate density. If, indeed, BHs are born
without a kick then the models are overproducing the systems going trough and
surviving CE and therefore the population of BBH, as we expect the mixture of
all formation channels to contribute on the whole population, which is what
recent studies suggests (e.g. Kruckow et al., 2016; Pavlovskii et al., 2017;
Klencki et al., 2020).

##  4 Discussion and conclusions

In this work, we present population synthesis models exploring the role of
core-collapse prescriptions, including results of detailed 1D core-collapse
simulations, and the associated natal kicks, on the observability of BHNS
systems, events like GW190426_152155, and their EMCs. A parametric code has
been used to model binary stars from ZAMS until the first-born compact object
strips its companion star. From that stage, the evolution of the system until
the formation of the second-born compact object was followed by detailed
binary-evolution models, allowing for accurate predictions of the spin of the
second-born compact object. The latter is critical in determining whether the
NS will be tidally disrupted by the BH, as spinning BHs increase significantly
the probability of EMCs produced by BHNS mergers.

We find that the N20 SN engine predicts BHNS merger rates higher by an order
of magnitude compare to the “rapid” and “delayed” mechanisms. This is a
consequence of the formation of low-mass BHs by direct collapse that are
predicted by SN engines from detailed core-collapse simulations, such as N20.
In addition, the N20 engine predicts higher rates of EMCs while being
consistent with the lack of observations of such events to date.

Our models show that future, more stringent constraints on the NS equation of
state will allow to distinguish between formation sub-channels of BHNSs. A
mean NS radius closer to $11$ km would indicate that the information from
future observations of BHNS EMCs is linked to systems where the NS is the
first-born compact object; as compact NSs are harder to disrupt by non-
spinning BHs. In contrast, evidence of a larger radii (such as Abbott et al.,
2020; Riley et al., 2019) will help us to infer information of low-mass BHs
and possible natal kicks linked to their formation.

Finally, we find that the synthesis of events like GW190426_152155 is favored
by SN engines that produce low-mass BHs with small or no SN kicks. The N20
engine predicts a rate of $\sim 0.6$ yr-1 for such event, one order of
magnitude higher than the predictions by the “rapid” and “delayed” mechanisms
in the STANDARD model, and 30% larger than the result from “delayed” with the
NO-BH-KICK model. However, such mechanisms overestimate the BBH local merger
density rate, with most of the predicted merging BBHs going through the CE
evolution channel. Detailed binary evolution calculations suggest that the
source of this apparent discrepancy may be the parametrizations of mass-
transfer stability criteria and envelope binding energy estimates, which, as
implemented in most binary population synthesis codes may, severely over-
predict the number of BBH progenitor systems going through and surviving the
CE phase.

The authors thank Vicky Kalogera and Christopher Berry for their thoughtful
comments. This work was supported by the Swiss National Science Foundation
Professorship grant (project number PP00P2 176868; PI Tassos Fragos). JRG is
supported by UNIGE, JJA and SC are supported by CIERA and AD, JGSP, and KAR
are supported by the Gordon and Betty Moore Foundation through grant GBMF8477.
KK received funding from the  European Research Council under the European
Union’s Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement n.
617001. YQ acknowledges funding from the Swiss National Science Foundation
under grant P2GEP2_188242. The computations were performed in part at the
University of Geneva on the Baobab and Lesta computer clusters and at
Northwestern University on the Trident computer cluster (the latter funded by
the GBMF8477 grant). All figures were made with the free Python modules
Matplotlib (Hunter, 2007). This research made use of
Astropy,222http://www.astropy.org a community-developed core Python package
for Astronomy (Robitaille et al., 2013; Price-Whelan et al., 2018).

## References

  * Abbott et al. (2020) Abbott, B., Abbott, R., Abbott, T., et al. 2020, The Astrophysical Journal Letters, 892, L3
  * Abbott et al. (2019) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, Physical Review X, 9, 031040, doi: 10.1103/PhysRevX.9.031040
  * Abbott et al. (2017a) —. 2017a, ApJ, 848, L13, doi: 10.3847/2041-8213/aa920c
  * Abbott et al. (2017b) Abbott, B. P., Abbott, R., Abbott, T. D., Acernese, F., & Ackley, K. e. a. 2017b, ApJ, 848, L12, doi: 10.3847/2041-8213/aa91c9
  * Abbott et al. (2018) Abbott, B. P., Abbott, R., Abbott, T., et al. 2018, Living Reviews in Relativity, 21, 3
  * Abbott et al. (2020a) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020a, arXiv e-prints, arXiv:2010.14527.  https://arxiv.org/abs/2010.14527
  * Abbott et al. (2020b) —. 2020b, ApJ, 896, L44, doi: 10.3847/2041-8213/ab960f
  * Abbott et al. (2020c) —. 2020c, arXiv e-prints, arXiv:2010.14533.  https://arxiv.org/abs/2010.14533
  * Barrett et al. (2018) Barrett, J. W., Gaebel, S. M., Neijssel, C. J., et al. 2018, Monthly Notices of the Royal Astronomical Society, 477, 4685
  * Bavera et al. (2020a) Bavera, S. S., Fragos, T., Qin, Y., et al. 2020a, A&A, 635, A97, doi: 10.1051/0004-6361/201936204
  * Bavera et al. (2020b) Bavera, S. S., Fragos, T., Zevin, M., et al. 2020b, arXiv e-prints, arXiv:2010.16333.  https://arxiv.org/abs/2010.16333
  * Belczynski et al. (2008) Belczynski, K., Kalogera, V., Rasio, F. A., et al. 2008, ApJS, 174, 223, doi: 10.1086/521026
  * Belczynski et al. (2007) Belczynski, K., Taam, R. E., Kalogera, V., Rasio, F. A., & Bulik, T. 2007, ApJ, 662, 504, doi: 10.1086/513562
  * Belczynski et al. (2020) Belczynski, K., Klencki, J., Fields, C. E., et al. 2020, A&A, 636, A104, doi: 10.1051/0004-6361/201936528
  * Breivik et al. (2019) Breivik, K., Coughlin, S. C., Zevin, M., et al. 2019, arXiv preprint arXiv:1911.00903
  * Capano et al. (2020) Capano, C. D., Tews, I., Brown, S. M., et al. 2020, Nature Astronomy, 1
  * Chatziioannou (2020) Chatziioannou, K. 2020, arXiv preprint arXiv:2006.03168
  * Claeys et al. (2014) Claeys, J., Pols, O., Izzard, R., Vink, J., & Verbunt, F. 2014, Astronomy & Astrophysics, 563, A83
  * Drozda et al. (2020) Drozda, P., Belczynski, K., O’Shaughnessy, R., Bulik, T., & Fryer, C. 2020, arXiv preprint arXiv:2009.06655
  * Ertl et al. (2020) Ertl, T., Woosley, S. E., Sukhbold, T., & Janka, H.-T. 2020, The Astrophysical Journal, 890, 51
  * Foucart et al. (2018) Foucart, F., Hinderer, T., & Nissanke, S. 2018, Physical Review D, 98, 081501
  * Fragos & McClintock (2015) Fragos, T., & McClintock, J. E. 2015, ApJ, 800, 17, doi: 10.1088/0004-637X/800/1/17
  * Fryer et al. (2012) Fryer, C. L., Belczynski, K., Wiktorowicz, G., et al. 2012, The Astrophysical Journal, 749, 91
  * Giacobbo & Mapelli (2018) Giacobbo, N., & Mapelli, M. 2018, Monthly Notices of the Royal Astronomical Society, 480, 2011
  * Giacobbo & Mapelli (2020) —. 2020, The Astrophysical Journal, 891, 141
  * Herant et al. (1994) Herant, M., Benz, W., Hix, W. R., Fryer, C. L., & Colgate, S. A. 1994, The Astrophysical Journal, 435, 339
  * Hobbs et al. (2005) Hobbs, G., Lorimer, D., Lyne, A., & Kramer, M. 2005, Monthly Notices of the Royal Astronomical Society, 360, 974
  * Hunter (2007) Hunter, J. D. 2007, Computing in science & engineering, 9, 90
  * Ivanova & Taam (2004) Ivanova, N., & Taam, R. E. 2004, ApJ, 601, 1058, doi: 10.1086/380561
  * Kalogera (1996) Kalogera, V. 1996, The Astrophysical Journal, 471, 352
  * Klencki et al. (2020) Klencki, J., Nelemans, G., Istrate, A. G., & Chruslinska, M. 2020, arXiv e-prints, arXiv:2006.11286.  https://arxiv.org/abs/2006.11286
  * Kroupa (2001) Kroupa, P. 2001, Monthly Notices of the Royal Astronomical Society, 322, 231
  * Kruckow et al. (2016) Kruckow, M. U., Tauris, T. M., Langer, N., et al. 2016, A&A, 596, A58, doi: 10.1051/0004-6361/201629420
  * Madau & Fragos (2017) Madau, P., & Fragos, T. 2017, The Astrophysical Journal, 840, 39
  * Marchant et al. (2019) Marchant, P., Renzo, M., Farmer, R., et al. 2019, The Astrophysical Journal, 882, 36
  * Neijssel et al. (2019) Neijssel, C. J., Vigna-Gómez, A., Stevenson, S., et al. 2019, MNRAS, 490, 3740, doi: 10.1093/mnras/stz2840
  * O’Connor & Ott (2011) O’Connor, E., & Ott, C. D. 2011, The Astrophysical Journal, 730, 70
  * Patton & Sukhbold (2020) Patton, R. A., & Sukhbold, T. 2020, arXiv preprint arXiv:2005.03055
  * Pavlovskii et al. (2017) Pavlovskii, K., Ivanova, N., Belczynski, K., & Van, K. X. 2017, MNRAS, 465, 2092, doi: 10.1093/mnras/stw2786
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3, doi: 10.1088/0067-0049/192/1/3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4, doi: 10.1088/0067-0049/208/1/4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15, doi: 10.1088/0067-0049/220/1/15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34, doi: 10.3847/1538-4365/aaa5a8
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10, doi: 10.3847/1538-4365/ab2241
  * Pejcha & Thompson (2015) Pejcha, O., & Thompson, T. A. 2015, The Astrophysical Journal, 801, 90
  * Peters (1964) Peters, P. C. 1964, Physical Review, 136, B1224
  * Podsiadlowski et al. (2004) Podsiadlowski, P., Langer, N., Poelarends, A., et al. 2004, The Astrophysical Journal, 612, 1044
  * Price-Whelan et al. (2018) Price-Whelan, A. M., Sipőcz, B., Günther, H., et al. 2018, The Astronomical Journal, 156, 123
  * Qin et al. (2018) Qin, Y., Fragos, T., Meynet, G., et al. 2018, Astronomy & Astrophysics, 616, A28
  * Riley et al. (2019) Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, The Astrophysical Journal Letters, 887, L21
  * Robitaille et al. (2013) Robitaille, T. P., Tollerud, E. J., Greenfield, P., et al. 2013, Astronomy & Astrophysics, 558, A33
  * Sana et al. (2012) Sana, H., De Mink, S., de Koter, A., et al. 2012, Science, 337, 444
  * Schneider et al. (2020) Schneider, F., Podsiadlowski, P., & Müller, B. 2020.  https://arxiv.org/abs/2008.08599
  * Sukhbold et al. (2016) Sukhbold, T., Ertl, T., Woosley, S., Brown, J. M., & Janka, H.-T. 2016, The Astrophysical Journal, 821, 38
  * Sukhbold et al. (2018) Sukhbold, T., Woosley, S. E., & Heger, A. 2018, ApJ, 860, 93, doi: 10.3847/1538-4357/aac2da
  * Ugliano et al. (2012) Ugliano, M., Janka, H., Arcones, A., & Marek, A. 2012, in Astronomical Society of the Pacific Conference Series, Vol. 453, Advances in Computational Astrophysics: Methods, Tools, and Outcome, ed. R. Capuzzo-Dolcetta, M. Limongi, & A. Tornambè, 91
  * Van den Heuvel & Yoon (2007) Van den Heuvel, E., & Yoon, S.-C. 2007, Astrophysics and Space Science, 311, 177
  * van den Heuvel (1976) van den Heuvel, E. P. J. 1976, in IAU Symposium, Vol. 73, Structure and Evolution of Close Binary Systems, ed. P. Eggleton, S. Mitton, & J. Whelan, 35
  * Webbink (1984) Webbink, R. F. 1984, ApJ, 277, 355, doi: 10.1086/161701
  * Zevin et al. (2020) Zevin, M., Spera, M., Berry, C. P., & Kalogera, V. 2020, The Astrophysical Journal Letters, 899, L1

##  Appendix A Implementation of the N20 engine

We consider the stellar models by Sukhbold et al. (2016), evolved in solar
metallicity, to describe the baryonic mass of the final remnant. We apply the
prescription deduced from their results to stars with different metallicity as
we aim to explore the effect of the average non-monotonic explodability trends
rather than predicting accurately the compactness of the core. Such non-
monotonic trend is preserved for different metallicities (Patton & Sukhbold,
2020).

Fig. 3 shows the baryonic remnant mass as a function of
$M_{\mathrm{He,\,core}}$, for the N20 engine (Sukhbold et al., 2016), the
“rapid” and “delayed” mechanisms (Fryer et al., 2012). On the same figure,
each stellar model, collapsed with N20, is classified as a successful
explosion or a direct collapse depending if the supernova shock is reactivated
by the neutrino flux or not, respectively. In the case of the N20 engine, to
predict if a star will undergo a successful explosion or direct collapse, we
extract the result in terms of the nearest-neighbor of the star’s
$M_{\mathrm{He,\,core}}$ with respect to the results from Sukhbold et al.
(2016). If the star is classified as a progenitor of a successful explosion,
its $M_{\mathrm{rem,\,bar}}$ will have the same value of the remnant baryonic
mass associated to the point with the nearest value for
$M_{\mathrm{He,\,core}}$. Otherwise, if the star is classified as a progenitor
of a direct collapse, then $M_{\mathrm{rem,\,bar}}$ is equal to the pre-SN
mass of the star.

To determine if the remnant is a BH or a NS, we calculate the remnant
gravitational mass considering the neutrino loss as in Zevin et al. (2020),
where, here, the maximum mass loss by neutrinos is considered to be
$0.5\,\mathrm{M}_{\odot}$. If the remnant gravitational mass is larger than
$2.5\,\mathrm{M}_{\odot}$ we assume that the compact object is a BH, else a
NS. For this work, the successful explosions that produce massive remnants
that will end up as BHs were not considered, as such cases are rare (from
Sukhbold et al., 2016, only 5 models from 105 successful explosions form a BH
for the N20 engine).

Figure 3: Baryonic remnant mass, $M_{\mathrm{rem,\,bar}}$, as a function of
the He-core mass at the pre-supernova phase, $M_{\mathrm{He,\,core}}$, from
N20 engine of Sukhbold et al. (2016) (circles); and for the “rapid” and
“delayed” (red and turquoise thick lines) mechanisms from Fryer et al. (2012).
For illustrative purposes we assume that the mass of the carbon-oxygen core is
$0.76M_{\mathrm{He,\,core}}$. Each model exploded with the N20 engine is
labeled as a successful explosion or direct collapse (i.e. failed explosion),
with purple and black colors, respectively. The black dashed line represents
the maximum NS baryonic mass, in our model, which corresponds to a maximum NS
gravitational mass of $2.5\,\mathrm{M}_{\odot}$.

##  Appendix B Grids of detailed NS/BH–He-star models

Figure 4: Two 2D slices of the 4D MESA grid for $Z=0.00312$, and
$m_{\mathrm{NS}}=1.44\,\mathrm{M}_{\odot}$ and
$m_{\mathrm{BH}}=8.57\,\mathrm{M}_{\odot}$, respectively. The final He-star
mass and angular momentum pre-SN values are colored for each successful track
according to each color bar. Each successful run stopped because of carbon
depletion or off-centre neon ignition (square markers), while other
termination flags are shown in the bottom legend. For visualisation purposes,
the models at $p=0.04$ days were excluded from the figure which for these
slices all runs resulted in ZAHeMS L1 or L2 Roche-lobe overflow.  Figure 5:
Median relative error as a function of metallicity of the log-transformed and
re-scaled five quantities $A_{i}$: the He-star mass (blue), angular momentum
(orange) and its carbon-oxygen core mass (green) before the supernova, the
orbital period before the supernova (red) and the lifetime of the binary
(green).

We use subset of detailed MESA models of BH–He-star systems of Bavera et al.
(2020b) which we extend to cover the NS star and low He-star masses. These
models treat the compact object as a point mass, hence, they can be applied to
simulate NS–He-star systems. For simplicity, we assumed that all physical
assumptions made in the BH–He-star regime applies also to the NS–He-star
regime too, including an Eddington mass-accretion rate limit of
$\dot{M}_{\mathrm{Edd}}=7.36\times
10^{-8}\left(M_{\mathrm{BH/NS}}/\mathrm{M}_{\odot}\right)\,\mathrm{M}_{\odot}{}\mathrm{yr}^{-1}$.
The selected subset of the original grid covers the initial parameter space of
10 metallicities, $Z$, in the log-range $[0.0001,0.0174]$ in log-steps of
$\Delta\log_{10}(Z)\simeq 0.25$, 11 BH masses in the log-range
$[2.5,54.4]\,\mathrm{M}_{\odot}$, 17 He-star masses in the log-range
$[8,80]\,\mathrm{M}_{\odot}$ and $20$ binary periods in the log-range
$[0.09,8]\,\mathrm{days}$. We extend this dataset to cover 10 NS masses in the
log-range $[1,2.28]\,\mathrm{M}_{\odot}$, 26 He-star masses where 20 are in
the log-range $[3,12]\,\mathrm{M}_{\odot}$, 5 in the range
$[3,7]\,\mathrm{M}_{\odot}$ (only for BH masses) and one at
$2.5\,\mathrm{M}_{\odot}$ and an extra period at 0.04 days. The smallest He-
star mass is chosen to guarantee coverage of the parameter space down to white
dwarf formation, while the maximum He-star mass and smallest orbital period
were chosen to include the full range of compact object – He-star systems
produced by our COSMIC models. The wide orbital period range ensures that we
cover the parameter space well past the point where any BHNS or BBH system
will merge within the age of the Universe. The original grid subset consisting
of 37,400 models was therefore extended to a total of 172,740 MESA models. The
fraction of failed MESA runs vary from $0.6\%$ to $1.5\%$ depending on
metallicity. In Figure 4 we show two 2D slices of the 4D parameter space
sliced at $Z=0.00312$ and, $m_{\mathrm{NS}}=1.44\,\mathrm{M}_{\odot}$ and
$m_{\mathrm{BH}}=8.57\,\mathrm{M}_{\odot}$, respectively, where we indicate
with a color the final He-star mass and angular momentum given initial orbital
separation and He-star masses of the detailed simulations.

These grids were used to determine the final outcomes and final parameters of
the late-end evolution stage of the binary systems through linear
interpolation. Each metallicity is interpolated separately. We want to
interpolate five quantities $A_{i}$: the He-star mass, angular momentum and
its carbon-oxygen core mass, orbital period before the supernova, as well as
the lifetime of the binary. Before interpolating each quantity, we log-
transformed it and re-scale it to the interval $[-1,1]$ to assign equal weight
to each dimension during the interpolation. The interpolation itself relies on
building a Delaunay triangulation of the input data points followed by
barycentric linear interpolation over the vertices of the (hyper)triangle
containing the location of interest. We test the accuracy of the interpolation
computing relative errors of a test grid which is constituted of an arbitrary
fraction (5%) of runs which we excluded from the train sample. To obtain a
consistent estimate, we repeat this experiment 10 times for each metallicity
and each interpolated quantity. If a point of an nonphysical region of the
parameter space (e.g. zero age He main sequence (ZAHeMS) overflow, max MT or
L2 overflow) is correctly interpolated to NaN (not a number) by the algorithm,
we consider it to have a zero relative error. On the other hand, if a point is
wrongly interpolated to NaN, we consider it to have a relative error of 1. In
Fig. 5 we report the median relative error of each transformed and rescaled
quantity $X_{i}\equiv\log_{10}(A_{i})^{[-1,1]}$ as
$\Delta_{i}=|X_{\mathrm{true,i}}-X_{\mathrm{interp,i}}|/X_{\mathrm{true,i}}$.
Because of the large sample of data-points we find small interpolation errors.
Most of the quantities show an increase of median relative error as a function
of metallicity. This is caused by the fact that at high metallicity the grids
show a less-linear behaviour than at low metallicity. This non-linearity is a
direct consequence of He-star stellar winds which, in our models, scale as
$(Z/Z_{\odot})^{0.85}$. In these systems, the He-stars lose a non-negligible
amount of mass and the orbits widen considerably. Moreover, NS–He-star in
tight orbits have higher probability to initiate a mass transfer case due to
the He-star tendency to expand more than at low metallicity.

##  Appendix C BHNS electromagnetic counterpart condition

We consider the model of Foucart et al. (2018) to determine the mass of the NS
that remains outside the BH ISCO after the tidal disruption, as

$\left[\max\left(\alpha\frac{1-2C_{\mathrm{NS}}}{\eta^{1/3}}{-}\beta\hat{R}_{\mathrm{ISCO}}\frac{C_{\mathrm{NS}}}{\eta}{+}\gamma,0\right)\right]^{\delta}M_{\mathrm{NS}}$
(C1)

where $\alpha{=}0.406$, $\beta{=}0.139$, $\gamma{=}0.255$, $\delta{=}1.761$,
$\hat{R}_{\mathrm{ISCO}}{\equiv}\frac{R_{\mathrm{ISCO}}c^{2}}{GM_{\mathrm{BH}}}$,
$C_{\mathrm{NS}}{=}Q\frac{R_{\mathrm{ISCO}}}{\hat{R}_{\mathrm{ISCO}}{R}_{\mathrm{NS}}}$,
$Q{=}\frac{M_{\mathrm{BH}}}{M_{\mathrm{NS}}}$, $\eta=Q/(1+Q)^{2}$. Here,
$M_{\mathrm{BH}}$ is the mass of the BH, $M_{\mathrm{NS}}$ the mass of the NS,
and $R_{\mathrm{NS}}$ is the radius of the NS. In this work we explore three
values for $R_{\mathrm{NS}}$, as 11, 12, 13, km.

Table 3:  Predicted local rate density $\mathcal{R}$ (in units of $\mathrm{Gpc}^{-1}{}\mathrm{yr}^{-1}$). Results here are shown for the common envelope (CE) and stable-mass transfer (SMT) channels sepratatly and combined for, both, the populations of BBH and BHNS for each model.  |  |  | STANDARD |  | FULL-ECSN-KICK |  | NO-BH-KICK
---|---|---|---|---|---|---|---
| Channel |  | $\mathcal{R}_{\mathrm{N20}}$ | $\mathcal{R}_{\mathrm{Rapid}}$ | $\mathcal{R}_{\mathrm{Delayed}}$ |  | $\mathcal{R}_{\mathrm{N20}}$ | $\mathcal{R}_{\mathrm{Rapid}}$ | $\mathcal{R}_{\mathrm{Delayed}}$ |  | $\mathcal{R}_{\mathrm{N20}}$ | $\mathcal{R}_{\mathrm{Rapid}}$ | $\mathcal{R}_{\mathrm{Delayed}}$
BBH | CE |  | 162.78 | 81.13 | 39.13 |  | 162.06 | 81.64 | 38.55 |  | 162.78 | 129.41 | 167.31
SMT |  | 39.73 | 33.05 | 31.27 |  | 41.43 | 33.6 | 31.08 |  | 39.73 | 31.38 | 26.36
CE + SMT |  | 202.51 | 114.18 | 70.4 |  | 203.49 | 115.24 | 69.63 |  | 202.51 | 160.79 | 193.67
BHNS | CE |  | 74.96 | 7.28 | 2.59 |  | 66.36 | 7.59 | 2.93 |  | 74.96 | 14.66 | 30.18
SMT |  | 2.28 | 2.86 | 3.06 |  | 2.99 | 3.12 | 2.88 |  | 2.28 | 1.88 | 2.49
CE + SMT |  | 77.24 | 10.14 | 5.65 |  | 69.35 | 10.71 | 5.81 |  | 77.24 | 16.54 | 32.67

##  Appendix D BBH and BHNS merger rate densities

To better understand the implication of the model physical modelling
assumptions done in this study, we also calculate the BBH merger rate
densities. The latter are much better constrained by GWTC-2, which finds the
local merger rate density of BBH to be
$23.9^{+14.9}_{-8.6}\,\mathrm{Gpc}^{-3}\mathrm{yr}^{-1}$ (Abbott et al.,
2020c). The BBH rate densities rising from the common envelope (CE) and stable
mass transfer (SMT) channels given the same set of assumption made here
(STANDARD-Delayed) are presented in Bavera et al. (2020b). For an one-to-one
comparison with BHNS rate densities of the different core-collapse and kick
prescriptions considered in this work, we summarise the rate densities of BBHs
and BHNS for the CE and SMT channels as well as their combination in Table 3.

In contrast to the original study done by Bavera et al. (2020b) the new models
differ in the following ways. First, we simulated a (i) metallicity range with
one third the resolution of Bavera et al. (2020b) but verified that this does
not have a noticeable impact on the rate estimates by reanalysing Bavera et
al. (2020b) models with the same metallicity sample resolution. Second, (ii)
the core-collapse of the secondary is assumed to be direct where in the
original study we followed the core-collapse of the MESA He-star profile at SN
accounting for disk formation. When an accretion disk is formed only a
fraction of its mass falls to the hole (see Appendix D in Bavera et al.,
2020b) which, in practice, means that here tidally spin up highly spinning BHs
are slightly more massive compared to Bavera et al. (2020b). Moreover, we only
(iii) interpolate binary properties before the supernova while in the original
work, which only investigated the delayed collapse mechanism, also
interpolated the second born compact object mass and spin (the former has, on
average, a larger interpolation error compared to the pre SN mass, see Fig.
E.1 of Bavera et al. (2020b)). Finally, we also updated the condition which
determines BH formation. We assume (iv) a BH is formed if the compact object
gravitational mass is larger than the maximum NS mass
($2.5\,\mathrm{M}_{\odot}$) while in the previous work we assumed that a
collapsing star leading to the formation of a BH had to have at least a
carbon-oxygen core mass and a remnant baryonic mass of $3\,\mathrm{M}_{\odot}$
in order to form a BH (as $0.5\,\mathrm{M}_{\odot}$ where assumed to be lost
because of neutrinos during the collapse of the proto-NS). We verified that
this change has no impact on the results.

In Table 3 we can see that the STANDARD-delayed model predicts a comparable
rate density $\sim{}39\,\mathrm{Gpc}^{-3}\mathrm{yr}^{-1}$ and
$\sim{}31\,\mathrm{Gpc}^{-3}\mathrm{yr}^{-1}$ for CE and SMT channels,
respectively. These values are in agreement with Bavera et al. (2020b) where
the small deviation in the numbers is given by the changes explained in the
previous paragraph. On the other hand, STANDARD-N20 overpredicts the CE+SMT
rates compared to the observations by a factor of at least 5. This is because
in the N20 engine all BHs are formed thorough direct collapse, without a kick.
When assuming no natal kicks (other than the readjustment of the orbits
because of neutrino mass loss) the NO-BH-KICK models with rapid and delayed
predict similar rates to the N20 engine meaning that the discrepancy is a
direct product of no BH kicks. In fact, if nature would to agree with the N20
engine, then the formation of merging BBH thorough the CE+SMT over-predicts
the systems surviving these channels. When looking more carefully at the rate
densities of these models we see that is the CE channel to over-predicting the
constrains from merging BBH in the local universe. As recent studies have
shown, the classical $\alpha_{\mathrm{CE}}-\lambda$ parameterization of CE and
mass transfer stability parameterization ($q_{\mathrm{crit}}$) might over-
predict the number of systems going thorugh and surviving this evolutionary
phase (e.g. Kruckow et al., 2016; Pavlovskii et al., 2017; Klencki et al.,
2020).


# Modelling Bitcoin Returns Volatility using GARCH models

Samuel Asante Gyamerah  Pan African University, Institute for Basic Sciences,
Technology, and Innovation, Kenya

###### Abstract

In this paper, an application of three GARCH-type models (sGARCH, iGARCH, and
tGARCH) with Student t-distribution, Generalized Error distribution (GED), and
Normal Inverse Gaussian (NIG) distribution are examined. The new development
allows for the modeling of volatility clustering effects, the leptokurtic and
the skewed distributions in the return series of Bitcoin. Comparative to the
two distributions, the normal inverse Gaussian distribution captured
adequately the fat tails and skewness in all the GARCH type models. The tGARCH
model was the best model as it described the asymmetric occurrence of shocks
in the Bitcoin market. That is, the response of investors to the same amount
of good and bad news are distinct. From the empirical results, it can be
concluded that tGARCH-NIG was the best model to estimate the volatility in the
return series of Bitcoin. Generally, it would be optimal to use the NIG
distribution in GARCH type models since time series of most cryptocurrency are
leptokurtic.

###### keywords:

Bitcoin , Volatility , tGARCH , Normal Inverse Gaussian

label1label1footnotetext: correspondence<EMAIL_ADDRESS>
##  1 Introduction

The price of Bitcoin are known to be very volatile. Market players and
investors are however interested in estimating accurately the volatility of
Bitcoin. This is as a result of the correlation between volatility and returns
on investment. It is notable that volatility are not directly observable and
consequently the need for efficient model that can capture the volatility in
the prices of Bitcoin. Estimating the volatility of Bitcoin is of high
importance since Bitcoin has the highest market capitalization in the
cryptocurrency market1. The cryptocurrency market is relatively new, hence
there have not been enough literatures on the different models used in
estimating the volatility in the market.
[1] modelled Bitcoin price data using an autoregressive jump-intensity GARCH
model and a standard GARCH model. Their results indicated that the
autoregressive jump-intensity GARCH model performed better in fitting the
Bitcoin price data than the standard GARCH model. [1] used multiple threshold-
GARCH and Asymmetric-power GARCH to measure and estimate the volatility in the
price of Bitcoin. They used Akaike information criterion (AIC), Bayesian
information criterion (BIC) and Hannan-Quinn information criterion (HQC) to
capture the leverage and regime switching features of the conditional variance
and concluded that Bitcoin market is not yet mature. Using the price of
Bitcoin, [2] studied the best conditional heteroscedasticity model in relation
to goodness-of-fit. They revealed that the optimal model was the the AR-cGARCH
model. [3] fitted 12 GARCH-type models to the log returns of the exchange
rates of seven different cryptocurrencies including Bitcoin. IGARCH (1, 1)
with normal innovations was found to be the optimal model using information
criterion. Their conclusion however differs from the known stylized
fact111returns have fatty tails and high kurtosis of financial time series
data.
Most financial time series data have shown that the conditional distribution
of returns series exhibit several stylized features such as excess kurtosis,
negative skewness, and temporal persistence in conditional movements. These
stylized facts and properties of the returns seris have significant
ramification for financial models especially volatility prediction, risk-
scenario distribution, and in the case of financial crisis. To capture these
stylized facts in these datas, econometricians have developed different tools
to model and estimate volatility.
This paper study the estimation capacity of different nonparametric GARCH-type
models on volatility of the return series of Bitcoin for the period between 01
January, 2014 to 16 August, 2019. Nonparametric methods are applied to the
GARCH-type models because they do not assume any distributional assumptions
and can capture the kurtosis and fatty tails of the return series of Bitcoin.
Hence, we employ the Students-t distribution, Generalized Error Distribution,
and the Normal Inverse Gaussian distribution.

##  2 Materials and Methodology

###  2.1 Data

In this study the closing price is selected as the price of a cryptocurrency
because it reflects all the activities of the cryptocurrency of the day.
Historical daily closing price of Bitcoin was extracted from
https://coinmarketcap.com/. The data extracted spanned from 01/01/2014 to
16/08/2019, totaling 2054 trading days. Bitcoin is a consensus network that
enables a new payment system and a completely digital currency.
Assume $P_{t}$ and $P_{t-1}$ represents the current day and previous day price
of a cryptocurrency, then the return series/log returms $(R_{t})$ is
calculated as

$R_{t}=log(P_{t})-log(P_{t-1})$ (1)

###  2.2 Test of Normality

The residuals data of the three cryptocurrency are checked for normality. If
the dataset is normally distributed, then a parametric statistic like the
Normal distribution can be assumed. However, if the residuals data are not
normally distributed, a non-parametric statistic will be used. If the
normality test fails, it is important to consider the histogram and the normal
probability plot to check the presence of outliers in the data set. We
employed the Jacque-Bera (JB) test statistics and the Anderson-Darling (AD)
test statistics for the normality test. The JB and the AD test statistics are
defined as in equation 2 and 3 respectively:

$JB=n\Big{(}\dfrac{S^{2}}{6}+\dfrac{(K-3)^{2}}{24}\Big{)},$ (2)

where $n$ is the sample size, $S$ and $K$ are the sample skewness and kurtosis
respectively.

$AD=-n-\dfrac{1}{n}\sum_{i=1}^{n}(2i-1)[InF(Y_{i})+In(1-F(Y_{n-i+1}))],$ (3)

where $F(\cdot)$ is the cumulative distribution function (CDF) for the
specified distribution and $i$ is the $i_{th}$ sample.
The JB test [4] uses the Lagrange multiplier approach on the Pearson family of
distributions to derive tests for normality if observations and regression
residuals. It is centered on sample skewness and kurtosis. The AD test [5] was
developed by Anderson and Darling in 1954. It gives more weight to the tails
of distribution. AD test is based on the empirical distribution function. The
JB and AD test compare the scores in the sample to a normally distributed set
of scores with same standard deviation and mean. The null hypothesis of these
test is "sample distribution is normally distributed". The distribution is not
normally distributed if the test is significant.

###  2.3 Testing for auto regressive conditional heteroscedasticity (ARCH)
Effects

To apply GARCH models to the Bitcoin returns series, we test for the presence
of stationarity and ARCH effects in the residual return series. The Ljung-box
and Lagrange multiplier (LM) test [6] are used to test for the presence of
ARCH effects in the data.

####  2.3.1 Ljung Box test

The Ljung Box test [7] is a method used in testing the absence of serial
autocorrelation up to a specified lag $k$. It can be defined as
$H_{0}$: The observed data do not exhibit serial correlation
$H_{1}$: The observed data exhibit serial correlation
The test statistics is given as:

$Q=N(N+2)\sum_{k=1}^{h}\dfrac{\hat{\rho}_{k}^{2}}{N-k}$ (4)

where $N$ is the sample size, $h$ is the time lag, $\hat{\rho}_{k}^{2}$ is the
estimated autocorrelation of the series at lag $k$.

####  2.3.2 Lagrange multiplier (LM) test

To apply GARCH models to the returns series data, it is important to test the
residuals for the presence of AutoRegressive Conditional Heteroscedasticity
(ARCH) effects. The LM test was used to test for the presence of ARCH effects.
LM is defined as;
$H_{0}$: No ARCH effect
$H_{1}$: Presence of ARCH effects.
The LM test statistics is given as:

$LM=N.R^{2}\sim\chi^{2}(q),$ (5)

where $N$ is the total number of observations, $R^{2}$ forms the regression,
and $q$ is the number of restrictions.

###  2.4 Conditional variance equation

####  2.4.1 Standard Generalized AutoRegressive Conditional Heteroskedasticity
(sGARCH)

GARCH is an extension of the ARCH model that integrates a moving average (MA)
part together with the autoregressive (AR) part.
Define a GARCH (p,q) model as

$\displaystyle x_{t}$ $\displaystyle=\sigma_{t}\epsilon_{t}$ (6)
$\displaystyle\sigma_{t}^{2}$
$\displaystyle=\omega+\sum_{i=1}^{p}\alpha_{i}x_{t-i}^{2}+\sum_{j=1}^{q}\beta_{j}\sigma_{t-j}^{2}$

where $\omega>0$, $\alpha_{i}>0$, $\beta_{j}>0$,
$\sum_{i=1}^{p}\alpha_{i}+\sum_{j=1}^{q}\beta_{j}<1$. $\epsilon_{t}$ is the is
a sequence of independent and identically distributed random variables. For
the sake of tractability, the order of all the GARCH models used will be
restricted to one.
The standard GARCH model [8] represented as sGARCH(1,1) is given as

$\sigma_{t}^{2}=\omega+\alpha_{1}\epsilon_{t-1}^{2}+\beta_{1}\sigma_{t-1}^{2}$

####  2.4.2 Integrated GARCH (iGARCH)

Suppose the sum of the coefficients in equation 6 is equal to 1, the GARCH
process is called an Integrated GARCH (iGARCH). That is, the iGARCH models are
unit-root GARCH models. A basic characteristic of the iGARCH model is that the
effect of past squared shocks $\eta_{t-i}=x_{t-1}^{2}-\sigma_{t-i}^{2}$ for
$i>0$ on $x_{t}^{2}$ is persistent. An iGARCH(p,q) model assumes the form

$\displaystyle x_{t}$ $\displaystyle=\sigma_{t}\epsilon_{t}$ (7)
$\displaystyle\sigma_{t}^{2}$
$\displaystyle=\omega+\sum_{i=1}^{p}\beta_{i}\sigma_{t-i}^{2}+\sum_{j=1}^{q}(1-\beta_{j})\alpha_{t-j}^{2}$

Define an iGARCH (1,1) model as

$\displaystyle x_{t}$ $\displaystyle=\sigma_{t}\epsilon_{t}$ (8)
$\displaystyle\sigma_{t}^{2}$
$\displaystyle=\omega+\beta_{1}\sigma_{t-1}^{2}+(1-\beta_{1})x_{t-1}^{2},$

where $1>\beta_{1}>0$.

####  2.4.3 Threshold GARCH (tGARCH)

The tGARCH model was developed independently by [9] and [10]. A tGARCH(p,q)
model assumes the form

$\sigma_{t}^{2}=\
+\sum_{i=1}^{q}(\alpha_{i}+\lambda_{i}\mathds{1}_{t-i})x_{t-i}^{2}+\sum_{j=1}^{p}\beta_{j}\sigma_{t-j}^{2}$
(9)

where

$\mathds{1}_{t-1}=\begin{cases}1,\qquad\mbox{if}\quad\epsilon_{t-i}<0,\qquad\mbox{bad
news}\\\ 0,\qquad\,\,\mbox{if}\quad\epsilon_{t-i}\geq 0,\qquad\mbox{good
news}\end{cases}$

$\alpha_{i}>0$, $\beta_{j}>0$, and $\lambda_{i}>0$.
The generalized version of the tGARCH TGARCH(1,1) is given as:

$\sigma_{t}^{2}=\omega+\alpha x_{t-1}^{2}+\beta\sigma_{t-1}^{2}+\lambda
x_{t-1}^{2}\mathds{1}_{t-1}$ (10)

where $\alpha$ and $\alpha+\lambda$ denote the effect of good news and bad
news respectively. A $\lambda>0$ is an evidence that bad news upsurge
volatility in the Bitcoin market. This indicates the existence of leverage
effects of the first order. For an asymmetric news effect, $\lambda\neq 0$.

###  2.5 Distribution Assumptions of the Error ($\epsilon_{t}$)

From the test of normality for the returns residuals in table 3 and figure 4,
it is clear that the distribution of the residuals returns are not normally
distributed. There is the presence of excess kurtosis and heavy-tails in the
distribution of the residuals returns. To account for the excess kurtosis and
fat tails that are present in the residuals of the returns series, we model
the error term in the GARCH models with a Student-t distribution, Generalized
Error Distribution (GED), and a Normal Inverse Gaussian (NIG) types of
distributions. These distributions are appropriate to capture the excess
kurtosis and the skewness in the residuals return series.

###  2.6 Student’s t-distribution

The Student’s t-distribution or the t-distribution is a sub-class of the
continuous probability distributions. It is used when the sample size is small
and the population standard deviation is unknown. It was proposed by [11]. The
probability density function (pdf) of the t-distribution is defined as

$f_{t}(x;v)=\dfrac{\Gamma(\frac{v+1}{2})}{\sqrt{v\pi}\Gamma(\frac{v}{2})}\Big{(}1+\dfrac{x^{2}}{v}\Big{)}^{-\frac{v+1}{2}},\qquad
x\in(-\infty,\infty)$ (11)

where $v$ is the number of degrees of freedom and $\Gamma(\cdot)$ is the gamma
function.
The log likelihood function is defined as below

$l(x;v)=N\,\mbox{In}\,\Bigg{(}\dfrac{\Gamma(\frac{v+1}{2})}{\sqrt{v\pi}\Gamma(\frac{v}{2})}\Bigg{)}-\dfrac{1}{2}\sum_{t=1}^{N}\Bigg{(}\,\mbox{In}\,\sigma_{t}^{2}+(v+1)\,\mbox{In}\,\Bigg{(}1+\dfrac{x_{t}^{2}}{\sigma_{t}^{2}(v-2)}\Bigg{)}\Bigg{]}$
(12)

####  2.6.1 Generalized Error Distribution (GED)

The standardized GED or the error distribution is a symmetrical unimodal sub-
class of the exponential family with shape parameter $v$. The pdf of the GED
is defined as in equation 13,

$f_{GED}(x;v)=k(v)\exp\left\\{-\dfrac{1}{2}\Big{|}\dfrac{x}{\lambda(v)}\Big{|}^{v}\right\\},\qquad
x\in(-\infty,\infty)$ (13)

where $v$ determines the tail weight with larger value of $v$ giving lesser
tai weight, $k(v)$, and $\lambda(v)$ are constant and are defined as

$\displaystyle k(v)$
$\displaystyle=\dfrac{v}{\lambda(v)2^{1+\frac{1}{v}}\Gamma\Big{(}\dfrac{1}{v}\Big{)}}$
$\displaystyle\lambda(v)$
$\displaystyle=\left\\{\dfrac{2^{-\frac{2}{v}\Gamma(\frac{1}{v})}}{\Gamma(\frac{3}{v})}\right\\}^{1/2}$

$\Gamma(x)=\int_{0}^{\infty}t^{x-1}e^{-t}dt,\,\,x>0$ is the Euler gamma
function.
The log likelihood function is defined as

$l(x;v)=\sum_{t=1}^{N}\Bigg{[}\mbox{In}\Big{(}\dfrac{v}{\lambda}\Big{)}-\dfrac{1}{2}\Big{|}\dfrac{x}{\sigma_{t}^{2}v}\big{|}^{v}-(1+v^{-1})\,\mbox{In}\,2-\mbox{In}\,\Gamma\Big{(}\dfrac{1}{v}\Big{)}-\dfrac{1}{2}\,\mbox{In}\,\sigma_{t}^{2}\Bigg{]}$
(14)

####  2.6.2 Normal Inverse Gaussian (NIG) Distribution

The NIG distribution was proposed by [12] and it is a sub-class of the
generalized hyperbolic distribution. The pdf of an NIG distribution for a
random variable $X$ is defined as

$f_{NIG}(x;\alpha,\kappa,\mu,\delta))=\dfrac{\alpha\delta}{\pi}\exp\\{\delta\sqrt{\alpha^{2}-\beta^{2}}+\beta(x-\mu)\\}K_{1}[\alpha
q(x)]$ (15)

where $q(x)=((x-\mu)^{2}+\delta^{2})^{1/2}$, $\alpha>0$ determines the shape,
$0\leq\mid\kappa\mid\leq\alpha$ determines the skewness, $\delta>0$ is the
scaling, $\mu\in\mathbb{R}$ is the location parameter and $x\in\mathbb{R}$.
$k_{1}(x)$ is the modified Bessel function of the second kind with index 1.
The NIG distribution has a heavier tail than the normal distribution and can
take different kinds of shapes. From equation 15, the shape of an NIG density
is described using a four dimensional parameters ($\alpha,\kappa,\mu,\delta$).
The NIG distribution is appropriate for capturing data sets with extremal
observations, skewness, and fat tails or semi-heavy tails. The log-likelihood
of the NIG distribution is given as

$l(x;\alpha,\kappa,\mu,\delta)=n\log\Big{(}\dfrac{\alpha\delta}{\pi}\Big{)}+n\delta\sqrt{\alpha^{2}-\kappa^{2}}+\kappa\sum_{t=1}^{N}(x_{t}-\mu)+\sum_{t=1}^{N}\log
K_{1}(x_{t};\alpha,\delta,\mu)$ (16)

###  2.7 Model selection

The best model was selected using two information criteria: the Akaike
information criteria (AIC), Bayesian information criteria (BIC). AIC and BIC
considers the accuracy of the model fit and the number of parameters in the
model; rewarding a better fit and penalizing an increased number of parameters
in the return series data. The optimal model that will be selected is the
model with the minimum AIC and BIC values. The GARCH models used in this study
were fitted using the maximum likelihood method.

$AIC=-2\,\,\mbox{In}\,\,\mbox{L}(\hat{\Theta})+2k,$
$BIC=-2\,\,\mbox{In}\,\,\mbox{L}+k\,\,\mbox{In}\,\,\mbox{L}(\hat{\Theta}),$

where $n,k,n,\Theta,\hat{\Theta}$ represent the number of observations, number
of unknown parameters, vector of the unknown parameters, and the maximum
likelihood estimates of the vector of the unknown parameters respectively.
The smaller the AIC and BIC values of a model, the better the fit of that
model.

##  3 Results and Discussion

###  3.1 Descriptive Statistics

The visualization in figure 1 revealed that the entire cryptocurrency market
is propped primarily by Bitcoin. Table 1 is the statistics of Bitcoin on the
cryptocurrency market and table 2 presents the descriptive statistics of the
return series of the three cryptocurrencies. The return series of Bitcoin is
positively skewed. These indicate that the returns of Bitcoin is non-
symmetric. The positive value of the skewness indicates that the distribution
of Bitcoin return series is skewed to the right or positively skewed. The
positive excess kurtosis ($216.7461$) of the Bitcoin returns series indicates
that the returns are leptokurtic. That is, the returns series has a fatty
tail. Figure 2 shows the time series plot of Bitcoin price (left figure) and
the return series (right figure) of Bitcoin for the time period. Figure 3 is
the histogram and the normal quantile-quantile (q-q) plot of the return series
for the same time period. The residuals plot is presented in figure 4. There
is the presence of fatty tails and skewness in the residuals of the return
series. This is confirmed from the Jacque-Bera and Anderson-Darling test
statistics is table 3. Both tests rejected the normality at $5\%$ significance
level.

Figure 1: Cryptocurrencies by market capitalization Table 1: Bitcoin
statistics, as at 01/09/2019

Cryptocurrency |  Return on investment (ROI) |  Market Capitalization |  All time high |  All time low
---|---|---|---|---
Bitcoin | 7,005.34% | $172,355,391, 698  | $20,089.00 | $65.53
Accessed at https://coinmarketcap.com/ on 01/09/2019 |

Table 2: Descriptive statistics of the three major cryptocurrency by market
capitalization

Cryptocurrency | Maximum | Minimum | Mean | Standard deviation | Skewness | Kurtosis
---|---|---|---|---|---|---
Bitcoin (USD) | -0.8488 | 1.4744 | 0.0012 | 0.0600 | 216.7471 | 6.4820

Table 3: Jacque-Bera and Anderson-Darling test of normality for the residuals of the returns series | Jacque-Bera test | Anderson-Darling test
---|---|---
P-value | $<2.2e-16$ | $<2.2e-16$
Figure 2: Closing price and return series of Bitcoin Figure 3: Histogram and
normal q-q plot of return series of Bitcoin Figure 4: Histogram and normal q-q
plot of return series of Bitcoin

The Augmented Dickey Fuller (ADF) test [13] is used to test for stationarity.
From table 4, the null hypothesis of statonarity is rejected at 5%
$\alpha$-level of significance. Hence, there is no need to difference the
return series. The Ljung-box and LM test are presented in table 4. From the
Ljung box test, the null hypothesis of ’no autocorrelation’ in the squared
residuals is rejected at 5% significance level. That is, there is dependency
in the squared returns series of Bitcoin. Using the Lagrange multiplier (LM)
test, the null hypothesis of ’no ARCH effects’ is rejected at 5% significance
level. From the Ljung box and LM test, it can be concluded that the volatility
ARCH effect is very much present in the return series. Hence, the GARCH models
to model the returns series data.

Table 4: Test of Auto Regressive Conditional Heteroscedasticity (ARCH) effect  | Ljung box test | LM test | ADF test
---|---|---|---
P-value | $7.835e-05$ | $<2.2e-16$ | 0.01

###  3.2 Estimated Volatility

Table 5 shows the maximum likelihood estimate (MLE) results of sGARCH(1,1),
iGARCH(1,1), and tGARCH(1,1) models for Bitcoin returns using student
t-distribution. From the table, the log-likelihood value $(4181.104)$ is
maximum for tGARCH(1,1) model. The values of the two information criterion
(AIC=$-4.0663$, BIC=$-4.0526$) of tGARCH(1,1) are minimum as compared to
sGARCH(1,1) and iGARCH(1,1). The visual QQ plot in Figure 5c is consistent
with the AIC, BIC, and log likelihood values of the tGARCH(1,1). These results
indicate that tGARCH(1,1) is the optimal model to describe the volatility of
the return series of Bitcoin using the student t-distribution. Table 6
presents MLE results of sGARCH(1,1), iGARCH(1,1), and tGARCH(1,1) models for
Bitcoin returns using the generalized error distribution. Compared to the
other models, the log-likelihood value is maximised under the tGARCH(1,1). The
AIC and BIC values of tGARCH(1,1) is the minimal as compared to the other two
models. Hence, the tGARCH(1,1) model is the best model relative to the
sGARCH(1,1) and iGARCH(1,1). From Table 7, the tGARCH(1,1) recorded the
maximum log likelihood value (4196.681). The AIC and BIC values (-4.0805 and
-4.0641 respectively) were also the lowest as compared to sGARCH(1,1) and
iGARCH(1,1). This indicates that tGARCH(1,1) is the best model for the
volatility of the return series of Bitcoin using the Normal Inverse Gaussian
distribution.
From the three selected models (tGARCH(1,1) using student t-distribution,
tGARCH(1,1) using generalized error distribution, and tGARCH(1,1) using normal
inverse Gaussian distribution), the tGARCH(1,1) using Normal Inverse Gaussian
distribution is the most optimal model for the volatility of Bitcoin return
series. This is confirmed from the maximum log likelihood value and the
minimum AIC and BIC values. This is clear from the visual plot in figure 8c.
It is therefore evident that optimal model in terms of the information
criterion, log likelihood, and QQ plot is the tGARCH(1,1) using the Normal
Inverse Gaussian distribution.
The volatility estimates obtained from tGARCH(1,1)-NIG model is displayed in
Figure 9a. It is evident that volatility moves through time. The density of
tGARCH(1,1)-NIG is shown in Figure 9b. Clearly, the NIG distribution was able
to capture the fat tails and skewness in the distribution. This confirms the
reliability and efficiency in using the NIG distribution for modelling the
volatility of Bitcoin return series.

Table 5: MLE results of sGARCH(1,1), iGARCH(1,1), and tGARCH(1,1) models for
Bitcoin returns with conditionally t-distributed errors

Model |  |  |  |  |  |
---|---|---|---|---|---|---
| $\hat{\omega}$ | $\hat{\alpha}$ | $\hat{\beta}$ | $1/v$ | $\hat{\lambda}$ |  | AIC | BIC | LogLikelihood
sGARCH(1,1) |  | 0.000036 | 0.192227 | 0.806773 | 3.200514 | — |  | -4.0367 | -4.0258 | 4149.72
|  | (0.000012) | (0.023060) | (0.023828) | (0.169435) | — |  |  |  |
iGARCH(1,1) |  | 0.000035 | 0.192695 | 0.807305 | 3.194742 | — |  | -4.0379 | -4.0297 | 4149.946
|  | (0.000010) | (0.021763) | — | (0.135971) | — |  |  |  |
tGARCH(1,1) |  | 0.001019 | 0.291106 | 0.835797 | 2.472141 | 0.007515 |  | -4.0663 | -4.0526 | 4181.104
|  | (0.000353) | (0.048010) | (0.018121) | (0.157392) | (0.059760) |  |  |  |
Note: Standard errors of estimates are reported in parentheses. |  |  |  |  |  |  |
Note 2: $\hat{\omega}$:The reaction of conditional variance, $\hat{\alpha}$:ARCH effect, $\hat{\beta}$:GARCH effect, $\hat{\lambda}$:Leverage effect |  |  |  |  |  |  |

Table 6: MLE results of sGARCH(1,1), iGARCH(1,1), and tGARCH(1,1) models for
Bitcoin returns with conditionally generalized error distribution errors

Model | Estimated parameters |  |
---|---|---|---
| $\hat{\omega}$ | $\hat{\alpha}$ | $\hat{\beta}$ | $1/v$ | $\hat{\lambda}$ |  | AIC | BIC | LogLikelihood
sGARCH(1,1) |  | 0.000039 | 0.192872 | 0.806128 | 0.860227 | — |  | -4.0495 | -4.0386 | 4162.87
|  | (0.000011) | (0.025780) | (0.022277) | (0.030944) | — |  |  |  |
iGARCH(1,1) |  | 0.000039 | 0.193739 | 0.806261 | 0.859610 | — |  | -4.0506 | -4.0424 | 4162.95
|  | (0.000010) | (0.022141) | — | (0.029442) | — |  |  |  |
tGARCH(1,1) |  | 0.001298 | 0.212642 | 0.826352 | 0.851479 | 0.055244 |  | -4.0632 | -4.0495 | 4177.926
|  | (0.000337) | (0.025368) | (0.021071) | (0.031225) | (0.061816) |  |  |  |
Note1: Standard errors of estimates are reported in parentheses. |  |  |  |  |  |  |

Table 7: MLE results of sGARCH(1,1), iGARCH(1,1), and tGARCH(1,1) models for
Bitcoin returns with conditionally normal inverse gaussian distribution errors

Model | Estimated parameters |  |
---|---|---|---
| $\hat{\omega}$ | $\hat{\alpha}$ | $\hat{\beta}$ | $1/v$ | $\hat{\kappa}$ | $\hat{\lambda}$ |  | AIC | BIC | LogLikelihood
sGARCH(1,1) |  | 0.00004 | 0.19748 | 0.80152 | 0.39161 | -0.13314 | — |  | -4.0599 | -4.0462 | 4174.485
|  | (0.000012) | (0.024326) | (0.022933) | (0.046759) | (0.032272) | — |  |  |  |
iGARCH(1,1) |  | 0.000039 | 0.198112 | 0.801888 | 0.390387 | -0.133277 | — |  | -4.061 | -4.050 | 4174.653
|  | (0.000010) | (0.022055) | — | (0.042123) | (0.032213) | — |  |  |  |
tGARCH(1,1) |  | 0.001048 | 0.233978 | 0.828983 | 0.319214 | -0.142017 | -0.005427 |  | -4.0805 | -4.0641 | 4196.681
|  | (0.000302) | (0.027732) | (0.019678) | (0.042593) | (0.033599) | (0.061556) |  |  |  |
Note: Standard errors of estimates are reported in parentheses. |  |  |  |  |  |  |  |

(a) sGARCH-Student t

(b) iGARCH-Student t

(c) tGARCH-Student t

Figure 5: QQ plot of sGARCH, iGARCH, and tGARCH using the Student
t-distribution (Student t)

(a) sGARCH-GED

(b) iGARCH-GED

(c) tGARCH-GED

Figure 6: QQ plot of the sGARCH, iGARCH, and tGARCH using the Generalized
Error Distribution (GED)

(a) sGARCH-NIG

(b) iGARCH-NIG

(c) TGARCH-NIG

Figure 7: QQ plot of sGARCH, iGARCH, and tGARCH using the Normal inverse
Gaussian (NIG) distribution

(a) TGARCH-t

(b) TGARCH-QED

(c) TGARCH-NIG

Figure 8: QQ plot using of the best distribution from sGARCH, iGARCH, and
TGARCH

(a) Conditional volatility

(b) Density

Figure 9: Conditional volatility and density of tGARCH(1,1)-NIG

##  4 Conclusion

This paper studied the volatility of daily return series of Bitcoin from
01/01/2014 to 16/08/2019. The results of the statistical properties revealed
that just like other financial time series data, the return series of Bitcoin
are leptokurtic. Different GARCH type models (sGARCH, iGARCH, tGARCH) were
compared and the tGARCH model was identified to be the most appropriate model
for estimating the time-varying volatility in Bitcoin return series. To
account for the skewness and fat tails, all the GARCH type models were
compared based on Student-t, Generalized Error Distribution (GED), and Normal
Inverse Gaussian (NIG) distribution. The NIG distribution performed better in
capturing the fat tail and skewness in the return series distribution. Hence,
the tGARCH-NIG model was the optimal model for modeling and estimating the
volatility in the return series of Bitcoin.

## Conflict of Interest

The author declare that there are no conflicts of interest regarding the
publication of this paper.

## Data Availability

Data for this work are available from the author upon request.

## References

  * Gronwald [2014] M. Gronwald,  The economics of bitcoins–market characteristics and price jumps (2014).
  * Katsiampa [2017] P. Katsiampa,  Volatility estimation for bitcoin: A comparison of garch models,  Economics Letters 158 (2017) 3–6.
  * Chu et al. [2017] J. Chu, S. Chan, S. Nadarajah, J. Osterrieder,  Garch modelling of cryptocurrencies,  Journal of Risk and Financial Management 10 (2017) 17.
  * Jarque and Bera [1987] C. M. Jarque, A. K. Bera,  A test for normality of observations and regression residuals,  International Statistical Review/Revue Internationale de Statistique (1987) 163–172.
  * Anderson and Darling [1954] T. W. Anderson, D. A. Darling,  A test of goodness of fit,  Journal of the American statistical association 49 (1954) 765–769.
  * Engle [2001] R. Engle,  Garch 101: The use of arch/garch models in applied econometrics,  Journal of economic perspectives 15 (2001) 157–168.
  * Box and Pierce [1970] G. E. Box, D. A. Pierce,  Distribution of residual autocorrelations in autoregressive-integrated moving average time series models,  Journal of the American statistical Association 65 (1970) 1509–1526.
  * Bollerslev [1986] T. Bollerslev,  Generalized autoregressive conditional heteroskedasticity,  Journal of econometrics 31 (1986) 307–327.
  * Ravichandran et al. [1989] K. Ravichandran, S. Bose, V. Akgiray, R. Al-Khouri, M. Ajlouni, M. Atmeh, I. Dobbs, G. Bekaert, G. Wu, F. Black, et al.,  Threshold generalized autoregressive conditional heteroskedasticity models.,  Research Journal of Business Management 6 (1989) 55–80.
  * Glosten et al. [1993] L. R. Glosten, R. Jagannathan, D. E. Runkle,  On the relation between the expected value and the volatility of the nominal excess return on stocks,  The journal of finance 48 (1993) 1779–1801.
  * Student [1908] Student,  The probable error of a mean,  Biometrika (1908) 1–25.
  * Barndorff-Nielsen [1977] O. Barndorff-Nielsen,  Exponentially decreasing distributions for the logarithm of particle size,  Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 353 (1977) 401–419.
  * Dickey and Fuller [1979] D. A. Dickey, W. A. Fuller,  Distribution of the estimators for autoregressive time series with a unit root,  Journal of the American statistical association 74 (1979) 427–431.