# How Much is 131 Million Dollars? Putting Numbers in Perspective with
Compositional Descriptions

Arun Tejasvi Chaganty
Computer Science Department
Stanford University
<EMAIL_ADDRESS>&Percy Liang
Computer Science Department
Stanford University
<EMAIL_ADDRESS>
###### Abstract

How much is 131 million US dollars? To help readers put such numbers in
context, we propose a new task of automatically generating short descriptions
known as perspectives, e.g. “$131 million is about the cost to employ everyone
in Texas over a lunch period”. First, we collect a dataset of numeric mentions
in news articles, where each mention is labeled with a set of rated
perspectives. We then propose a system to generate these descriptions
consisting of two steps: formula construction and description generation. In
construction, we compose formulae from numeric facts in a knowledge base and
rank the resulting formulas based on familiarity, numeric proximity and
semantic compatibility. In generation, we convert a formula into natural
language using a sequence-to-sequence recurrent neural network. Our system
obtains a $15.2\%$ F1 improvement over a non-compositional baseline at formula
construction and a 12.5 BLEU point improvement over a baseline description
generation.

##  1 Introduction

Figure 1:  An overview of the perspective generation task: given a numeric
mention, generate a short description (a perspective) that allows the reader
to appreciate the scale of the mentioned number. In our system, we first
construct a formula over facts in our knowledge base and then generate a
description of that formula.

When posed with a mention of a number, such as “Cristiano Ronaldo, the player
who Madrid acquired for […] a $131 million” (Figure 1), it is often difficult
to comprehend the scale of large (or small) absolute values like $131 million
(Paulos, 1988; Seife, 2010). Studies have shown that providing relative
comparisons, or _perspectives_ , such as “about the cost to employ everyone in
Texas over a lunch period” significantly improves comprehension when measured
in terms of memory retention or outlier detection (Barrio et al., 2016).

Previous work in the HCI community has relied on either manually generated
perspectives (Barrio et al., 2016) or present a fact as is from a knowledge
base (Chiacchieri, 2013). As a result, these approaches are limited to
contexts in which a relevant perspective already exists.

In this paper, we generate perspectives by composing facts from a knowledge
base. For example, we might describe $100,000 to be “about twice the median
income for a year”, and describe $5 million to be the “about how much the
average person makes over their lifetime”. Leveraging compositionality allows
us to achieve broad coverage of numbers from a relatively small collection of
familiar facts, e.g. median income and a person’s lifetime.

Using compositionality in perspectives is also concordant with our
understanding of how people learn to appreciate scale. Jones and Taylor (2009)
find that students learning to appreciate scale do so mainly by anchoring with
familiar concepts, e.g. $50,000 is slightly less than the median income in the
US, and by unitization, i.e. improvising a system of units that is more
relatable, e.g. using the Earth as a measure of mass when describing the mass
of Jupiter to be that of 97 Earths. Here, compositionality naturally unitizes
the constituent facts: in the examples above, money was unitized in terms of
median income, and time was unitized in a person’s lifetime. Unitization and
anchoring have also been proposed by Chevalier et al. (2013) as the basis of a
design methodology for constructing visual perspectives called concrete
scales.

When generating compositional perspectives, we must address two key
challenges: constructing familiar, relevant and meaningful formulas and
generating easy-to-understand descriptions or perspectives. We tackle the
first challenge using an overgenerate-and-rank paradigm, selecting formulas
using signals from familiarity, compositionality, numeric proximity and
semantic similarity. We treat the second problem of generation as a
translation problem and use a sequence-to-sequence recurrent neural network
(RNN) to generate perspectives from a formula.

We evaluate individual components of our system quantitatively on a dataset
collected using crowdsourcing. Our formula construction method improves on F1
over a non-compositional baseline by about 17.8%. Our generation method
improves over a simple baseline by 12.5 BLEU points.

##  2 Problem statement

The input to the _perspective generation_ task is a sentence $s$ containing a
numeric mention $x$: a span of tokens within the sentence which describes a
quantity with value $x.\mathrm{value}$ and of unit $x.\mathrm{unit}$. In
Figure 1, the numeric mention $x$ is “$131 million”, $x.\mathrm{value}=1.31e8$
and $x.\mathrm{unit}=\text{\$}$. The output is a description $y$ that puts $x$
in perspective.

We have access to a knowledge base $\mathcal{K}$ with numeric tuples
$t=(t.\mathrm{value},t.\mathrm{unit},t.\mathrm{description})$. Table 1 has a
few examples of tuples in our knowledge base. Units (e.g. $/per/yr) are
fractions composed either of fundamental units (length, area, volume, mass,
time) or of ordinal units (e.g. cars, people, etc.).

Description |  Value Unit
---|---
cost of an employee |  71e3 $/year/person
population of Texas |  27e3 person
number of employees at Google |  57e3 person
average household size |  2.54 person
time taken for a basketball game |  60 minute
average lifetime for a person |  79 year
a week |  1 week
time taken for lunch |  30 minute
cost of property in the Bay area |  1e3 $/ft2
area of a city block |  10e3 $m^{2}$
Table 1:  A subset of our knowledge base of numeric tuples. Tuples with
fractional units (e.g. $/ft2) can be combined with other tuples to create
formulas.

The first step of our task, described in Section 4, is to construct a formula
$f$ over numeric tuples in $\mathcal{K}$ that has the same value and unit as
the numeric mention $x$. A valid formula comprises of an arbitrary multiplier
$f.m$ and a sequence of tuples $f.\mathrm{tuples}$. The value of a formula,
$f.\mathrm{value}$, is simply the product of the multiplier and the values of
the tuples, and the unit of the formula, $f.\mathrm{unit}$, is the product of
the units of the tuples. In Figure 1, the formula has a multiplier of $1$ and
is composed of tuples 1, 2 and 3; it has a value of $1.3e8$ and a unit of $.

The second step of our task, described in Section 5, is to generate a
perspective $y$, a short noun phrase that realizes $f$. Typically, the
utterance will be formed using variations of the descriptions of the tuples in
$f.\mathrm{tuples}$.

##  3 Dataset construction

We break our data collection task into two steps, mirroring formula selection
and description generation: first, we collect descriptions of formulas
constructed exhaustively from our knowledge base (for generation), and then we
use these descriptions to collect preferences for perspectives (for
construction).

### Collecting the knowledge base.

We manually constructed a knowledge base with 142 tuples and 9 fundamental
units111Namely, length, area, volume, time, weight, money, people, cars and
guns. These units were chosen because they were well represented in the
corpus. from the United States Bureau of Statistics, the orders of magnitude
topic on Wikipedia and other Wikipedia pages. The facts chosen are somewhat
crude; for example, though “the cost of an employee” is a very context
dependent quantity, we take its value to be the median cost for an employer in
the United States, $71,000. Presenting facts at a coarse level of granularity
makes them more familiar to the general reader while still being appropriate
for perspective generation: the intention is to convey the right scale, not
necessarily the precise quantity.

### Collecting numeric mentions.

We collected 53,946 sentences containing numeric mentions from the newswire
section of LDC2011T07 using simple regular expression patterns like
$([0-9]+(,[0-9]+)*(.[0-9]+)? ((hundred)|(thousand)|(million)|
(billion)|(trillion))). The values and units of the numeric mentions in each
sentence were normalized and converted to fundamental units (e.g. from miles
to length). We then randomly selected up to 200 mentions of each of the 9
types in bins with boundaries
${10}^{-3},1,{10}^{3},{10}^{6},{10}^{9},{10}^{12}$ leading to 4,931 mentions
that are stratified by unit and magnitude.222Some types had fewer than 200
mentions for some bins. Finally, we chose mentions which could be described by
at least one numeric expression, resulting in the 2,041 mentions that we use
in our experiments (Figure 2). We note that there is a slight bias towards
mentions of money and people because these are more common in the news corpus.

Figure 2:  A histogram of the absolute values of numeric mentions by type.
There are 100–300 mentions of each unit.

### Generating formulas.

Next, we exhaustively generate valid formulas from our knowledge base. We
represent the knowledge base as a graph over units with vertices and edges
annotated with tuples (Figure 3). Every vertex in this graph is labeled with a
unit $u$ and contains the set of tuples with this unit:
$\\{t\in\mathcal{K}:t.\mathrm{unit}=u\\}$. Additionally, for every vertex in
the graph with a unit of the form $u_{1}/u_{2}$, where $u_{2}$ has no
denominator, we add an edge from $u_{1}/u_{2}$ to $u_{1}$, annotated with all
tuples of type $u_{2}$: in Figure 3 we add an edge from money/person to money
annotated with the three person tuples in Table 1. The set of formulas with
unit $u$ is obtained by enumerating all paths in the graph which terminate at
the vertex $u$. The multiplier of the formula is set so that the value of the
formula matches the value of the mention. For example, the formula in Figure 1
was constructed by traversing the graph from money/time/person to money: we
start with a tuple in money/time/person (cost of an employee) and then
multiply by a tuple with unit time (time for lunch) and then by unit person
(population of Texas), thus traversing two edges to arrive at money.

Using the 142 tuples in our knowledge base, we generate a total of 1,124
formulas sans multiplier.

Figure 3:  The graph over tuples generated from the knowledge base subset in
Table 1.

### Collecting descriptions of formulas.

The main goal of collecting descriptions of formulas is to train a language
generation system, though these descriptions will also be useful while
collecting training data for formula selection. For every unit in our
knowledge base and every value in the set
$\\{10^{-7},10^{-6}\ldots,10^{10}\\}$, we generated all valid formulas. We
further restricted this set to formulas with a multiplier between $1/100$ and
$100$, based on the rationale that human cognition of scale sharply drops
beyond an order of magnitude (Tretter et al., 2006). In total, 5000 formulas
were presented to crowdworkers on Amazon Mechanical Turk, with a prompt asking
them to rephrase the formula as an English expression (Figure
4).333Crowdworkers were paid $0.08 per description. We obtained 5–7
descriptions of each formula, leading to a total of 31,244 unique
descriptions.

Figure 4:  A screenshot of the crowdsourced task to generate natural language
descriptions, or perspectives, from formulas.

### Collecting data on formula preference.

Finally, given a numeric mention, we ask crowdworkers which perspectives from
the description dataset they prefer. Note that formulas generated for a
particular mention may differ in multiplier with a formula in the description
dataset. We thus relax our constraints on factual accuracy while collecting
this formula preference dataset: for each mention $x$, we choose a random
perspective from the description dataset described above corresponding to a
formula whose value is within a factor of 2 from the mention’s value,
$x.\mathrm{value}$. A smaller factor led to too many mentions without a valid
comparison, while a larger one led to blatant factual inaccuracies. The
perspectives were partitioned into sets of four and displayed to crowdworkers
along with a “None of the above” option with the following prompt: “We would
like you to pick up to two of these descriptions that are useful in
understanding the scale of the highlighted number” (Figure 5). A formula is
rated to be useful by simple majority.444Crowdworkers were paid $0.06 to vote
on each set of perspectives.

Figure 5:  A screenshot of the crowdsourced task to identify which formulas
are useful to crowdworkers in understanding the highlighted mentioned number.
Figure 6:  A histogram comparing formula length to ratings of usefulness
(clipped for readability). Non-compositional perspectives with a single tuple
are broadly useful. Useful compositional perspectives tend to be more context-
specific than non-compositional ones, and many of the formulas that can be
generated from the knowledge base are spurious.

Figure 6 provides a summary of the dataset collected, visualizing how many
formulas are useful, controlling for the size of the formula. The exhaustive
generation procedure produces a large number of spurious formulas like “20
$\times$ trash generated in the US $\times$ a minute $\times$ number of
employees on Medicare”. Nonetheless, compositional formulas are quite useful
in the appropriate context; Table 2 presents some mentions with highly rated
perspectives and formulas.

Sentence |  That’s about … |  Formula
---|---|---
The Billings-based Stillwater Mining produced 601,000 ounces of platinum. |  4 times the weight of an elephant. |  4 $\times$ weight of an elephant.
Authorities estimate there are about 60 million guns in Yemen. |  twice the gun ownership of the population of Texas |  2 $\times$ gun ownership $\times$ population of Texas
Water is flowing into Taihu lake at a rate of 150 cubic meters per second. |  how much water would flow from a tap left on for a week. |  rate of flow of water from tap $\times$ a week
The bank had held auctions, selling around US$1 billion worth of three-month bills. |  half the cost of employing the population of Texas for a work day. |  1/2 $\times$ cost of an employee $\times$ time taken for a work day $\times$ population of Texas
The government[s] have promised to rent about 1.2 million sq. feet. |  the area of forest logged in a single minute |  90 $\times$ area of forest logged $\times$ a minute
Table 2:  Examples of numeric mentions, perspectives and their corresponding
formulas in the dataset. All the examples except the last one are rated to be
useful by crowdworkers.

##  4 Formula selection

We now turn to the first half of our task: given a numeric mention $x$ and a
knowledge base $\mathcal{K}$, select a formula $f$ over $\mathcal{K}$ with the
same value and unit as the mention. It is easy to generate a very large number
of formulas for any mention. For the example, “Cristiano Ronaldo, the player
who Madrid acquired for […] $131 million.”, the small knowledge base in Table
1 can generate the 12 different formulas,555The full knowledge base described
in Section 3 can generate 242 formulas with the unit money (sans multiplier).
including the following:

  1. 1.

$1$ $\times$ the cost of an employee $\times$ the population of Texas $\times$
the time taken for lunch.

  2. 2.

$400$ $\times$ the cost of an employee $\times$ average household size
$\times$ a week.

  3. 3.

$1$ $\times$ the cost of an employee $\times$ number of employees at Google
$\times$ a week.

  4. 4.

$1$ $\times$ cost of property in the Bay Area $\times$ area of a city block.

Some of the formulas above are clearly worse than others: the key challenge is
picking a formula that will lead to a meaningful and relevant perspective.

### Criteria for ranking formulas.

Type |  Features | #
---|---|---
Proximity |  $\operatorname*{sign}(\log(f.m))$, $|\log(f.m)|$ | 1
Familiarity |  $\mathbb{I}[t]$ | 142
Compatibility |  $\mathbb{I}[t,t^{\prime}]$ | 20022
Similarity |  ${\operatorname{wvec}(s)}^{\top}$ $\operatorname{wvec}(t.\mathrm{description})$ | 1
Table 3:  Feature templates used to score a formulas $f$ and their counts (#),
where $f.m$ is the formula’s multiplier and $t,t^{\prime}\in
f.\mathrm{tuples}$ are tuples in the formula.

We posit the following principles to guide our choice in features (Table 3).

Proximity: A numeric perspective should be within an order of magnitude of the
mentioned value. Conception of scale quickly fails with quantities that exceed
“human scales” (Tretter et al., 2006): numbers that are significantly away
from $1/10$ and $10$. We use this principle to prune formulas with multipliers
not in the range $[1/100,100]$ (e.g. example 2 above) and introduce features
for numeric proximity.

Familiarity: A numeric perspective should be composed of concepts familiar to
the reader. The most common technique cited by those who do well at scale
cognition tests is reasoning in terms of familiar objects (Tretter et al.,
2006; Jones and Taylor, 2009; Chevalier et al., 2013). Intuitively, the
average American reader may not know exactly how many people are in Texas, but
is familiar enough with the quantity to effectively reason using Texas’
population as a unit. On the other hand, it is less likely that the same
reader is familiar with even the concept of Angola’s population.

Of course, because it is so personal, familiarity is difficult to capture.
With additional information about the reader, e.g. their location, it is
possible to personalize the chosen tuples (Kim et al., 2016). Without this
information, we back off to a global preference on tuples by using indicator
features for each tuple in the formula.

Compatibility: Similarly, some tuple combinations are more natural (“median
income $\times$ a month”) while others are less so (“weight of a person
$\times$ population of Texas”). We model compatibility between tuples in a
formula using an indicator feature.

Similarity: A numeric perspective should be relevant to the context. Apart
from helping with scale cognition, a perspective should also place the
mentioned quantity in appropriate context: for example, NASA’s budget of $17
billion could be described as 0.1% of the United States’ budget or the amount
of money it could cost to feed Los Angeles for a year. While both perspectives
are appropriate, the former is more relevant than the latter.

We model context relevance using word vector similarity between the tuples of
the formula and the sentence containing the mention as a proxy for semantic
similarity. Word vectors for a sentence or tuple description are computed by
taking the mean of the word vectors for every non-stop-word token. The word
vectors at the token level are computed using word2vec (Mikolov et al., 2013).

Formula |  Score
---|---
Studies estimate 36,000 people die on average each year from seasonal flu.
1/4 $\times$ global death rate $\times$ a day |  0.67
5 $\times$ death rate in the US $\times$ a day |  0.64
1/3 $\times$ number of employees at Microsoft |  0.60
Gazprom’s exports to Europe […] will total 60 billion cubic meters …
oil produced by the US $\times$ average lifetime |  0.78
average coffee consumption $\times$ population of the world $\times$ average lifetime |  0.78
2 $\times$ average coffee consumption $\times$ population of Asia $\times$ average lifetime |  0.73
Table 4:  The top three examples outputted by the ranking system with the
scores reported by the system.

### Evaluation.

We train a logistic regression classifier using the features described in
Table 3 using the perspective ratings collected in Section 3. Recall that the
formula for each perspective in the dataset is assigned a positive (“useful”)
label if it was labeled to be useful to the majority of the workers. Table
7(a) presents results on classifying formulas as useful with a feature
ablation.666Significance results are computed by the bootstrap test as
described in Berg-Kirkpatrick et al. (2012) using the output of classifiers
trained on the entire training set.

Familiarity and compatibility are the most useful features when selecting
formulas, each having a significant increase in F1 over the proximity
baseline. There are minor gains from combining these two features. On the
other hand, semantic similarity does not affect performance relative to the
baseline. We find that this is mainly due to the disproportionate number of
unfamiliar formulas present in the dataset that drown out any signal. Table 4
presents two examples of the system’s ranking of formulas.

Feature set | Train | Dev
---|---|---
| P | R | F1 | P | R | F1
Proximity | 56.4 | 48.7 | 52.2 | 56.3 | 48.8 | 52.3
Similarity | 65.1 | 34.9 | 45.4 | 65.1 | 34.9 | 45.4
Familiarity∗ | 70.5 | 63.5 | 66.8 | 69.6 | 62.9 | 66.1
Compatibility+ | 66.9 | 74.4 | 70.4 | 65.4 | 73.1 | 69.0
F + C† | 73.8 | 70.3 | 72.1 | 71.5 | 68.9 | 70.1
F + C + P† | 73.8 | 70.3 | 72.1 | 71.5 | 68.9 | 70.1
F + C + P + S† | 73.8 | 70.3 | 72.0 | 71.4 | 68.6 | 69.9
(a)  the formula construction system. Precision, Recall and F1 are cross-
validated on 10-folds. ∗significant F1 versus P and S with $p<0.01$.
+significant F1 versus P, S and F with $p<0.01$. †significant F1 versus P, S,
F and C with $p<0.05$.

System |  Train BLEU |  Test BLEU
---|---|---
Baseline |  65.00 |  57.32
RNN∗ |  81.50 |  69.79
(b)  the description generation system. ∗significant BLEU score versus the
baseline with $p<0.01$.

Table 5:  Evaluation of perspective generation subsystems.

##  5 Perspective generation

Our next goal is to generate natural language descriptions, also known as
perspectives, given a formula. Our approach models the task as a sequence-to-
sequence translation task from formulas to natural language. We first describe
a rule-based baseline and then describe a recurrent neural network (RNN) with
an attention-based copying mechanism (Jia and Liang, 2016).

### Baseline.

As a simple approach to generate perspectives, we just combine tuples in the
formula with the neutral prepositions of and for, e.g. “1/5th of the cost of
an employee for the population of Texas for the time taken for lunch.”

### Sequence-to-sequence RNN.

We use formula-perspective pairs from the dataset to create a sequence-to-
sequence task: the input is composed using the formula’s multiplier and
descriptions of its tuples connected with the symbol ‘*’; the output is the
perspective (Figure 7).

Figure 7:  We model description generation as a sequence transduction task,
with input as formulas (at bottom) and output as perspectives (at top). We use
a RNN with an attention-based copying mechanism.

Our system is based on the model described in Jia and Liang (2016). Given a
sequence of input tokens ($\mathbf{x}=(x_{i})$), the model computes a context-
dependent vector ($\mathbf{b}=(b_{i})$) for each token using a bidirectional
RNN with LSTM units. We then generate the output sequence ($y_{j}$) left to
right as follows. At each output position, we have a hidden state vector
($s_{j}$) which is used to produce an “attention” distribution
($\mbox{\boldmath$\alpha$}_{j}=(\alpha_{ji})$) over input tokens:
$\alpha_{ji}=\operatorname{Attend}(s_{j},b_{i})$. This distribution is used to
generate the output token and update the hidden state vector. To generate the
token, we either sample a word from the current state or copy a word from the
input using attention. Allowing our model to copy from the input is helpful
for our task, since many of the entities are repeated verbatim in both input
and output. We refer the reader to Jia and Liang (2016) for more details.

### Evaluation.

Input formula |  Generated perspective
---|---
7 $\times$ the cost of an employee $\times$ a week |  7 times the cost of employing one person for one week
1/10 $\times$ the cost of an employee $\times$ the population of California $\times$ the time taken for a football game |  one tenth the cost of an employee during a football game by the population of California
1 $\times$ coffee consumption $\times$ a minute $\times$ population of the world |  the amount of coffee consumed in one minute on the world
6 $\times$ weight of a person $\times$ population of California |  six times the weight of the people who is worth
Table 6:  Examples of perspectives generated by the sequence-to-sequence RNN.
The model is able to capture rephrasings of fact descriptions and reordering
of the facts. However, it often confuses prepositions and, very rarely, can
produce nonsensical utterances.

We split the perspective description dataset into a training and test set such
that no formula in the test set contains the same set of tuples as a formula
in the training set.777Note that formulas with the same set of tuples can
occur multiple times in the either the training or test set with different
multipliers. Table 7(b) compares the performance of the baseline and sequence-
to-sequence RNN using BLEU. The sequence-to-sequence RNN performs
significantly better than the baseline, producing more natural rephrasings.
Table 6 shows some output generated by the system (see Table 6).

##  6 Human evaluation

lr+rnn | baseline |
---|---|---
perspective rated useful? | #
Yes | Yes | 31
Yes | No | 63
No | Yes | 61
No | No | 56

(a)  A summary of the number of times the perspective generated by lr+rnn or
baseline was rated useful by a majority of crowdworkers.

Cause of error | #
---|---
Proximity | 9
Familiarity | 6
Compatibility | 8
Similarity | 49
Generation | 24
Other | 14
Invalid mention | 7
Total | 117

(b)  An analysis of errors produced by lr+rnn when its perspectives were not
rated useful. Errors caused by poor formula selection are further categorized
by selection criteria violated.

Cat. |  Mention
---|---
|  lr+rnn perspective (vs. baseline)
Prox. |  …ready to ship about 2,300 miles across the Pacific to the mainland …
|  three times the distance from San Francisco to Los Angeles (vs. the
distance from San Francisco to Dallas TX).
Sim. |  China had disposed of about 100,000 tons of CFCs” …
|  one fifth of the weight of garbage produced in the United States by the
population of Texas in one week. (vs. the average food wasted every year).
Fam. |  …the project could save New England ratepayers $4.6 billion in energy costs over 25 years.
|  one eighth the cost of employing the population of Asia for one hour. (vs.
the construction cost of The Cosmopolitan in Las Vegas.)
Comp. |  Hominids started shaping stone tools about 2.6 million years ago.
|  5 times the total time taken to build the number of cars registered. (vs.
17000 times the average lifetime for a tree).
(c)  Examples of errors categorized by the criteria defined in Section 4.

Table 7:  Results of an end-to-end human evaluation of the output produced by our perspective generation system (lr+rnn) and a baseline (baseline) that picks the numerically closest tuple in the knowledge base for each mention. |  Mention |  Perspective (that’s about…)
---|---|---
+ |  In 2007, Turkmenistan exported 50 billion cubic meters of gas to Russia. |  the amount of oil produced by the US during a lifetime
+ |  It can carry up to 10 nuclear warheads and has a range of 8,000 km. |  the distance from San Francisco to Beijing
- |  the 2.7 million square feet that Mission Bay’s largest developer is entitled to build |  twice the area of forest logged in a minute
- |  Las Vegas Sands claims the 10.5 million square feet is the largest building in Asia. |  one half of an area of an average farm
Table 8:  Examples of perspectives generated by our system that frame the
mentioned quantity to be larger or smaller (top to bottom) than initially the
authors thought.

In addition to the automatic evaluations for each component of the system, we
also ran an end-to-end human evaluation on an independent set of 211 mentions
collected using the same methodology described in Section 3. Crowdworkers were
asked to choose between perspectives generated by our full system (lr+rnn) and
those generated by the baseline of picking the numerically closest tuple in
the knowledge base (baseline). They could also indicate if either both or none
of the shown perspectives appeared useful.888Crowdworkers were paid $0.06 per
to choose a perspective for each mention. Each mention and set of perspectives
were presented to 5 crowdworkers.

Table 7 summarizes the results of the evaluation and an error analysis
conducted by the authors. Errors were characterized as either being errors in
generation (e.g. Table 6) or violations of the criteria in selecting good
formulas described in Section 4 (Table 7). The other category mostly contains
cases where the output generated by lr+rnn appears reasonable by the above
criteria but was not chosen by a majority of workers. A few of the mentions
shown did not properly describe a numeric quantity, e.g. “…claimed
responsibility for a 2009 gun massacre …” and were labeled invalid mentions.
The most common error is the selection of a formula that is not contextually
relevant to the mentioned text because no such formula exists within the
knowledge base (within an order of magnitude of the mentioned value): a larger
knowledge base would significantly decrease these errors.

##  7 Related work and discussion

We have proposed a new task of perspective generation. Compositionality is the
key ingredient of our approach, which allows us synthesize information across
multiple sources of information. At the same time, compositionality also poses
problems for both formula selection and description generation.

On the formula selection side, we must compose facts that make sense. For
semantic compatibility between the mention and description, we have relied on
simple word vectors (Mikolov et al., 2013), but more sophisticated forms of
semantic relations on larger units of text might yield better results (Bowman
et al., 2015).

On the description generation side, there is a long line of work in generating
natural language descriptions of structured data or logical forms Wong and
Mooney (2007); Chen and Mooney (2008); Lu and Ng (2012); Angeli et al. (2010).
We lean on the recent developments of neural sequence-to-sequence models
(Sutskever et al., 2014; Bahdanau et al., 2014; Luong et al., 2015). Our
problem bears some similarity to the semantic parsing work of Wang et al.
(2015), who connect generated canonical utterances (representing logical
forms) to real utterances.

If we return to our initial goal of helping people understand numbers, there
are two important directions to explore. First, we have used a small knowledge
base, which limits the coverage of perspectives we can generate. Using
Freebase (Bollacker et al., 2008) or even open information extraction (Fader
et al., 2011) would dramatically increase the number of facts and therefore
the scope of possible perspectives.

Second, while we have focused mostly on basic compatibility, it would be
interesting to explore more deeply how the juxtaposition of facts affects
framing. Table 8 presents several examples generated by our system that frame
the mentioned quantities to be larger or smaller than the authors originally
thought. We think perspective generation is an exciting setting to study
aspects of numeric framing (Teigen, 2015).

### Reproducibility

All code, data, and experiments for this paper are available on the CodaLab
platform at
https://worksheets.codalab.org/worksheets/0x243284b4d81d4590b46030cdd3b72633/.

## Acknowledgments

We would like to thank Glen Chiacchieri for providing us information about the
Dictionary of Numbers, Maneesh Agarwala for useful discussions and references,
Robin Jia for sharing code for the sequence-to-sequence RNN, and the anonymous
reviewers for their constructive feedback. This work was partially supported
by the Sloan Research fellowship to the second author.

## References

  * Angeli et al. (2010) G. Angeli, P. Liang, and D. Klein. 2010.  A simple domain-independent probabilistic approach to generation.  In Empirical Methods in Natural Language Processing (EMNLP).
  * Bahdanau et al. (2014) D. Bahdanau, K. Cho, and Y. Bengio. 2014.  Neural machine translation by jointly learning to align and translate.  arXiv preprint arXiv:1409.0473 .
  * Barrio et al. (2016) P. J. Barrio, D. G. Goldstein, and J. M. Hofman. 2016.  Improving the comprehension of numbers in the news.  In Conference on Human Factors in Computing Systems (CHI).
  * Berg-Kirkpatrick et al. (2012) T. Berg-Kirkpatrick, D. Burkett, and D. Klein. 2012.  An empirical investigation of statistical significance in NLP.  In Empirical Methods in Natural Language Processing (EMNLP). pages 995–1005.
  * Bollacker et al. (2008) K. Bollacker, C. Evans, P. Paritosh, T. Sturge, and J. Taylor. 2008.  Freebase: a collaboratively created graph database for structuring human knowledge.  In International Conference on Management of Data (SIGMOD). pages 1247–1250.
  * Bowman et al. (2015) S. Bowman, G. Angeli, C. Potts, and C. D. Manning. 2015.  A large annotated corpus for learning natural language inference.  In Empirical Methods in Natural Language Processing (EMNLP).
  * Chen and Mooney (2008) D. L. Chen and R. J. Mooney. 2008.  Learning to sportscast: A test of grounded language acquisition.  In International Conference on Machine Learning (ICML). pages 128–135.
  * Chevalier et al. (2013) F. Chevalier, R. Vuillemot, and G. Gali. 2013.  Using concrete scales: A practical framework for effective visual depiction of complex measures.  IEEE Transactions on Visualization and Computer Graphics 19:2426–2435.
  * Chiacchieri (2013) G. Chiacchieri. 2013.  Dictionary of numbers.  http://www.dictionaryofnumbers.com/.
  * Fader et al. (2011) A. Fader, S. Soderland, and O. Etzioni. 2011.  Identifying relations for open information extraction.  In Empirical Methods in Natural Language Processing (EMNLP).
  * Jia and Liang (2016) R. Jia and P. Liang. 2016.  Data recombination for neural semantic parsing.  In Association for Computational Linguistics (ACL).
  * Jones and Taylor (2009) M. G. Jones and A. R. Taylor. 2009.  Developing a sense of scale: Looking backward.  Journal of Research in Science Teaching 46:460–475.
  * Kim et al. (2016) Y. Kim, J. Hullman, and M. Agarwala. 2016.  Generating personalized spatial analogies for distances and areas.  In Conference on Human Factors in Computing Systems (CHI).
  * Lu and Ng (2012) W. Lu and H. T. Ng. 2012.  A probabilistic forest-to-string model for language generation from typed lambda calculus expressions.  In Empirical Methods in Natural Language Processing (EMNLP). pages 1611–1622.
  * Luong et al. (2015) M. Luong, H. Pham, and C. D. Manning. 2015.  Effective approaches to attention-based neural machine translation.  In Empirical Methods in Natural Language Processing (EMNLP). pages 1412–1421.
  * Mikolov et al. (2013) T. Mikolov, K. Chen, G. Corrado, and Jeffrey. 2013.  Efficient estimation of word representations in vector space.  arXiv .
  * Paulos (1988) J. A. Paulos. 1988.  Innumeracy: Mathematical illiteracy and its consequences.  Macmillan.
  * Seife (2010) C. Seife. 2010.  Proofiness: How you’re being fooled by the numbers.  Penguin.
  * Sutskever et al. (2014) I. Sutskever, O. Vinyals, and Q. V. Le. 2014.  Sequence to sequence learning with neural networks.  In Advances in Neural Information Processing Systems (NIPS). pages 3104–3112.
  * Teigen (2015) K. H. Teigen. 2015.  Framing of numeric quantities.  The Wiley Blackwell Handbook of Judgment and Decision Making pages 568–589.
  * Tretter et al. (2006) T. R. Tretter, M. G. Jones, and J. Minogue. 2006.  Accuracy of scale conceptions in science: Mental maneuverings across many orders of spatial magnitude.  Journal of Research in Science Teaching 43:1061–1085.
  * Wang et al. (2015) Y. Wang, J. Berant, and P. Liang. 2015.  Building a semantic parser overnight.  In Association for Computational Linguistics (ACL).
  * Wong and Mooney (2007) Y. W. Wong and R. J. Mooney. 2007.  Generation by inverting a semantic parser that uses statistical machine translation.  In Human Language Technology and North American Association for Computational Linguistics (HLT/NAACL). pages 172–179.


11institutetext: Department of Astrophysics/IMAPP, Radboud University, P.O.
Box 9010, 6500 GL Nijmegen, The Netherlands 22institutetext: SRON, Netherlands
Institute for Space Research, Niels Bohrweg 4, 2333 CA Leiden, The Netherlands
33institutetext: Institute of Astronomy, KU Leuven, Celestijnenlaan 200D,
B-3001 Leuven, Belgium

# Failed supernovae as a natural explanation for the binary black hole mass
distribution

P. Disberg<EMAIL_ADDRESS>G. Nelemans 112233

###### Abstract

Context. The more gravitational wave sources are detected, the better the mass
distribution of binary black holes (BBHs) becomes known. This “stellar
graveyard” shows several features, including an apparent mass gap which makes
the distribution bimodal. The observed chirp mass distribution, in turn,
appears to be trimodal.

Aims. We aim to investigate to which extend we can explain the observed mass
distribution with stellar evolution, specifically with the hypothesis that the
mass gap is caused by the difference between successful and failed supernovae
(SNe).

Methods. We pose a hypothetical remnant function, based on literature of
stellar evolution simulations, which relates initial mass to remnant mass,
includes a “black hole island” and produces a bimodal remnant distribution.
Moreover, we look at observed type II SN rates in an attempt to detect the
effect of failed SNe. Finally, using a simplified estimation of binary
evolution, we determine the remnant distribution resulting from our remnant
function and compare it to observation.

Results. We find that failed SNe lower type II SN rates by approximately
$25\%$, but the inferred rate from SN surveys is not accurate enough to
confirm this. Furthermore, our estimation based on the remnant function
produces a mass distribution that matches the general shape of the observed
distributions of individual as well as chirp masses.

Conclusions. Based on our research, we conclude that the failed SNe mechanism
and the presence of the black hole island are a natural hypothesis for
explaining the individual BBH mass distribution and chirp mass distribution.
However, for a more firm conclusion more detailed simulations are needed.

###### Key Words.:

gravitational waves – gravitational lensing: strong – supernovae: general

##  1 Introduction

Figure 1: GW data from Abbott et al. (2021a). The top panel shows the mass
distribution of the individual BBHs (represented by a histogram of the mean
values, with a bin-width of $2M_{\odot}$). We can also represent the
datapoints as asymmetric Gaussians, as described in appendix A. The sum of
these Gaussians, $\sum_{i}dN_{i}/dm$, is given as well (solid line), together
with a cumulative distribution function (CDF, dotted line on a different
axis). The bottom panel shows the chirp mass distribution (grey lines) and is
made to resemble the distribution from Abbott et al. (2021b). In order to do
this, we use asymmetric Gaussians as well, which take into account the $90\%$
credible interval. We also show an adjustable-width kernel density estimation
(AWKDE, dashed line), from Shimazaki & Shinomoto (2010); Shimazaki et al.
(2018). The Gaussians are plotted on a logarithmic axis, while the sum and the
AWKDE use a linear axis. Similarly to Abbott et al. (2021b), we omit GW190814
from our analysis.

The more gravitational wave (GW) sources are detected, the better the
distribution of the masses of the detected remnants becomes known (e.g. Abbott
et al., 2021c). These remnants are often black holes (BHs) in a binary system,
which merge with each other after their binary evolution, through the emission
of GWs (e.g. Abbott et al., 2016). This evolution includes the collapse into
BHs, potentially accompanied by supernovae (SNe) of both stars, which means
the SN mechanism can be important in determining the final remnant mass. We
are interested in the distribution of these final masses.
The top panel of fig. 1 shows the distribution of the individual masses of the
detected GW sources, from GWTC 1, 2 and 3 (Abbott et al., 2021a). The binary
black holes (BBHs) seem to show a bimodal distribution, with a gap between
$14M_{\odot}$ and $22M_{\odot}$, and an additional peak at $35M_{\odot}$,
although this peak seems to disappear when taking into account the uncertainty
of the data. The bottom panel shows the distribution of the corresponding
chirp masses of the black hole binaries (BHBs), represented by asymmetric
Gaussians (as described in appendix A). This distribution appears to be
trimodal (Abbott et al., 2021b), with one mode below $10M_{\odot}$, one
between $10M_{\odot}$ and $20M_{\odot}$, and one above $20M_{\odot}$. Here,
there is a gap as well: at $11M_{\odot}$.
One hypothethical explanation for this is given by Broadhurst et al. (2018,
2022), who theorize that the gap, specifically the gap in the individual mass
distribution, is real and caused by gravitational lensing, because of which
the distance to some of the GW sources is underestimated and therefore the
chirp mass overestimated. The gap is then located between the non-lensed and
lensed sources. In order to produce a distribution similar to the top panel of
fig. 1, Broadhurst et al. (2018, 2022) pose a merger rate which has high
values at high redshift and low values at low redshift, since this produces
the desired number of lensed and unlensed sources. We estimate, however, that
the merger rate value they pose at high redshift implies a merger fraction
which is approximately $50$ times higher than the value Mandel & Farmer (2022)
estimate through their Drake-like approach. Even if we make optimistic
assumptions about some of the factors in this Drake-like equation, we still
find a factor of $14$ unaccounted for (as shown in appendix B). Although this
factor is not large enough to completely dismiss the lensing hypothesis of
Broadhurst et al. (2018, 2022), we find it sufficient to deem their hypothesis
improbable.
We aim to explore whether the supernova (SN) mechanism could provide an
alternative, more natural explanation for the observed gap. Literature
suggests that the shock-wave in a core-collapse SN (CCSN) can be stalled and
not cause a successful explosion (Mazurek, 1982). The star can then collapse
in its totality and form a stellar remnant in a direct collapse, or “failed”,
SN. Failed SNe are difficult to detect, since they do not cause a bright
explosion but instead make a massive star seemingly disappear (Kochanek et
al., 2008). Because of these difficulties, the properties of failed SNe are
not well-known, even though there are several potential candidates (e.g. Gerke
et al., 2015; Adams et al., 2017a, b; Basinger et al., 2021; Neustadt et al.,
2021). In addition to the search for failed SNe, others have also looked at
the implications for remnant mass distributions, e.g. Kochanek (2014a, b) who
attempts to estimate the low-mass BH distribution based on the failed SN
mechanism. Overall, we state that the difference between failed SNe and
successful ones could potentially create a gap in the mass distribution: some
sources have a successful SN and lose mass before forming a remnant and other
sources have a failed SN, conserve mass and form a more massive remnant. This
means that the parameter which is important for the remnant mass distribution
is the mass limit above which SNe fail.
In this work we investigate the SN mechanism in close binaries and how it
affects the BBH mass distribution (when we mention SNe we refer to CCSNe,
unless specified otherwise). We start by investigating this SN mechanism in
section 2, wherein we describe the transitions between successful and failed
SNe through an initial-final mass relation based on the works of Schneider et
al. (2021) and Marchant et al. (2019) (section 2.1). We also attempt to
compare this relation to observed SN rates (section 2.2). Then, in section 3,
we show that this model naturally produces a bimodal mass distribution and a
trimodal chirp mass distribution, similar to fig. 1. While our paper was under
review, Schneider et al. (2023) published a paper with a similar approach. In
the discussion (section 4) we will make a comparison. Finally, in section 5,
we will summarize our conclusions.

##  2 Failed supernovae

###  2.1 Remnant function

Figure 2: Remnant masses ($m_{\text{RN}}$), as function of ZAMS mass
($m_{\text{ZAMS}}$). The blue curve shows the remnant function for case A and
case B mass transfer, based on the simulations by Schneider et al. (2021) and
Marchant et al. (2019). The Marchant et al. (2019) results are shown as well
(black dots). The triangle shows the transition point between the part of our
function which is based on Schneider et al. (2021) and the part which is based
on Marchant et al. (2019). The red curve shows the corresponding remnant
function for case C mass transfer, where the part after the transition point
is the PPISN fit translated to lower masses by $42M_{\odot}$. The dashed lines
show $m_{\text{RN}}/m_{\text{ZAMS}}=1,\ 0.65$ and $0.35$, since the latter two
are used as approximation in estimating the amount of mass transfer in our
estimation. The shaded area shows the values of the remnant mass gap.
Furthermore, we note that this remnant function is limited to BHs, for our
purposes we are not interested in the neutron star (NS) masses. The remnant
functions are explicitly given in appendix C.

We are interested in the differences between successful and failed SNe, and
the relation between zero-age main-sequence (ZAMS) mass and the masses of the
remnants they form. SN explodability simulations indicate that there may not
be only one mass range for successful SNe and one range for failed SNe (Ertl
et al., 2016; Müller et al., 2016; Kresse et al., 2021). Instead, they find
that there may be a failed SN range somewhere between $20M_{\odot}$ and
$24M_{\odot}$, above which stars can have a successful SN again up until about
$27M_{\odot}$, where the failed SNe take over again, supported by other
studies of pre-SN compactness (e.g. Sukhbold & Woosley, 2014; Sukhbold et al.,
2016). We combine this with the results of Schneider et al. (2021), who
perform a simulation and investigate stellar evolution and remnant masses.
They also find that there are two ranges for failed SNe: one small range for
lower masses, which they call the “BH island”, and a larger failed SN range
for higher masses. They note that the two points of transition from successful
to failed SNe coincide with the masses at which the stellar core change from
convective to radiative carbon and neon burning, respectively. This could
explain why there is a BH island in the first place. Schneider et al. (2021)
also note that their relation between ZAMS mass and remnant mass gives rise to
a bimodal mass distribution.
We construct a possible remnant function, starting by describing the BH
island, based on Schneider et al. (2021), both for their models which concern
case A/B mass transfer as well as their case C models. They describe the BH
island as a function which produces remnants between approximately
$8M_{\odot}$ and $10M_{\odot}$. In order to reproduce the observed BBH
distribution, however, we change this to a linear function which produces
remnants between $8M_{\odot}$ and $14M_{\odot}$. We give the exact definition
of our remnant function, including these mass ranges, in appendix C. For the
other failed SNe, we use the fits given by Schneider et al. (2021), but
because these fits start at a remnant mass of about $16M_{\odot}$, we shift
them upwards to $22M_{\odot}$ instead. We motivate this by general uncertainty
in the models plus the fact that rather than complete collapse or complete
ejection, there could be partial fallback, so that some stars which are in the
successful SN domain will in fact be able to form BHs. Finally, at the highest
masses, we include the effect of pair-instability SNe (PISN) that leave behind
no stellar remnant (Fraley, 1968) and pulsational pair-instability supernova
(PPISN) (Woosley et al., 2007; Woosley, 2017), at slightly lower masses that
can form the transition between the remnants of failed SNe and the remnantless
PISNe. For this, we make a polynomial fit to the PPISN results from Marchant
et al. (2019). In order to create one coherent remnant function, we increase
the slope of the case A/B function, connecting it to this PPISN fit, and also
shift the PPISN fit for case C to lower ZAMS masses, connecting it to the case
C function. The first adjustment does not influence our results significantly
and the second adjustment is justified because we expect case C stars to have
a more massive helium core at the end of their evolution than case A/B stars
have. This is due to the fact that case C stars are able to grow more massive
cores, meaning we expect case C stars to have a lower threshold for PPISNe. We
show our complete remnant functions, one for case A/B and one for case C, in
fig. 2. Here, we neglect influences of aspects such as metallicity, which we
discuss in section 4. Nevertheless, this remnant function suffices for our
purpose: showing that the bimodal BBH distribution is to be expected based on
stellar evolution (which we do in section 3).

###  2.2 Supernova rates

Figure 3: Data from SN surveys (Graur et al., 2015, 2017; Botticella et al.,
2017). The top panel shows the type II SN rates per unit mass
($R_{\text{SNuM}}$) versus the SFR per unit mass (sSFR), together with a fit
for all the data-points (dashed line) and a fit which excludes the two
outliers at $\text{sSFR}\approx 1.5\cdot 10^{-12}\text{yr}^{-1}\text{ and
}1.5\cdot 10^{-9}\text{yr}^{-1}$ (dotted line). The fits use the least-squares
method on a first degree polynomial in logarithmic space and the error bars
are one standard deviation. These data are corrected for the Chabrier (2003)
IMF, since the original data follows the Salpeter (1955) IMF. For a constant
SFR, $L=\kappa\psi\int\tau\phi L(m)dm$, where $L(m)\propto m^{3.5}$ is the
luminosity of a star with mass $m$. However, since $\tau(m)\cdot L(m)\propto
m$, it cancels out with $\kappa$ and makes the SFR independent of the IMF.
However, the SN rate is proportional to $\kappa$, so it is multiplied by a
factor $\kappa_{\text{Chabrier}}/\kappa_{\text{Salpeter}}\approx 1.38$,
meaning that this correction also multiplies $R_{\text{SNuM}}/\text{sSFR}$ by
a factor of $1.38$. The bottom panel shows a histogram of the
$R_{\text{SNuM}}/\text{sSFR}$ ratio (with a bin-width of $2\cdot
10^{-3}M_{\odot}^{-1}$), and a kernel density estimation (KDE).

Before we turn to the mass distribution resulting from the remnant function,
we try to detect failed SNe by investigating how they reduce SN rates. After
all, since failed SNe do not produce the bright signal successful SNe do, an
increase in amount of failed SNe would reduce the total SN rate. In order to
make a rough estimate of the effect of failed SNe on SN rates, we assume one
mass range of stars which have a successful SN, from $m_{l}$ to $m_{u}$. Above
$m_{u}$, stars are too heavy to explode and have a failed SN, meaning they do
not contribute to the SN rate. This model differs from our remnant function,
since here we assume only one mass range of failed SNe, but for our purposes
this suffices. In this estimation, we neglect the influences of binary
interaction, effectively limiting our scope to type II SNe.
We start by describing the SN rates. Using the initial mass function (IMF),
$N(m)dm=\kappa M\phi(m)dm$ where $\kappa^{-1}=\int m\phi(m)dm$ and $\phi(m)$
the Chabrier (2003) IMF shape, we express the total mass $M$ in terms of the
star formation rate (SFR) $\psi(t)$. This comes down to the total mass of the
stars of mass $m$ which are created a lifetime $\tau$ ago and equals
$M=\psi(t-\tau(m))dt$. Therefore, the SN rate is simply the amount of stars
within a certain mass range which reach the end of their stellar lifetime at
time $t$ and is defined as:

$R_{\text{SN}}(t,m_{l},m_{u})=\kappa\int\dfrac{M}{dt}\phi(m)dm=\kappa\int_{m_{l}}^{m_{u}}\psi(t-\tau)\phi(m)dm\quad,$
(1)

where we use $\tau=\tau(m)=10^{10}\text{yr}\left(m/M_{\odot}\right)^{-2.5}$.
Since the precise shape of the SFR is often not determined accurately, it is
difficult to determine $R_{\text{SN}}$ using this equation. Because of this,
we look at the ratio between the SN rate and the SFR. The integrand of this
ratio has a factor $\psi(t-\tau)/\psi(t)$, which is why we make the
approximation $\psi(t-\tau)/\psi(t)\approx 1$. This approximation holds for
constant SFRs and decreases in accuracy for SFRs which evolve on a short
timescale. Also, we set $m_{l}=8M_{\odot}$, which is non-controversial (e.g.
Kochanek et al., 2008; Kochanek, 2014b). The ratio of SN rate to SFR becomes,
then:

$\dfrac{R_{\text{SN}}(t,8M_{\odot},m_{u})}{\psi(t)}\approx\dfrac{\int_{8M_{\odot}}^{m_{u}}\phi(m)dm}{\int
m\phi(m)dm}\equiv\Lambda(m_{u})\quad,$ (2)

where we used the definition of $\kappa$ and introduce $\Lambda(m_{u})$ as
shorthand for this function. The $\kappa$ integral is taken from
$m_{min}=10^{-1}M_{\odot}$ up to $m_{\max}=10^{2}M_{\odot}$, which means the
value which does not take failed SNe into account is simply
$\Lambda(m_{\max})$.
In order to determine $m_{u}$, we consider the SN simulations of Ertl et al.
(2016), Müller et al. (2016) and Kresse et al. (2021), together with our
remnant function based on Schneider et al. (2021). These simulations use
different models of neutrino engines, which represent the collapsed core, and
determine the mass dependent explodability per model. Müller et al. (2016)
find a mass dependent probability distribution, where approximately
$m<20.5M_{\odot}$ and $23.5M_{\odot}<m<27M_{\odot}$ have a high probability
for a successful SN. The simulations of Ertl et al. (2016) and Kresse et al.
(2021) find a similar distribution, with approximately $m<21.5M_{\odot}$ and
$25M_{\odot}<m<27.5M_{\odot}$ which go SN. These are in good agreement with
Schneider et al. (2021), who find that their case C results (as shown in fig.
2) are similar to the results for single stars. Keeping in mind the IMF, we
make a crude approximation and estimate, based on these simulations, that
$m_{u}\approx 22.5M_{\odot}$. This is comparable to the estimation of
$20M_{\odot}$ which Mashian & Loeb (2017) use, for instance.
The value of $m_{u}$ gives $\Lambda(22.5M_{\odot})=9.06\cdot
10^{-3}M_{\odot}^{-1}$. In contrast, if we neglect the influence of failed
SNe, this value becomes $\Lambda(m_{\max})=1.18\cdot 10^{-2}M_{\odot}^{-1}$.
This means we have found a failed SNe correction of
$1-\Lambda(22.5)/\Lambda(m_{\max})=0.232\approx 25\%$.
We now turn to SN surveys and use the works of Graur et al. (2015, 2017) and
Botticella et al. (2017) in order to detect the predicted effect of failed
SNe. As shown in the top panel of fig. 3: these data consist of the type II SN
rate per unit mass ($R_{\text{SNuM}}$) versus the specific SFR (sSFR), i.e.
the SFR per unit mass. The figure includes two fits, shaped as
$R_{\text{SNuM}}=a\cdot\text{sSFR}^{b}$, where the first one includes all the
data-points and the second one neglects two outliers. The fitted values of
$b$, $0.68$ and $1.04$, confirm that our assumption
$\psi(t-\tau)/\psi(t)\approx 1$ (which means $R_{\text{SN}}\propto\psi$ and
$b=1$) is appropriate. If we we assume $b\approx 1$, the two points which are
neglected for the second fit indeed appear to be outliers. The consequence of
this assumption is that $a=R_{\text{SNuM}}/\text{sSFR}=\Lambda$. Our estimate
$\Lambda(22.5M_{\odot})=9.06\cdot 10^{-3}M_{\odot}^{-1}$ differs significantly
from the fitted values of $a$: $4.0\cdot 10^{-6}M_{\odot}^{-1}$ and $1.9\cdot
10^{-2}M_{\odot}^{-1}$. If we set $b=1$, the fitted values of $a$ become
$\left(7.2\pm 1.7\right)\cdot 10^{-3}M_{\odot}^{-1}$ and $\left(7.3\pm
0.9\right)\cdot 10^{-3}M_{\odot}^{-1}$ respectively, which is consistent with
$\Lambda$. The bottom panel of fig. 3 shows the distribution of the individual
values of $R_{\text{SNuM}}/\text{sSFR}$, including a kernel density
estimation. The distribution shows an outlier at $3\cdot
10^{-2}M_{\odot}^{-1}$, which corresponds to the outlier at $1.5\cdot
10^{-12}\text{yr}^{-1}$ in the top panel which was neglected in the second
fit. Also, a peak is situated at approximately $0.8\cdot
10^{-2}M_{\odot}^{-1}$, which is approximately the value of the fitted $a$
with $b=1$.
The question arises whether uncertainty in other aspects of our model can
account for a similar deviation from the standard $\Lambda(m_{\max})$. This
deviation can, for instance, be accomplished by choosing $m_{l}$ to be
$9.72M_{\odot}$, which seems high compared to the literature value of
$8M_{\odot}$, although not unreasonably. Furthermore, uncertainties in the IMF
could also account for such a deviation. Not only does Chabrier (2003) give
different parameter values for binary systems (which would cause a deviating
value of $\Lambda$), but there is also uncertainty in the high-mass slope of
the IMF. Parravano et al. (2018) state that this slope has a value of
$\alpha=2.35_{-0.15}^{+0.35}$. A value of $2.45$ is enough to account for the
$25\%$ deviation, which is well within this range. This means that even if
$\Lambda(22.5M_{\odot})$ fits the data better than $\Lambda(m_{\max})$, which
the data may suggest, it would be difficult to attribute this to failed SNe.
Because of this, we cannot conclude that we find the effects of failed SNe in
the SN survey data and therefore cannot constrain the value of $m_{u}$ based
on this analysis. Despite not being able to confirm the failed SNe
observationally, we are interested in how the remnant function from section
2.1 shapes the BBH mass distribution.

##  3 Binary black hole mass distribution

Figure 4: BBH mass distribution resulting from our estimation, where we
simulate approximately $N=1.2\cdot 10^{6}$ binaries and set $\eta=0.5$ and
$\zeta=0.5$. The latter means that half of the binaries follow the case A/B
remnant function, and the other half follow the case C remnant curve (as shown
in fig. 2). The top row shows the resulting mass distribution of the
individual primary (red) and secondary (grey) masses, in a stacked histogram
with a bin-width of $1M_{\odot}$. The bottom row shows the chirp mass
distribution in a similar histogram, where we distinguish between four
possible configurations of the primary and secondary mass with respect to the
central value of the mass gap ($m_{G}$): either both BHs are below the gap
(dark blue), the secondary is below the gap while the primary is above it
(blue) or vice versa (beige), or both are above the gap (light blue). The left
column shows the exact results of our estimation, while the right column adds
a Gaussian uncertainty to the results in order to make the simulated
distribution look more similar to the observed one (as shown in fig. 1). We
base the Gaussian uncertainty on the average standard deviations of the GW
data, as described in appendix D.

We now want to show that our remnant function indeed produces a BBH
distribution similar to observation. In order to do this, we make a simple
estimation of binary evolution and determine the resulting remnant mass
distribution. In our estimation, we start with approximately $1.2\cdot 10^{6}$
binaries, in which the primary masses are distributed according to the
Chabrier (2003) IMF and the secondary masses are distributed uniformly. This
means that we define the primary star here to have the highest initial mass.
In our model, we also include mass transfer from the primary to the secondary
in a Roche-lobe overflow (RLO) phase, and neglect mass transfer from the
secondary to the primary. We state that the secondary accretes a fraction
$\eta$ of the mass expelled from the primary during the RLO phase. The precise
value of $\eta$ is uncertain (Dorozsmai & Toonen, 2022), so we will simply set
$\eta=0.5$, since it does not affect the general shape of the remnant
distribution. In order to determine the amount of mass which is expelled from
the primary, we cannot simply take the difference between ZAMS mass and
remnant mass, since in successful SNe there is mass loss outside of the RLO
phase. We therefore approximate the expelled mass of the primary in the RLO
phase by $0.35\cdot m_{\text{ZAMS}}$ for case A/B and $0.65\cdot
m_{\text{ZAMS}}$ for case C, since these curves approximate the failed SN
curves, as shown in fig. 2. The secondary, then, accretes a fraction $\eta$ of
this expelled mass and is rejuvenated, which means that we use this new mass
as if it were the ZAMS mass. Besides $\eta$ we also define the parameter
$\zeta$, which equals the fraction of binaries which follow the case C curve.
Since $\zeta$ also does not influence the general shape of the distribution
significantly, we set $\zeta=0.5$ as well.
The estimated mass distribution is shown in fig. 4. The left panels show the
exact results from our estimation and the right panels show the results with
an added Gaussian uncertainty, as described in appendix D. Unsurprisingly, the
remnant mass distribution is indeed bimodal: it shows one peak at lower
masses, caused by the BH island, and one peak at higher masses. The exact
results of our estimation, without the Gaussian uncertainty, also show a PPISN
peak around $44M_{\odot}$, caused by the flat PPISN distribution in fig. 2.
The general shape of the estimated distribution is similar to the observed
distribution, although it is difficult to compare the two in detail because
the relative heights of the peaks are influenced by multiple aspects. Not only
is there a bias towards high-mass mergers in observation, but the parameters
$\eta$ and $\zeta$ also influence the heights. The effects of these
parameters, including the fact that varying them does not influence our
conclusions about the general distribution, is shown in appendix E.
The resulting chirp mass distribution also shows interesting features. Since
the individual mass distribution is bimodal, with one peak at either side of
the mass gap at $m_{G}=17M_{\odot}$, the chirp mass distribution shows four
possible configurations, with the primary and secondary mass at the same and
opposite sides of the mass gap. These configurations can clearly be found in
the chirp mass distribution: there is one peak centered just below
$10M_{\odot}$ for binaries with both BHs below the mass gap, one peak centered
around $28M_{\odot}$ for binaries with both BHs above the mass gap, and also
two peaks centered around $15M_{\odot}$, representing a mixed population with
BHs at either side of the gap. The latter mostly consists of binaries where
the primary is above the gap and the secondary below, although, depending on
$\eta$, a small population of binaries where only the secondary is above the
gap can be found as well. These features can be linked to the observed chirp
mass distribution, because the observed distribution (fig. 1) shows a trimodal
distribution very similar to our results: one peak just below $10M_{\odot}$,
one around $28M_{\odot}$, a small group of binaries around $15M_{\odot}$ and a
gap just above $10M_{\odot}$. We also predict a gap just below $20M_{\odot}$,
which is somewhat visible in our AWKDE but even more visible in the estimation
by Abbott et al. (2021b). Our results are not identical to the observed chirp
mass distribution, especially above $20M_{\odot}$ they seem to differ, but the
general shapes of the distributions do seem to agree.

Figure 5: Primary versus secondary mass distributions, for the results of our
estimation (left) and the GW data (right), both in two-dimensional histograms
with a bin-width of $1M_{\odot}$. The left panel shows the upper right panel
from fig. 4. We use the notation $P^{\prime}$ and $S^{\prime}$ here to denote
the stars with the largest and smallest final masses, respectively, in
contrast to $P$ and $S$ from fig. 4, which concern the initial masses. The
right panel shows the cosmologically reweighted posterior distributions of the
GW data (LIGO Scientific Collaboration and Virgo Collaboration, 2021; LIGO
Scientific Collaboration and Virgo Collaboration and KAGRA Collaboration,
2021). Since the posteriors are not normalized, we rescaled them to
$\int_{0M_{\odot}}^{100M_{\odot}}\frac{dN}{dm}dm\approx 1.4\cdot 10^{5}$.
Then, after adding all the posteriors, we normalized the resulting
distribution in between $5M_{\odot}$ and $65M_{\odot}$. We omit GW200322 here,
since no valid posterior for this source has been given. The grey lines are
lines of equal chirp mass, with intervals of $5M_{\odot}$.

In a comparison between our simple estimation and the GW data, we examine the
primary versus secondary mass distribution. Fig. 5 shows this distribution,
together with the posterior distributions of the GW observations. Both the
simulated and the observed distributions clearly show the BH island and the
other failed SNe. Also, there seem to be some observed mergers which can be
linked to the mixed population. There are, however, also differences between
the two distributions. For example, there is a local maximum at
$m_{P^{\prime}}\approx m_{S^{\prime}}\approx 34M_{\odot}$, while the
corresponding maximum in our distribution is situated around
$m_{P^{\prime}}\approx m_{S^{\prime}}\approx 28M_{\odot}$. Also, the observed
distribution has relatively more BHs in the BH island, even though there is an
observational bias towards high-mass mergers. However, we are interested in
the general shape of the distribution and not normalisation and details. After
all, consideration of partial fallback as well as the parameters $\eta$ and
$\zeta$ influence the relative amount of BHs in the BH island. Another
interesting feature of fig. 5 is the fact that the observed mass gap seems to
be situated at lower masses than the simulated mass gap: while we set the
simulated mass gap at $m_{P^{\prime}}\approx m_{S^{\prime}}\approx
17M_{\odot}$, the observed gap appears to be around $m_{P^{\prime}}\approx
m_{S^{\prime}}\approx 14M_{\odot}$. This could mean that we did not have to
shift the results of Schneider et al. (2021) upwards in our remnant function
at all. The difference can be explained by the fact that applying the Gaussian
uncertainty to the individual masses results in a two-dimensional distribution
which is more symmetric than the observed posteriors. Finally, it is
interesting to compare fig. 5 with the rate densities determined by Abbott et
al. (2021b). Where we find three populations: a BH island, other failed SNe
and a mixed population, they find a similar distribution but for NS-NS, BH-BH
and BH-NS mergers, respectively.

##  4 Discussion

We have proposed failed SNe, and in particular the difference between failed
and successful SNe, as a natural explanation for the shape of the BBH mass
distribution. In order to do this, we posed a remnant function (section 2.1),
investigated type II SN rate data (section 2.2), and examined the BBH
distribution caused by this function (section 3). Our model is rather
simplistic, however, and we do not claim that our results describe the actual
binary mass distribution in detail. We are therefore careful not to draw too
strong conclusions in comparing our results from fig. 4 and fig. 5 to the data
as shown in fig. 1. Although, the simplicity of our model does indicate that
our conclusions are a robust feature of BBH distributions caused by remnant
functions similar to ours.
Furthermore, we argue that our remnant function is a more natural hypothetical
explanation for the BBH mass distribution than the gravitational lensing
hypothesis mentioned in the introduction, for two reasons. Firstly, although
they are difficult to compare, we estimate the difference of a factor 50
between the Broadhurst et al. (2018, 2022) merger rate and an optimistic BBH
fraction estimate, discussed in appendix B, to be greater than the difference
between our remnant function and the results of Schneider et al. (2021) and
Marchant et al. (2019). Secondly, the chirp mass distribution is difficult to
explain using the lensing hypothesis. After all, the two stars in the binary
are either both lensed or both non-lensed. This means that there should not
have been any binary where the primary is above the gap and the secondary
below. The lensing hypothesis, therefore, has trouble explaining the trimodal
chirp mass distribution, while our remnant function reproduces it.
As mentioned in the introduction: Schneider et al. (2023) published a paper
which uses a similar method and reaches comparable conclusions, while our
paper was in review. They use simulations equivalent to those of Schneider et
al. (2021), showing a BH island, which result in a bimodal distribution of the
individual BBHs. Their BBH distribution looks similar to the one we find,
except for the PPISN peak since they do not incorporate these. Moreover, they
argue that the bimodal mass distribution causes a trimodal chirp mass
distribution with the same argument we use: the mixing between BHs from the BH
island and the higher mass BHs. They also note that fallback after the SN can
extend the BH island to lower masses, because of the same reason we have
extended our BH island in fig. 2 to include lower masses. Furthermore,
according to Schneider et al. (2023), metallicity mainly affects the case A/B
curve in fig. 2, but because of the mass-loss through wind associated with
metallicity they state that an increase in metallicity can shift the resulting
mass distribution to lower masses. Overall, we find that their results are in
good agreement with our findings, which strengthens our argument.

##  5 Conclusions

After our analysis of the BBH mass distribution and chirp mass distribution
(section 3) resulting from our remnant function (section 2), we conclude the
following:

  * •

The mass distribution of BBHs appears to show a gap for approximately
$14M_{\odot}<m<22M_{\odot}$ (fig. 1). Gravitational lensing as an explanation
for this gap (Broadhurst et al., 2018, 2022) seems, initially, to assume an
unreasonably high merger rate for $z>2$ (fig. 2). We find that this merger
rate implies a BBH fraction which is approximately $50$ times larger than the
one described by Mandel & Farmer (2022). This is improbably large, but not
large enough to completely dismiss this possibility.

  * •

We investigate failed SNe as a more natural explanation for the gap. Firstly,
we approximate the results of the simulations by Schneider et al. (2021) and
Marchant et al. (2019) with a comprehensive remnant function which describes
the relation between ZAMS mass and remnant mass, for case A/B and case C
binaries (fig. 2). This function includes a BH island, which existence is
supported by SN explodability simulations and is the cause of the bimodality
in the BBH mass distribution.

  * •

We investigate whether we can confirm the mass range for failed SNe
observationally. We assume a single mass range for successful SNe, which turn
into failed SNe above a certain mass limit. We take an optimistic value of
$22.5M_{\odot}$ for this mass limit and find that this implies a
$R_{\text{SN}}/\text{SFR}$ reduction of approximately $25\%$, since failed SNe
are not detected and included in the type II SN rate. However, even though
this reduction is compatible with SN survey data to some degree, the data is
quite uncertain and small changes in other model parameters could also account
for a similar reduction in SN rate. We therefore cannot state that the failed
SNe mass range can be confirmed by the SN survey data.

  * •

Using this remnant function, a bimodal BBH distribution can be estimated which
look similar to observation. Also, our simplistic estimation produces a
trimodal chirp mass distribution, which corresponds to observation. The
primary versus secondary mass distribution resulting from our estimation
corresponds to some degree with the posterior samples of the observed GW
sources. This distribution clearly shows a mass gap, and a comparison of
estimation and observation indicates that our remnant function does not need
to deviate from literature values in the degree that it does, strengthening
our conclusions.

Based on our results, we conclude that failed SNe, and in particular the
relation between ZAMS mass and remnant mass which includes a BH island, can
provide a natural explanation for the apparent bimodality in the observed BBH
mass distribution. We therefore state that, when trying to explain the general
shape of the observed BBH distribution, a consideration of stellar evolution
is shown to be fruitful.
Future research could expand or improve multiple aspects of our work.
Reproducing the observed BBH distribution in detail, for instance, would be an
interesting topic of research. Not only could the remnant function perhaps
stay closer to literature values and still produce the desired mass gap, as
implied by fig. 5, but other aspects such as a consideration of partial
fallback, PPISNe and the values of $\eta$ and $\zeta$ also influence the
relative heights of the different peaks. It would be an interesting endeavour
to create a model which includes these aspects in such a way that it
reproduces the observed distribution in more detail. An interesting question
could be if the chirp mass distribution which corresponds to such a detailed
model also resembles observation in detail. And, finally, more (precise) SN
survey and GW data would improve the observational verification, which could
help in identifying the precise effects of failed SNe.

###### Acknowledgements.

This research has made use of data or software obtained from the Gravitational
Wave Open Science Center (gw-openscience.org), a service of LIGO Laboratory,
the LIGO Scientific Collaboration, the Virgo Collaboration, and KAGRA. LIGO
Laboratory and Advanced LIGO are funded by the United States National Science
Foundation (NSF) as well as the Science and Technology Facilities Council
(STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of
Niedersachsen/Germany for support of the construction of Advanced LIGO and
construction and operation of the GEO600 detector. Additional support for
Advanced LIGO was provided by the Australian Research Council. Virgo is
funded, through the European Gravitational Observatory (EGO), by the French
Centre National de Recherche Scientifique (CNRS), the Italian Istituto
Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions
by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan,
Monaco, Poland, Portugal, Spain. The construction and operation of KAGRA are
funded by Ministry of Education, Culture, Sports, Science and Technology
(MEXT), and Japan Society for the Promotion of Science (JSPS), National
Research Foundation (NRF) and Ministry of Science and ICT (MSIT) in Korea,
Academia Sinica (AS) and the Ministry of Science and Technology (MoST) in
Taiwan. This research is supported by the Netherlands Organisation for
Scientific Research (NWO). We thank both Onno Pols for useful discussions and
the anonymous referees who provided comments which helped to improve this
paper.

## References

  * Abbott et al. (2016) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, ApJ, 818, L22
  * Abbott et al. (2021a) Abbott, B. P. et al. 2021a, SoftwareX, 13
  * Abbott et al. (2021b) Abbott, B. P. et al. 2021b, arXiv preprint
  * Abbott et al. (2021c) Abbott, R., Abbott, T. D., Abraham, S., et al. 2021c, ApJ, 913, L7
  * Adams et al. (2017a) Adams, S. M., Kochanek, C. S., Gerke, J. R., & Stanek, K. Z. 2017a, MNRAS, 469, 1445
  * Adams et al. (2017b) Adams, S. M., Kochanek, C. S., Gerke, J. R., Stanek, K. Z., & Dai, X. 2017b, MNRAS, 468, 4968
  * Basinger et al. (2021) Basinger, C. M., Kochanek, C. S., Adams, S. M., Dai, X., & Stanek, K. Z. 2021, MNRAS, 508, 1156
  * Botticella et al. (2017) Botticella, M. T., Cappellaro, E., Greggio, L., et al. 2017, A&A, 598, A50
  * Broadhurst et al. (2018) Broadhurst, T., Diego, J. M., & Smoot, G. F. 2018, arXiv preprint
  * Broadhurst et al. (2022) Broadhurst, T., Diego, J. M., & Smoot, G. F. 2022, arXiv preprint
  * Chabrier (2003) Chabrier, G. 2003, PASP, 115, 763
  * Dorozsmai & Toonen (2022) Dorozsmai, A. & Toonen, S. 2022, arXiv preprint
  * Ertl et al. (2016) Ertl, T., Janka, H.-T., Woosley, S. E., Sukhbold, T., & Ugliano, M. 2016, ApJ, 818, 124
  * Fraley (1968) Fraley, G. S. 1968, Ap&SS, 2, 96
  * Gerke et al. (2015) Gerke, J. R., Kochanek, C. S., & Stanek, K. Z. 2015, MNRAS, 450, 3289
  * Graur et al. (2017) Graur, O., Bianco, F. B., Huang, S., et al. 2017, ApJ, 837, 120
  * Graur et al. (2015) Graur, O., Bianco, F. B., & Modjaz, M. 2015, MNRAS, 450, 905
  * Kochanek (2014a) Kochanek, C. S. 2014a, MNRAS, 446, 1213
  * Kochanek (2014b) Kochanek, C. S. 2014b, ApJ, 785, 28
  * Kochanek et al. (2008) Kochanek, C. S., Beacom, J. F., Kistler, M. D., et al. 2008, ApJ, 684, 1336
  * Kresse et al. (2021) Kresse, D., Ertl, T., & Janka, H.-T. 2021, ApJ, 909, 169
  * LIGO Scientific Collaboration and Virgo Collaboration (2021) LIGO Scientific Collaboration and Virgo Collaboration. 2021, GWTC-2.1: Deep Extended Catalog of Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run - Parameter Estimation Data Release, https://doi.org/10.5281/zenodo.5117702
  * LIGO Scientific Collaboration and Virgo Collaboration and KAGRA Collaboration (2021) LIGO Scientific Collaboration and Virgo Collaboration and KAGRA Collaboration. 2021, GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run — Parameter estimation data release, https://doi.org/10.5281/zenodo.5546663
  * Madau & Dickinson (2014) Madau, P. & Dickinson, M. 2014, ARA&A, 52, 415
  * Mandel & Farmer (2022) Mandel, I. & Farmer, A. 2022, PhR, 955, 1
  * Marchant et al. (2019) Marchant, P., Renzo, M., Farmer, R., et al. 2019, ApJ, 882, 36
  * Mashian & Loeb (2017) Mashian, N. & Loeb, A. 2017, MNRAS, 470, 2611
  * Mazurek (1982) Mazurek, T. 1982, ApJ, 259, L13
  * Müller et al. (2016) Müller, B., Heger, A., Liptai, D., & Cameron, J. B. 2016, MNRAS, 460, 742
  * Neustadt et al. (2021) Neustadt, J. M. M., Kochanek, C. S., Stanek, K. Z., et al. 2021, MNRAS, 508, 516
  * Parravano et al. (2018) Parravano, A., Hollenbach, D., & McKee, C. F. 2018, MNRAS, 480, 2449
  * Salpeter (1955) Salpeter, E. E. 1955, ApJ, 121, 161
  * Schneider et al. (2023) Schneider, F. R. N., Podsiadlowski, P., & Laplace, E. 2023, arXiv preprint
  * Schneider et al. (2021) Schneider, F. R. N., Podsiadlowski, P., & Müller, B. 2021, A&A, 645, A5
  * Shimazaki et al. (2018) Shimazaki, H., Cooper, L. A. D., & Ray, S. 2018, AdaptiveKDE, https://github.com/cooperlab/AdaptiveKDE
  * Shimazaki & Shinomoto (2010) Shimazaki, H. & Shinomoto, S. 2010, Journal of computational neuroscience, 29, 171
  * Sukhbold et al. (2016) Sukhbold, T., Ertl, T., Woosley, S. E., Brown, J. M., & Janka, H.-T. 2016, ApJ, 821, 38
  * Sukhbold & Woosley (2014) Sukhbold, T. & Woosley, S. E. 2014, ApJ, 783, 10
  * Woosley (2017) Woosley, S. E. 2017, ApJ, 836, 244
  * Woosley (2019) Woosley, S. E. 2019, ApJ, 878, 49
  * Woosley et al. (2007) Woosley, S. E., Blinnikov, S., & Heger, A. 2007, Nature, 450, 390

##  Appendix A Asymmetric Gaussians

In order to approximate the posterior chirp mass distribution from Abbott et
al. (2021b) in fig. 1, we use Gaussian probability density functions with two
different standard deviations. This is equivalent to using two halves of
different Gaussians which meet at the mean. These functions are described by
the standard Gaussian formula:

$f^{\pm}(x)=\dfrac{1}{\sigma_{\pm}\sqrt{2\pi}}\exp\left(-\dfrac{1}{2}\left(\dfrac{x-\mu}{\sigma_{\pm}}\right)^{2}\right)\quad,$
(3)

where $f^{-}(x)$ describes the function below the mean $\mu$ with standard
deviation $\sigma_{-}$ and $f^{+}(x)$ describes the function above the mean
with standard deviation $\sigma_{+}$. In other words,
$\sigma_{\pm}=\sigma_{-}$ if $x\leq\mu$ and $\sigma_{\pm}=\sigma_{+}$ if
$x>\mu$. Also, $\overline{\sigma}$ is simply the average standard deviation
$\left(\sigma_{-}+\sigma_{+}\right)/2$. The combination of these halves can
then be constructed as

$f(x)=k_{1}\cdot\left\\{\begin{matrix}f^{-}(x)\hfill&\text{for }x\leq\mu\\\
k_{2}f^{+}(x)&\text{for }x>\mu\end{matrix}\right.\quad,$ (4)

where the constant $k_{2}$ makes sure $f^{-}(\mu)=k_{2}f^{+}(\mu)$ and $k_{1}$
ensures normalization. From this description of $k_{2}$ and eq. 3 it
immediately follows that $k_{2}=f^{-}(\mu)/f^{+}(\mu)=\sigma_{+}/\sigma_{-}$.
Also, the normalization requirement can be used to determine $k_{1}$:
$\int_{-\infty}^{\infty}f(x)dx=1$ gives $k_{1}=\sigma_{-}/\overline{\sigma}$.
Filling in the values of $k_{1}$ and $k_{2}$ gives, then, the overall
asymmetric Gaussian

$f(x)=\dfrac{1}{\overline{\sigma}\sqrt{2\pi}}\exp\left(-\dfrac{1}{2}\left(\dfrac{x-\mu}{\sigma_{\pm}}\right)^{2}\right)\quad,$
(5)

which has the following antiderivative:

$F(x)=\dfrac{\sigma_{\pm}}{2\overline{\sigma}}\erf\left(\dfrac{x-\mu}{\sqrt{2}\sigma_{\pm}}\right)\quad.$
(6)

However, the asymmetry brings with it some ambiguity. When only $\mu$ and the
$90\%$ credible deviations ($L_{-}$ and $L_{+}$) are given, it is not
immediately obvious how to define a 90% interval in equation 5. As a first
option, it is possible to assume that the interval, $\mu-L_{-}\leq
x\leq\mu+L_{+}$, has borders proportional to the standard deviations with a
constant $a$ (i.e. $L_{\pm}=a\sigma_{\pm}$). This means
$\int_{\mu-a\sigma_{-}}^{\mu+a\sigma_{+}}f(x)dx=90\%$, or

$90\%=\erf\left(\dfrac{a}{\sqrt{2}}\right)\implies
a=\sqrt{2}\erf^{-1}\left(90\%\right)\approx 1.65\quad.$ (7)

Alternatively, it is possible to define the 90% interval as "cutting of" 5% on
either side. The relation between the interval borders and the standard
deviations then no longer has a single proportionality constant but two:
$b_{-}$ and $b_{+}$ (i.e. $L_{\pm}=b_{\pm}\sigma_{\pm}$). The fact that
$\int_{-\infty}^{\mu-L_{-}}f(x)dx=\int_{\mu+L_{+}}^{\infty}f(x)dx=5\%$
translates into

$5\%=\dfrac{\sigma_{\pm}}{2\overline{\sigma}}-\dfrac{\sigma_{\pm}}{2\overline{\sigma}}\erf\left(\dfrac{L_{\pm}}{\sqrt{2}\sigma_{\pm}}\right)=\dfrac{1-\erf\left(b_{\pm}/\sqrt{2}\right)}{1+\sigma_{\mp}/\sigma_{\pm}}\quad.$
(8)

This equation can be transformed, after substitution of
$\sigma_{\pm}=L_{\pm}/b_{\pm}$, into the two ($\pm$) equations:

$b_{\pm}=\sqrt{2}\erf^{-1}\left(95\%-5\%\dfrac{b_{\pm}}{b_{\mp}}\dfrac{L_{\mp}}{L_{\pm}}\right)\quad,$
(9)

which can be solved numerically to obtain $b_{-}$ and $b_{+}$. It is clear
that eq. 9 satisfies $\sigma_{-}=\sigma_{+}\rightarrow b_{-}=b_{+}=a$. Fig. 6
shows, similarly to fig. 1, the chirp mass distribution represented by both
methods of asymmetric Gaussians (to avoid a chaotic plot, only a subset of the
data is shown). It is clear that both methods, using either $a$ or $b_{\pm}$,
give almost indistinguishable results for the data-set. In fact, the total
distributions differ less than 5%. Because of this, we use the $a$-method (eq.
7) in this work, which is mathematically more straightforward.

Figure 6: A randomly selected subset of the chirp masses (Abbott et al.,
2021a), represented by asymmetric Gaussians (as in fig. 1) which are similar
to the posterior distributions of Abbott et al. (2021b). The $90\%$-interval
is determined by either the $a$-method from eq. 7 (top panel) or the
$b_{\pm}$-method from eq. 9 (bottom panel).

##  Appendix B Merger rate

According to Broadhurst et al. (2018, 2022), gravitational lensing can explain
the apparent bimodality in the BBH mass distribution, as shown in the top
panel of fig. 1. They note that gravitational lensing focuses the GW signal,
causing the distance of the merger to be underestimated by a factor equal to
the square root of the magnification. This means that the inferred chirp mass
is overestimated. After all, an increased distance means a higher redshift, or
a lower frequency, which gives a higher inferred chirp mass. They state that
the gap in BBH masses can be explained by the hypothesis that a significant
part of the observed BBHs is gravitationally lensed. In order to explain the
observed BBH mass distribution, their model needs a merger rate which has a
high value at high redshift, since there has to be a significant amount of
mergers to be lensed by massive objects at lower redshift. This merger rate,
then, has to decline rapidly towards lower redshift. They pose this
hypothetical merger rate, not as a result from a physical model, but because
from it the desired bimodal BBH distribution follows. They explicitly assume
all mergers with a chirp mass above $20M_{\odot}$ to be gravitationally lensed
(Broadhurst et al., 2022).
The merger rate they use takes the shape of $R_{\text{merger}}(z)=a\cdot
S(z)\cdot F(z)$, where the merger rate ($R_{\text{merger}}$) is given in terms
of a constant $a=6\cdot 10^{3}\ \text{yr}^{-1}\text{Gpc}^{-3}$, as well as a
function $S(z)$ which traces the Madau & Dickinson (2014) SFR and a modulation
function $F(z)$, which equals $1$ for $z>2$ and causes an exponential decline
for $2>z>0.5$. Fig. 7 shows both the merger rate and the Madau & Dickinson
(2014) SFR, together with the SN rate, which we define in section 2.2.
The merger rate seems to be extremely high for $z>2$. In order to investigate
if this is at all possible, we turn to the work of Mandel & Farmer (2022), who
use a “Drake-like” equation to estimate the fraction of binaries that produce
a merging BBH ($f_{\text{BBH}}$). In order to produce a BBH merger, the
following fractions need to be multiplied: the fraction of primary stars which
fall within the SN range ($f_{\text{prim}}$), the fraction of their secondary
stars which fall inside this range ($f_{\text{sec}}$), the fraction of these
binaries which have a separation suitable for merging ($f_{\text{init sep}}$),
the fraction of these binaries which survive the primary SN ($f_{\text{survive
SN1}}$) as well as the common envelope phase ($f_{\text{CE}}$) and the
secondary SN ($f_{\text{survive SN2}}$), and finally the fraction of these
binaries which successfully merge ($f_{\text{merge}}$). In total this gives:

$\displaystyle\begin{split}f_{\text{BBH}}&=f_{\text{prim}}\cdot
f_{\text{sec}}\cdot f_{\text{init sep}}\cdot f_{\text{survive SN1}}\cdot
f_{\text{CE}}\cdot f_{\text{survive SN2}}\cdot f_{\text{merge}}\\\ &\approx
0.001\cdot 0.5\cdot 0.5\cdot 1\cdot 0.1\cdot 1\cdot 0.2=5\cdot
10^{-6}\quad,\end{split}$ (10)

where we filled in the estimations of the fractions made by Mandel & Farmer
(2022).
In order to compare the merger rate inferred by Broadhurst et al. (2018, 2022)
to eq. 1, this merger rate has to be converted to a BBH fraction (which we
will call $f^{\prime}_{\text{BBH}}$), this equals the ratio between the number
of binary stars that merge and the total number of stars. The number of binary
stars that merge is calculated by taking twice the merger rate (since a merger
consists of two stellar remnants). The latter is determined through the SFR.
However, we need to consider the fact that the SFR is determined through the
observed luminosity. We define the binary to have a certain mass ratio
$q=m_{S}/m_{P}$, which means that the total luminosity of the binary can be
described by $L\propto
m_{P}^{3.5}+m_{S}^{3.5}=\left(1+q^{3.5}\right)m_{P}^{3.5}$. In other words,
the SFR implies $1+q^{3.5}$ stars while there are actually $2$. Therefore the
total number of stars has to be adjusted by a factor
$2/\left(1+q^{3.5}\right)$. Here, we make the approximation
$2/\left(1+q^{3.5}\right)\approx 2$, which is valid even if $q$ takes values
as high as $0.6$ or $0.7$.
We use the Chabrier (2003) IMF and compute the IMF integrals from
$m_{\min}=10^{-1}M_{\odot}$ to $m_{\max}=10^{2}M_{\odot}$. The total amount of
stars is, then, simply $\kappa M\int\phi(m)dm$. If we assume a Madau &
Dickinson (2014) SFR, which is the function $S(z)$ multiplied by a constant
$b=1.5\cdot 10^{7}M_{\odot}\text{yr}^{-1}\text{Gpc}^{-3}$, we get for $z>2$:

$\displaystyle\begin{split}f^{\prime}_{\text{BBH}}&=\dfrac{N_{\text{BBH}}}{N_{\text{tot}}}=\dfrac{2R_{\text{merger}}/dt}{2\kappa\psi(z)/dt\int\phi(m)dm}=\dfrac{aS(z)F(z)}{bS(z)}\dfrac{\int
m\phi(m)dm}{\int\phi(m)dm}\\\ &\approx\dfrac{6\cdot 10^{3}}{1.5\cdot
10^{7}}\cdot 0.67=2.7\cdot 10^{-4}\quad,\end{split}$ (11)

where $\psi(z)$ is the SFR at redshift $z$, which means $M=\psi(z)/dt$. We do
not take into account the fact that there is a time difference between the
stellar formation and the merger events, effectively assuming that the bulk of
the mergers occur on a relatively short timescale, compared to the SFR.
In order to analyze $f^{\prime}_{\text{BBH}}$, we use fractions similar to
those used in eq. 10. Leaving out $f_{\text{survive SN1}}$ and
$f_{\text{survive SN2}}$, since they equal $1$, the remaining factors have to
equal $f^{\prime}_{\text{BBH}}=2.7\cdot 10^{-4}$ and account for the deviation
from $f_{\text{BBH}}$ of $f^{\prime}_{\text{BBH}}/f_{\text{BBH}}=56\approx
50$. The question, now, is if this is possible. Using IMF integrals, combined
with the $2/\left(1+q^{3.5}\right)$ correction, we estimate that reasonable
adjustments to $f_{\text{prim}}$ and $f_{\text{sec}}$ could both account for a
factor of approximately $2$. This means that the remaining
$f^{\prime}_{\text{init sep}}$, $f^{\prime}_{\text{CE}}$ and
$f^{\prime}_{\text{merge}}$ have to account for a factor of $56/4\approx 14$,
or an average factor of $14^{1/3}\approx 2.4$ per fraction. This factor is
big, but not large enough to completely dismiss the Broadhurst et al. (2018,
2022) merger rate at $z>2$. After all, the factor of $14$ can be distributed
between the remaining fractions without one of them becoming larger than $1$.
In contrast, if we look at the merger rate at $z<0.5$, the fraction seems to
become quite small. This is due to the factor $F(z)$ in the merger rate. If we
approximate the merger rate below $z=0.5$ as a constant value (specifically as
the value at $z=0.2$), we find that this value equals $F(0.2)\cdot
f^{\prime}_{\text{BBH}}(z>2)\approx 1.4\cdot 10^{-3}\cdot 2.7\cdot
10^{-4}=3.8\cdot 10^{-7}$. This is approximately $8\%$ of the Mandel & Farmer
(2022) fraction (eq. 1). Since $f_{\text{BBH}}$ consists of $7$ fractions,
each has to be adjusted by a factor of $0.08^{1/7}\approx 0.7$, on average,
which is well within their uncertainty.

Figure 7: Madau & Dickinson (2014) star formation rate per $M_{\odot}$, the
corresponding SN rate with $m_{l}=8M_{\odot}$ and $m_{u}=30M_{\odot}$ (through
eq. 1) and the Broadhurst et al. (2018, 2022) BBH merger rate.

##  Appendix C Remnant function

The remnant function we pose in fig. 2 consists of the following equations.
Firstly, we follow Woosley (2019) for the relation between ZAMS mass and the
final mass of the helium core:

$f_{\text{He}}(m)=\left\\{\begin{matrix}0.0385m^{1.603}\hfill&\text{for }m\leq
30M_{\odot}\hfil\\\ 0.50m-5.87M_{\odot}\hfill&\text{for
}m>30M_{\odot}\end{matrix}\right.\quad.$ (12)

We use this relation to rewrite the results of Marchant et al. (2019) as a
function of ZAMS mass instead of helium core mass. This takes the form of two
fourth degree polynomials, shaped as

$f_{\text{PPISN}}(m)=\left\\{\begin{matrix}M_{\odot}\sum_{i=0}^{4}a_{i}\left(\dfrac{m}{M_{\odot}}\right)^{i}&\text{for
}m\leq 151.74M_{\odot}\hfill\\\
M_{\odot}\sum_{i=0}^{4}b_{i}\left(\dfrac{m}{M_{\odot}}\right)^{i}&\text{for
}m>151.74M_{\odot}\hfill\end{matrix}\right.\quad,$ (13)

where the coefficients $a_{i}$ are determined through a polynomial fit through
the data points at $f_{\text{He}}(m)/M_{\odot}=100,115.74,127.74$ and
$151.74$, and $b_{i}$ is given as a fit through the data points at
$f_{\text{He}}(m)/M_{\odot}=151.74,167.74,179.74$ and $189.74$, as given by
Marchant et al. (2019). The two polynomials connect at
$m=f^{-1}_{\text{He}}(70M_{\odot})=151.74M_{\odot}$, with the restriction that
$\frac{d}{dm}f_{\text{PPISN}}(m)=0$ here. The values for the coefficients are,
then:

$\begin{pmatrix}a_{0}\\\ a_{1}\\\ a_{2}\\\ a_{3}\\\
a_{4}\end{pmatrix}=\begin{pmatrix}-3.38751\cdot 10^{2}\\\ 1.01736\cdot
10^{1}\\\ -1.06666\cdot 10^{-1}\\\ 5.18494\cdot 10^{-4}\\\ -9.74407\cdot
10^{-7}\end{pmatrix}\text{ and }\begin{pmatrix}b_{0}\\\ b_{1}\\\ b_{2}\\\
b_{3}\\\ b_{4}\end{pmatrix}=\begin{pmatrix}-3.47418\cdot 10^{4}\\\
8.53251\cdot 10^{2}\\\ -7.84251\\\ 3.20206\cdot 10^{-2}\\\ -4.90175\cdot
10^{-5}\end{pmatrix}\quad,$ (14)

this fit is valid for at least $100M_{\odot}\leq m\leq 192M_{\odot}$. Now, in
order to describe the failed SN population, apart from the BH island, we use
the work of Schneider et al. (2021). Their results entail both a relation
between ZAMS mass and CO core mass and a fit for the failed SN population,
which depends on this relation. For the case A/B mass transfer, this comes
down to:

$f_{\text{FSNe}}^{\text{AB}}(m)=15.832M_{\odot}\exp(0.058\left[\dfrac{3}{20}\dfrac{m}{M_{\odot}}+4\right])-14.784M_{\odot}\quad,$
(15)

where the part in square brackets is a linear fit to the CO core mass relation
and the last term determines the upper border of the mass gap. Similarly, the
failed SN function for case C mass transfer equals

$f_{\text{FSNe}}^{\text{C}}(m)=9.727M_{\odot}\exp(0.041\left[\dfrac{16}{31}\dfrac{m}{M_{\odot}}-\dfrac{128}{31}\right])+4.594M_{\odot}\quad.$
(16)

Eq. 15 and eq. 16 determine the upper border of the mass gap, at
$22M_{\odot}$. Now, for the BH island we simply use a linear function which
produces remnants ranging from $8M_{\odot}$ up to the lower border of the mass
gap, at $14M_{\odot}$. These different functions are valid in ZAMS mass ranges
which are similar to the values used by Schneider et al. (2021) and Marchant
et al. (2019), although (small) deviations from these values should not
influence the general shape of the remnant distribution. For case A/B mass
transfer, this comes down to

$f_{\text{RN}}^{\text{AB}}(m)=\left\\{\begin{matrix}0&\text{for
}m<30.5M_{\odot}\hfill\\\ 1.2m-28.6M_{\odot}&\text{for }30.5M_{\odot}\leq
m\leq 35.5M_{\odot}\hfill\\\ 0&\text{for }35.5M_{\odot}<m<70M_{\odot}\hfill\\\
f_{\text{FSNe}}^{\text{AB}}(m)&\text{for }70M_{\odot}\leq m\leq
100M_{\odot}\hfill\\\ f_{\text{PPISNe}}(m)&\text{for }100M_{\odot}<m\leq
192M_{\odot}\hfill\\\ 0&\text{for
}m>192M_{\odot}\hfill\end{matrix}\right.\quad,$ (17)

where we assume that not only PISNe but also regular successful SNe do not
form remnants, since we are not interested in the distribution of neutron
stars. The function for case C is similar to eq. 17:

$f_{\text{RN}}^{\text{C}}(m)=\left\\{\begin{matrix}0&\text{for
}m<20.5M_{\odot}\hfill\\\ 2.0m-33.0M_{\odot}&\text{for }20.5M_{\odot}\leq
m\leq 23.5M_{\odot}\hfill\\\ 0&\text{for
}23.5M_{\odot}<m<35.5M_{\odot}\hfill\\\
f_{\text{FSNe}}^{\text{C}}(m)&\text{for }35.5M_{\odot}\leq m\leq
58M_{\odot}\hfill\\\ f_{\text{PPISNe}}(m+42M_{\odot})&\text{for
}58M_{\odot}<m\leq 150M_{\odot}\hfill\\\ 0&\text{for
}m>150M_{\odot}\hfill\end{matrix}\right.\quad,$ (18)

where the PPISN function is shifted to lower masses in order to connect with
the failed SN function. Eq. 17 and eq. 18 give the curves which are shown in
fig. 2.

##  Appendix D Model uncertainties

In order to apply a Gaussian uncertainty to the results of our estimation in
fig. 4, we make lognormal fits to the average relative standard deviations of
the GW data (Abbott et al., 2021a). Fig. 8 shows the data together with our
fits. Applying this to our model, we pick a standard deviation randomly from
the lognormal distributions and replace the result by a randomly chosen point
from a Gaussian with the result as mean and the chosen standard deviation.
However, the lognormal fits have a sharp enough peak that taking a single
value for the relative standard deviations, e.g. at the maximum, would not
influence the results significantly.

Figure 8: The relative average standard deviations in the GW data (Abbott et
al., 2021a) for both individual masses (left) and chirp masses (right).

##  Appendix E Mass transfer efficiency and type

Fig. 9 and fig. 10 show the results of our estimation (fig. 4), but for
varying $\eta$ and $\zeta$. The figures show that both parameters do not
greatly influence the final mass distributions. The strongest influence of
these factors concerns the peak at $m\approx 22M_{\odot}$ and the gap at
$\mathcal{M}\approx 19M_{\odot}$, but these effects are not big enough to
influence our conclusions.

Figure 9: Results of the estimation from fig. 4, for varying $\eta$ and
$\zeta=0.5$.

Figure 10: Results of the estimation from fig. 4, for varying $\zeta$ and
$\eta=0.5$.