# PICNN: A Pathway towards Interpretable Convolutional Neural Networks

Wengang Guo, Jiayi Yang, Huilin Yin, Qijun Chen, Wei Ye Corresponding author.

###### Abstract

Convolutional Neural Networks (CNNs) have exhibited great performance in
discriminative feature learning for complex visual tasks. Besides
discrimination power, interpretability is another important yet under-explored
property for CNNs. One difficulty in the CNN interpretability is that filters
and image classes are entangled. In this paper, we introduce a novel pathway
to alleviate the entanglement between filters and image classes. The proposed
pathway groups the filters in a late conv-layer of CNN into class-specific
clusters. Clusters and classes are in a one-to-one relationship. Specifically,
we use the Bernoulli sampling to generate the filter-cluster assignment matrix
from a learnable filter-class correspondence matrix. To enable end-to-end
optimization, we develop a novel reparameterization trick for handling the
non-differentiable Bernoulli sampling. We evaluate the effectiveness of our
method on ten widely used network architectures (including nine CNNs and a
ViT) and five benchmark datasets. Experimental results have demonstrated that
our method PICNN (the combination of standard CNNs with our proposed pathway)
exhibits greater interpretability than standard CNNs while achieving higher or
comparable discrimination power.

## Introduction

The remarkable discrimination power of convolutional neural networks (CNNs)
fosters great applications in numerous tasks. However, in safety-critical
domains such as autonomous vehicles (Zablocki et al. 2022) and healthcare
(D’Amour et al. 2022), interpretability is another crucial property that needs
to be considered.

The interpretability of CNNs has received growing interest in recent research.
Earlier posthoc explanation methods (Zeiler and Fergus 2014; Simonyan,
Vedaldi, and Zisserman 2013; Springenberg et al. 2014; Zhou et al. 2016;
Selvaraju et al. 2017; Bau et al. 2017) focus on generating offline
interpretations such as saliency maps (Simonyan, Vedaldi, and Zisserman 2013)
and class activation mapping (CAM) (Zhou et al. 2016) for well-trained CNNs.
But posthoc methods cannot improve the intrinsic interpretability of CNNs, as
they operate independently from the training process. Recently, the focus of
the community has shifted to train interpretable CNNs. For example, (Zhang,
Wu, and Zhu 2018; Shen et al. 2021) encourage each filter in a late conv-layer
to respond to only one object part. Nevertheless, these methods cannot
directly reveal concepts learned by filters, as they rely on some posthoc
methods (Bau et al. 2017; Zhou et al. 2015) to associate filters to the
predefined concepts.

Figure 1:  Top: comparison of Grad-CAM visualizations using different class-
specific clusters of filters. Bottom: the display of the learned filter-class
correspondence matrix $\mathbf{P}$. The dataset is PASCAL VOC Part (Chen et
al. 2014) with six animal classes, and the target late conv-layer consists of
20 class-specific filters.

In this paper, interpretability is characterized by the alignment between
filters and class-level concepts that are human-understandable. We propose a
novel method to train interpretable CNNs whose filters in the late conv-layer
can directly reveal class-level concepts without the help of any posthoc
methods. We focus on filters in the late conv-layer since these filters encode
class-level concepts rather than low-level primitives (Zeiler and Fergus
2014). (Fong and Vedaldi 2018) has shown that a CNN uses multiple filters in
the late conv-layers to collectively encode a class-level concept. Considering
a target late conv-layer with $N$ filters, there exist $2^{N}$ possible
combinations of filters, allowing these filters to encode up to $2^{N}$ class-
level concepts. This huge concept space prohibits the explanation of the
relations between filters and classes due to their complex many-to-many
correspondences, known as filter-class entanglement. As argued in prior work
(Liang et al. 2020), filter-class entanglement is one of the most critical
reasons that hamper the interpretability of CNNs.

To mitigate filter-class entanglement, we group the $N$ filters into $K$
class-specific clusters, where $K$ is the number of image classes. Clusters
and classes are in a one-to-one correspondence relation, i.e., each class-
specific cluster is encouraged to contain filters that encode the class-level
concepts of only one particular image class. A direct way to group filters is
to predefine a fixed filter-cluster assignment matrix before training.
However, this results in suboptimal performance, as the optimal number of
filters assigned to each class is unknown a priori. (Shen et al. 2021) employs
spectral clustering (Shi and Malik 2000) to group filters during each forward
propagation. Differing from our paper, their goal is to make each filter
represent a set of interpretable image features, such as a specific object
part or image region with a clear meaning.

In this paper, we develop a novel pathway that can group filters into class-
specific clusters and be applied to many network architectures. Specifically,
the pathway introduces a learnable matrix $\mathbf{P}\in\mathbb{R}^{K\times
N}$ whose element $p_{y,i}$ represents the probabilistic correspondence
between the $i$-th filter and $y$-th class. Based on $\mathbf{P}$, we can use
the argmax operation to generate a binary filter-cluster assignment matrix.
However, the argmax operation is non-differentiable, impeding end-to-end
optimization. One potential solution is using Gumbel-Softmax (Jang, Gu, and
Poole 2016) to approximate the argmax operation to enable gradient
propagation. Nevertheless, this approximation introduces bias, as Gumbel-
Softmax constitutes a biased estimate of the original argmax operation. Such
bias can accumulate during the training process and result in suboptimal
performance. To address this challenge, we consider $p_{y,j}$ as the parameter
of the Bernoulli distribution, based on which a filter-cluster assignment
matrix is sampled from $\mathbf{P}$. Since the sampling is also non-
differentiable, we propose a novel reparameterization trick to support end-to-
end training. Our filter-cluster assignment strategy yields superior
performance compared to the Gumbel-Softmax-based strategy. During end-to-end
training, filter grouping and filter learning are jointly optimized: Filter
grouping is conducted in the forward propagation, while filter learning is
conducted in the backward propagation. Figure 1 displays the learned
correspondence matrix $\mathbf{P}$ and Grad-CAM (Selvaraju et al. 2017)
visualizations for validating the learned filter-class correspondence.

The contributions of this paper are as follows:

  * •

We develop a novel pathway to transform a standard CNN into an interpretable
CNN. This pathway is versatile and can be flexibly combined with many CNN
architectures and even Transformer architectures.

  * •

We propose to group filters into class-specific clusters by the filter-cluster
assignment matrix, which is sampled from a learnable filter-class
correspondence matrix $\mathbf{P}$ by the Bernoulli distribution. Each element
of $\mathbf{P}$ can be used as the parameter of the Bernoulli distribution.

  * •

We propose a reparameterization trick for sampling by the Bernoulli
distribution to support end-to-end training.

  * •

We evaluate the effectiveness of our proposed PICNN using five benchmark
datasets. Experimental results demonstrate that PICNN achieves higher or
comparable discrimination power and better interpretability than backbone
models.

## Related Work

We briefly survey the literature on posthoc filter interpretability, learning
interpretable filters, and class-specific filters.

Posthoc filter interpretability has been widely studied, which aims to build
the mapping between the abstract patterns of filters in well-trained CNNs and
the human-understandable domains such as images (Wang et al. 2021). Earlier
works (Zeiler and Fergus 2014; Mahendran and Vedaldi 2015; Yosinski et al.
2015) inspect maximal activations of filters across different images. CAM
(Zhou et al. 2016) localizes image regions that are most important for a
target class. CAM has gained much popularity and promotes numerous further
studies (Selvaraju et al. 2017; Lee et al. 2021; Hasany, Petitjean, and
Mériaudeau 2023; Sarkar et al. 2023). Notably, Grad-CAM (Selvaraju et al.
2017) extends CAM by using gradients to weigh the contribution of each filter
to the target class. Other well-known posthoc methods include guided
backpropagation (Springenberg et al. 2014), saliency maps (Simonyan, Vedaldi,
and Zisserman 2013), and deconvolutional network (Zeiler and Fergus 2014;
Dosovitskiy and Brox 2016). Moreover, some works disentangle the filter
representations of a well-trained CNN into an explanatory graph (Zhang et al.
2018), a decision tree (Zhang et al. 2019), or textual descriptions (Hendricks
et al. 2016; Yang, Kim, and Joo 2022). However, these posthoc methods cannot
remove the existing filter-class entanglement in well-trained CNNs and may not
faithfully capture what the original CNNs compute (Rudin 2019). In contrast,
our work focuses on training interpretable CNNs, where the relations between
filters and classes are clearly revealed by the learned filter-class
correspondence matrix.

Learning interpretable filters has been conducted in ICNN (Zhang, Wu, and Zhu
2018). The learned each filter represents a specific object part, such as
animal eyes. Filters in ICNN could only represent object parts in ball-like
areas. To overcome this limitation, ICCNN (Shen et al. 2021) extends filter
interpretability to image regions with arbitrary shapes. However, the filters
learned in both ICNN and ICCNN are not class-specific, as they are trained in
a class-agnostic manner. Learning class-specific filters is first introduced
by Class-Specific Gate (CSG) (Liang et al. 2020). Specifically, CSG introduces
a binary filter-class gate matrix to represent correspondences between filters
and classes. The binary filter-class gate matrix is relaxed by the $l_{1}$
norm to weigh feature maps of filters and simultaneously forced to be sparse.
Differing from CSG, we consider the filter-class correspondences as
probability parameters and optimize these parameters using a probability
approach.

Class-specific filters have been widely utilized across various studies.
Class-specific filters enable the integration of direct supervision into the
hidden layers of CNNs, thereby alleviating gradient vanishing (Jiang et al.
2017). (Wang, Morariu, and Davis 2018; Martinez et al. 2019) improve fine-
grained recognition using a filter bank that captures class-specific
discriminative patches. Except for discriminative models, generative models
also benefit from class-specific formulation. (Tang et al. 2020; Li et al.
2021) create multiple class-specific generators within a GAN (Goodfellow et
al. 2014) to facilitate the generation of small objects. (Kweon et al. 2021)
proposes a class-specific adversarial erasing framework to generate more
precise CAM. Unlike these models that enhance the discrimination or generation
power of CNNs, we focus on filter interpretability. Additionally, these models
require predefined filter-class correspondences, whereas our work
automatically learns filter-class correspondences.

## Model PICNN

This paper aims to train interpretable CNNs whose filters are disentangled
from image classes and can directly reveal class-level concepts. To achieve
this, we propose to group filters into class-specific clusters. Each class-
specific cluster of filters exclusively activates for inputs from a particular
class and deactivates for inputs from other classes.

### Objective Function

We enrich the target conv-layer in a CNN, typically the last conv-layer, with
an extra pathway, i.e., the interpretation pathway to learn class-specific
filters (Figure 2). We preserve the original discrimination pathway to fit the
underlying task.

1) In the discrimination pathway (the blue shaded area in Figure 2), the
classifier (i.e., the multilayer perceptron (MLP)) receives complete feature
maps from the target conv-layer and outputs a softmax-normalized prediction.
The objective is to fit the underlying task by optimizing a loss function,
such as the cross-entropy loss for a classification task:

$\mathcal{L}_{\text{dis}}=H(\mathbf{y},\mathbf{y}_{1})$ (1)

where $\mathbf{y}\in\mathbb{R}^{K}$ is the one-hot encoding of the class
label111In the paper, integer $y\in\\{1,\cdots,K\\}$ represents the true class
label of an image, while vector $\mathbf{y}\in\mathbb{R}^{K}$ represents the
one-hot encoding of $y$., $\mathbf{y}_{1}\in\mathbb{R}^{K}$ is the prediction
from the discrimination pathway, and $K$ is the number of image classes.

2) In the interpretation pathway (the green shaded area in Figure 2), we input
into the same classifier used in the discrimination pathway the feature maps
of a class-specific cluster of filters. We postulate that each class-specific
cluster of filters specialized for a given class is essential for the network
to classify that class, whereas other clusters can be removed/masked without
compromising discrimination power. Then, the classifier outputs a prediction
$\mathbf{y}_{2}\in\mathbb{R}^{K}$. The cross-entropy loss in this
interpretation pathway is:

$\mathcal{L}_{\text{int}}=H(\mathbf{y},\mathbf{y}_{2})$ (2)

The interpretation pathway is optimized to achieve a comparable classification
performance as the discrimination pathway. The optimization of the
interpretation pathway encourages the learning of the class-specific filters
in the discrimination pathway. We simultaneously optimize the losses of these
two pathways as follows:

$\mathcal{L}=\mathcal{L}_{\text{dis}}+\lambda\mathcal{L}_{\text{int}}$ (3)

where $\lambda$ is a regularization parameter.

### Class-specific Grouping of Filters

#### Filter-class Correspondence Matrix

We introduce a learnable matrix $\mathbf{P}\in\mathbb{R}^{K\times N}$ to
indicate the probabilities of the filter-class correspondences, where $N$ is
the number of filters in the target conv-layer. A larger value of $p_{y,i}$
represents a stronger correspondence between the $i$-th filter and the $y$-th
class. With the correspondence matrix $\mathbf{P}$, we can generate the
filter-cluster assignment matrix $\mathbf{Z}$ by the argmax operation.
However, this operation is non-differentiable. Using Gumbel-Softmax (Jang, Gu,
and Poole 2016) for approximating the argmax operation induces approximation
bias and yields suboptimal performance. Thus, we propose to sample
$\mathbf{Z}$ from $\mathbf{P}$ by the Bernoulli distribution parameterized by
the element of $\mathbf{P}$.

#### Bernoulli Sampling

Supposing the input belongs to the $y$-th class, we conduct $N$ independent
Bernoulli trials to generate a binary filter-cluster assignment vector
$\mathbf{z}=[z_{1},...,z_{N}]\in\\{0,1\\}^{N}$ from
$\mathbf{p}_{y}=[p_{y,1},...,p_{y,N}]$, i.e.,
$z_{i}\sim\mathbf{Ber}(p_{y,i}),i=1,2,\ldots,N$, where $\mathbf{p}_{y}$ is the
$y$-th row of matrix $\mathbf{P}$. The $i$-th element $z_{i}$ determines
whether the $i$-th filter belongs to (if $z_{i}=1$) the class-specific cluster
of the $y$-th class or not (if $z_{i}=0$). After that, we use $\mathbf{z}$ to
mask the complete feature maps $\mathbf{H}=\\{\mathbf{H}_{i}\\}_{i=1}^{N}$
(each $\mathbf{H}_{i}$ is a 2D matrix) of all the filters to generate masked
feature maps
$\widetilde{\mathbf{H}}=\\{\widetilde{\mathbf{H}}_{i}\\}_{i=1}^{N}$ by the
Hadamard product:

$\widetilde{\mathbf{H}}_{i}=\mathbf{H}_{i}\odot z_{i}\quad i=1,\cdots,N$ (4)

where $z_{i}$ is broadcasted along the spatial dimensions of $\mathbf{H}_{i}$
for compatible product operation. If $z_{i}=1$, $\widetilde{\mathbf{H}}_{i}$
is the original feature maps of the $i$-th filter; If $z_{i}=0$,
$\widetilde{\mathbf{H}}_{i}$ is set to a matrix of all zero elements. Finally,
we feed both the complete feature maps $\mathbf{H}$ and the masked ones
$\widetilde{\mathbf{H}}$ into the classifier to compute predictions
$\mathbf{y}_{1}$ and $\mathbf{y}_{2}$, respectively.

### Optimization

Now we elaborate on how to optimize the correspondence matrix $\mathbf{P}$
together with other network parameters using stochastic gradient descent
(SGD). The challenge is that sampling the filter-cluster assignment matrix
$\mathbf{Z}$ from $\mathbf{P}$ by the Bernoulli trials is non-differentiable,
which impedes SGD optimization of $\mathbf{P}$. To deal with this, we resort
to the reparameterization trick (Kingma and Welling 2014).

Figure 2: The pipeline of our model PICNN, which consists of a) the
discrimination pathway to fit the underlying classification task and b) the
interpretation pathway to group filters into class-specific clusters, whose
optimization encourages the discrimination pathway to learn class-specific
filters.

To avoid computing the partial derivative for sampling, we propose to
reparameterize the assignment vector $\mathbf{z}$, which is initially drawn
from the Bernoulli distribution with the probability parameter
$\mathbf{p}_{y}$. Specifically, we redefine $z_{i}\sim\mathbf{Ber}(p_{y,i})$
as:

$z_{i}=p_{y,i}+\eta_{i}$ (5)

where $\eta_{i}$ is a normalization scalar defined as:

$\displaystyle\eta_{i}=\begin{cases}1-p_{y,i},&p_{y,i}\geq\epsilon_{i}\\\
-p_{y,i},&p_{y,i}<\epsilon_{i}\end{cases}\quad i=1,\cdots,N$ (6)

where $\epsilon_{i}$ is sampled from the continuous uniform distribution
$\mathbf{\mathcal{U}}(0,1)$, i.e.,
$\epsilon_{i}\sim\mathbf{\mathcal{U}}(0,1)$. This reparameterization trick
converts the sampling of $z_{i}$ from a Bernoulli distribution, which relies
on trainable parameters, into sampling $\epsilon_{i}$ from a fixed uniform
distribution. This trick enables us to compute the partial derivative
$\frac{\partial\mathbf{z}}{\partial\mathbf{P}}$.

######  Theorem 1.

$z_{i}\sim\mathbf{Ber}(p_{y,i})$.

###### Proof.

Substituting Equation 6 into Equation 5, we obtain:

$\displaystyle z_{i}=\begin{cases}1,&p_{y,i}\geq\epsilon_{i}\\\
0,&p_{y,i}<\epsilon_{i}\end{cases}\quad i=1,\cdots,N$ (7)

Since $\epsilon_{i}\sim\mathbf{\mathcal{U}}(0,1)$, the probability
$p(p_{y,i}\geq\epsilon_{i})$ equals the ratio between the length of the
interval $[0,p_{y,i}]$ and the total interval $[0,1]$. Thus, we have
$p(z_{i}=1)=p(p_{y,i}\geq\epsilon_{i})=p_{y,i}$ and
$p(z_{i}=0)=1-p(z_{i}=1)=1-p_{y,i}$, meaning $z_{i}$ follows the Bernoulli
distribution $\mathbf{Ber}(p_{y,i})$. ∎

Our model may suffer from a trivial solution, as the prediction
$\mathbf{y}_{2}$ is indirectly dependent on the label $y$ (from which a red
arrow flows in Figure 3). As a result, the network can exploit this dependency
and effortlessly minimize the loss $\mathcal{L}_{\text{int}}$ (Figure 4). To
prevent this, we propose to replace the label $y$ with a pseudo-label
$\widetilde{y}$ for indexing the correspondence matrix $\mathbf{P}$. The
pseudo-label $\widetilde{y}$ is drawn from a categorical distribution
parameterized by the softmax-normalized prediction $\mathbf{y}_{1}$ from the
discrimination pathway, i.e.,
$\widetilde{y}\sim\mathbf{Cat}(K,\mathbf{y}_{1})$. This intuitive idea is that
the pseudo-label $\widetilde{y}$ gradually approaches the label $y$ as the
training unfolds. To support SGD optimization, we utilize the
reparameterization trick once again to redefine $\widetilde{y}$. Specifically,
we reparameterize the pseudo-label as:

$\widetilde{\mathbf{y}}=\mathbf{y}_{1}+\boldsymbol{\tau}$ (8)

where $\boldsymbol{\tau}$ is a constant vector defined as:

$\displaystyle\boldsymbol{\tau}[k]=\begin{cases}1-\mathbf{y}_{1}[k],&\sum\limits_{j=0}^{k-1}\mathbf{y}_{1}[j]<\xi\leq\sum\limits_{j=0}^{k}\mathbf{y}_{1}[j]\\\
-\mathbf{y}_{1}[k],&\text{otherwise}\end{cases}$ (9) $\displaystyle
k=1,\cdots,K$

where $\xi\sim\boldsymbol{\mathcal{U}}(0,1)$, $[k]$ indexes the $k$-th entry
of a vector and we append 0 before the first element of $\mathbf{y}_{1}$ and 1
after the last element of $\mathbf{y}_{1}$.

######  Theorem 2.

$\widetilde{y}\sim\mathbf{Cat}(K,\mathbf{y}_{1})$.

###### Proof.

Substituting Equation 9 into Equation 8, we obtain:

$\displaystyle\widetilde{\mathbf{y}}[k]=\begin{cases}1,&\sum\limits_{j=0}^{k-1}\mathbf{y}_{1}[j]<\xi\leq\sum\limits_{j=0}^{k}\mathbf{y}_{1}[j]\\\
0,&\text{otherwise}\end{cases}$ (10)

Since $\xi\sim\mathbf{\mathcal{U}}(0,1)$, the probability
$p(\sum_{i=0}^{k-1}\mathbf{y}_{1}[i]<\xi\leq\sum_{i=0}^{k}\mathbf{y}_{1}[i])$
equals to the ratio between the length of the interval
$[\sum_{i=0}^{k-1}\mathbf{y}_{1}[i],\sum_{i=0}^{k}\mathbf{y}_{1}[i]]$ and the
total interval $[0,1]$. Thus, we have
$p(\widetilde{{y}}=k)=\mathbf{y}_{1}[k]$, i.e.,
$\widetilde{{y}}\sim\mathbf{Cat}(K,\mathbf{y}_{1})$.

∎

Figure 3: The interpretation pathway encounters a label leakage issue, as the
label $y$ is employed to index the correspondence matrix $\mathbf{P}$. To deal
with this issue, we replace $y$ with $\widetilde{y}$ sampled from the
categorical distribution parameterized by the softmax-normalized prediction of
the $y$-th class in the discrimination pathway.

Figure 4:  Trivial solution. a) Using the label $y$ as index vector, the loss
$\mathcal{L}_{\text{int}}$ quickly becomes zero, leading to a trivial
solution. b) After replacing the label $y$ with the pseudo-label
$\widetilde{y}$, our model shows a stable evolution.

Figure 4 demonstrates that this categorical sampling effectively circumvents
the trivial solution.

### Time Complexity Analysis

Most clustering methods incur a high time complexity between
$\mathcal{O}(M^{2})$ and $\mathcal{O}(M^{3})$ if the number of input images is
$M$. In PICNN, the interpretation pathway plays a role that is equivalent to a
clustering component but with a lower time complexity of
$\mathcal{O}(M(2K+2K^{2}+2KN+Nd^{2}+2KN^{2}))$, where $K$ is the number of
image classes, $N$ is the number of filters and $d$ is the size of the feature
map. And the space complexity is $\mathcal{O}(NK)$. We report the running time
of PICNN for large-scale architectures in Table 3. PICNN adds at most 5%
computational overhead to backbone models.

|  | CIFAR-10 | PASCAL VOC Part | STL-10
---|---|---|---|---
Backbone | Method | ACC1$\uparrow$ | ACC2$\uparrow$ | ACC3$\downarrow$ | MIS$\uparrow$ | ACC1$\uparrow$ | ACC2$\uparrow$ | ACC3$\downarrow$ | MIS$\uparrow$ | ACC1$\uparrow$ | ACC2$\uparrow$ | ACC3$\downarrow$ | MIS$\uparrow$
| STD | 0.952 | 0.588 | 0.950 | 0.191 | 0.782 | 0.769 | 0.781 | 0.203 | 0.895 | 0.778 | 0.894 | 0.148
| ICCNN | 0.955 | N/A | N/A | 0.067 | 0.754 | N/A | N/A | 0.140 | 0.871 | N/A | N/A | 0.052
| CSG | 0.830 | 0.505 | 0.659 | 0.132 | 0.781 | 0.686 | 0.666 | 0.148 | 0.888 | 0.373 | 0.699 | 0.115
ResNet-18 | PICNN | 0.951 | 0.948 | 0.333 | 0.243 | 0.800 | 0.798 | 0.639 | 0.246 | 0.896 | 0.893 | 0.657 | 0.224
| STD | 0.918 | 0.401 | 0.913 | 0.128 | 0.803 | 0.629 | 0.799 | 0.138 | 0.702 | 0.684 | 0.697 | 0.142
VGG-11 | PICNN | 0.915 | 0.910 | 0.192 | 0.177 | 0.840 | 0.837 | 0.346 | 0.166 | 0.716 | 0.709 | 0.316 | 0.131
| STD | 0.904 | 0.308 | 0.893 | 0.091 | 0.838 | 0.733 | 0.837 | 0.122 | 0.646 | 0.618 | 0.644 | 0.115
AlexNet | PICNN | 0.905 | 0.900 | 0.188 | 0.191 | 0.821 | 0.817 | 0.432 | 0.169 | 0.645 | 0.641 | 0.205 | 0.118
| STD | 0.958 | 0.438 | 0.954 | 0.163 | 0.918 | 0.917 | 0.714 | 0.148 | 0.829 | 0.791 | 0.824 | 0.200
DenseNet-121 | PICNN | 0.957 | 0.957 | 0.462 | 0.257 | 0.917 | 0.915 | 0.763 | 0.229 | 0.832 | 0.831 | 0.734 | 0.266
| STD | 0.950 | 0.540 | 0.948 | 0.208 | 0.930 | 0.844 | 0.928 | 0.176 | 0.829 | 0.818 | 0.828 | 0.210
MobileNetV2 | PICNN | 0.956 | 0.952 | 0.403 | 0.261 | 0.928 | 0.925 | 0.844 | 0.251 | 0.839 | 0.834 | 0.746 | 0.258
| STD | 0.955 | 0.351 | 0.951 | 0.174 | 0.942 | 0.889 | 0.940 | 0.165 | 0.836 | 0.820 | 0.828 | 0.229
EfficientNet-B0 | PICNN | 0.956 | 0.955 | 0.380 | 0.253 | 0.940 | 0.937 | 0.873 | 0.253 | 0.842 | 0.833 | 0.764 | 0.275
Table 1: Experimental results on three benchmark datasets. We conduct
comparison experiments between standard methods and our method on six CNN
architectures. STD stands for standard CNNs, ICCNN stands for Compositional
CNN, CSG stands for Class-Specific Gate CNN, and PICNN is our method. A higher
value is better for ACC1, ACC2, and MIS, while a lower value is better for
ACC3.

## Experiments

### Experimental Settings

#### Evaluation Metric

We use classification accuracy and mutual information score (MIS) (Liang et
al. 2020) as evaluation metrics to assess the discrimination power and
interpretability. We report the following four metrics: 1) classification
accuracy of the discrimination pathway (ACC1$\uparrow$). ACC1 evaluates the
discrimination power. 2) classification accuracy of the interpretation pathway
(ACC2$\uparrow$). Since the input to the interpretation pathway is the feature
maps of a class-specific cluster of filters, ACC2 is consistently lower than
ACC1. 3) classification accuracy using all the filters except a class-specific
cluster of filters (ACC3$\downarrow$). Since the filters used in the
computation of ACC3 and ACC2 complement each other, the lower ACC3 is, the
better the model. 4) MIS (MIS $\uparrow$) measures the mutual information
between filter activation and prediction on classes. MIS is defined as
$\text{MIS}=\mathbf{mean}_{i}(\mathbf{max}_{y}(m_{yi}))$, where
${m}_{yi}=\mathbf{MI}(\mathbf{H}_{i}\|\mathbf{y})$. The metrics ACC2, ACC3,
and MIS quantify the degree of filter-class entanglement from different
aspects, thus they can evaluate the interpretability.

#### Datasets and Types of Backbone CNNs

We use three benchmark classification datasets in Table 1, including CIFAR-10
(Krizhevsky, Hinton et al. 2009), STL-10 (Coates, Ng, and Lee 2011), and
PASCAL VOC Part (Chen et al. 2014). CIFAR-10 consists of 50,000 training
images and 10,000 test images in 10 classes. STL-10 contains 5,000 training
images and 8,000 test images in 10 classes. Following (Liang et al. 2020), we
select six animal classes from PASCAL VOC Part with a 70%/30% training/test
split. To further evaluate the efficacy and effectiveness of PICNN on large
datasets with more classes, we use CIFAR-100 (Krizhevsky, Hinton et al. 2009)
and TinyImageNet (Deng et al. 2009) in Table 2. Like CIFAR-10, CIFAR-100 also
consists of 50,000 training images and 10,000 test images evenly distributed
into 100 classes. TinyImageNet is a scaled-down version of the original
ImageNet involving 200 classes with a 10:1 ratio of training images to test
images. We use the official training/test data split, except for PASCAL VOC
Part.

We combine our interpretation pathway with six typical CNN architectures,
three large CNN architectures, and a transformer architecture (Dosovitskiy et
al. 2020) (ViT-b-12) to make them interpretable. The six CNN architectures are
VGG-11 (Simonyan and Zisserman 2014), AlexNet (Krizhevsky 2014), ResNet-18 (He
et al. 2016), DenseNet-121 (Huang et al. 2017), MobileNetV2 (Sandler et al.
2018), and EfficientNet-B0 (Tan and Le 2019). The three deeper and wider CNN
architectures include ResNet-50 (He et al. 2016), ResNet-152 (He et al. 2016),
and Wide-ResNet (Zerhouni et al. 2017).

#### Implementation Details

Our code is based on the PyTorch (Paszke et al. 2019) toolbox and publicly
available at Github222https://github.com/spdj2271/PICNN. We make filters in
the target conv-layer class-specific by optimizing our loss function. The
regularization parameter $\lambda$ is set to 2 and the effect of the $\lambda$
values is discussed later. Other default settings include: a batch size of
128; the Adam optimizer with an initial learning rate of 0.001; pretrained
weights from ImageNet (Deng et al. 2009); and a total of 200 training epochs.
All metrics presented in this paper are computed on the test set. The
experiments are carried out on a server with an Xeon(R) Platinum 8352V CPU and
one Nvidia RTX 4090 GPU.

### Experimental Results

#### Comparison with Existing Methods

We first compare our PICNN with a Standard CNN (STD) (ResNet-18, $\lambda=0$),
Interpretable Compositional CNNs (ICCNN) (Shen et al. 2021) and Class-Specific
Gate (CSG) (Liang et al. 2020). All the networks use the same backbone
ResNet-18. For STD, we randomly initialize the correspondence matrix
$\mathbf{P}$ and use our Bernoulli sampling to generate the filter-cluster
assignment matrix. As shown in the upper part of Table 1, PICNN significantly
outperforms the comparison methods in terms of ACC2, ACC3, and MIS across all
datasets. Since the filters learned by ICCNN are not class-specific, ACC2 and
ACC3 are not suitable to evaluate its performance. The corresponding result is
denoted as N/A. Additionally, PICNN is better than STD on PASCAL VOC Part and
STL-10 datasets in terms of ACC1 and achieves comparable results to STD on
CIFAR-10. CSG shows a noticeable decrease in ACC1 on CIFAR-10, and ICCNN shows
relatively poorer performance on PASCAL VOC Part and STL-10. The results on
large datasets are shown in Table 2 and demonstrate once again that PICNN has
a significant improvement in ACC2, ACC3, and MIS over STD.

All these results indicate that PICNN achieves better interpretability and
maintains comparable discrimination power simultaneously. In contrast, the
comparison methods are unsatisfactory in both the discrimination power and
interpretability.

Dataset | Class | Model | ACC1$\uparrow$ | ACC2$\uparrow$ | ACC3$\downarrow$ | MIS$\uparrow$
---|---|---|---|---|---|---
CIFAR-100 | 100 | STD | 0.755 | 0.039 | 0.753 | 0.020
PICNN | 0.724 | 0.707 | 0.131 | 0.026
TinyImageNet | 200 | STD | 0.619 | 0.010 | 0.619 | 0.010
PICNN | 0.555 | 0.514 | 0.082 | 0.017
Table 2: Performance on large multi-class datasets with ResNet-18 as the
backbone.

#### Various Network Architectures

One advantage of our method lies in its versatility to combine with various
network architectures. As shown in Table 1 and Table 3, our novel pathway
significantly boosts the interpretability of all the ten network
architectures, which include six typical CNN architectures, three large CNN
architectures, and a transformer architecture. Except for using DenseNet-121
as the backbone on the PASCAL VOC Part dataset, PICNN has the best performance
in terms of all three interpretability metrics ACC2, ACC3, and MIS. Compared
with STD, the improvement is significant. PICNN does not work well when
DenseNet-121 is the backbone. One reason might be: each layer in DenseNet-121
is connected to all other layers. As a result, the target conv-layer may
include too many low-level concepts from early conv-layers, which leads to
difficulties in grouping filters into class-specific clusters. We can also see
that PICNN achieves higher ACC1 values in 14 out of 21 experiments. Besides,
the performance of PICNN with ViT-b-12 as the backbone demonstrates that our
approach is not only applicable to CNN architectures but also to transformer
architectures.

Backbone | Metric | ACC1$\uparrow$ | ACC2$\uparrow$ | ACC3$\downarrow$ | MIS$\uparrow$ | Train |  Infer
---|---|---|---|---|---|---|---
ResNet | STD | 0.951 | 0.457 | 0.948 | 0.115 | 18302 | 1564
-50 | PICNN | 0.955 | 0.950 | 0.318 | 0.231 | 18484 | 1592
ResNet | STD | 0.957 | 0.481 | 0.954 | 0.150 | 27339 | 2403
-152 | PICNN | 0.958 | 0.954 | 0.309 | 0.230 | 27377 | 2514
Wide- | STD | 0.940 | 0.295 | 0.930 | 0.117 | 76572 | 4733
ResNet | PICNN | 0.940 | 0.933 | 0.137 | 0.204 | 77212 | 4831
ViT | STD | 0.799 | 0.242 | 0.794 | 0.066 | 37245 | 3234
-b-12 | PICNN | 0.802 | 0.746 | 0.039 | 0.096 | 37615 | 3257
Table 3: Evaluation for large CNNs and Vision Transformer architecture on the
CIFAR-10 dataset. The final two columns display the training and inference
time (in seconds).

#### Effectiveness of the Bernoulli Sampling

Instead of using our Bernoulli sampling to generate the filter-cluster
assignment matrix, we can apply Gumbel-Softmax to the correspondence matrix
$\mathbf{P}$. We conducted an ablation study using Gumbel-Softmax (temperature
is set to 0.01) to replace the Bernoulli sampling in PICNN. It can be seen
from Table 4 that the results of the Bernoulli sampling are better than those
of Gumbel-Softmax, especially using the metric ACC3. One possible explanation
is that Gumbel-Softmax introduces random Gumbel noise to approximate the
argmax operation for enabling gradient backpropagation. The approximation
affects the results.

Model | ACC1$\uparrow$ | ACC2$\uparrow$ | ACC3$\downarrow$ | MIS$\uparrow$
---|---|---|---|---
Gumbel-Softmax | 0.950 | 0.946 | 0.910 | 0.240
Bernoulli-Samping | 0.951 | 0.948 | 0.333 | 0.243
Table 4: Performance comparison of PICNN using different strategies for
generating filter-cluster assignment matrix on the CIFAR-10 dataset. ResNet-18
is used as the backbone of PICNN.

#### Effect of the filter-to-class ratio

For a target conv-layer in a given CNN, we set the filter-to-class ratio to
12.8. This means that for CIFAR-10, we choose the number of filters to be 128.
To exploit the effect of the filter-to-class ratio $r$, we vary $r=N/K$ on the
CIFAR-10 dataset, where $N$ is the number of filters in the target conv-layer
and $K$ is the number of image classes. As shown in Figure 5, PICNN shows
consistent advantages over STD (ResNet-18) across different $r$ values. This
indicates that PICNN is robust to the filter-to-class ratio.

Figure 5: Effect of the filter-to-class ratio $r={N}/{K}$ on the performance.
The backbone (STD) of PICNN is ResNet-18.

#### Effect of Regularization Parameter $\lambda$

Figure 6 plots the curves of ACC1, ACC2, ACC3, and MIS of PICNN using
ResNet-18 as the backbone on the CIFAR-10 dataset when varying the values of
$\lambda$ from 0 to 10. ACC3 decreases steadily, MIS increases fluctuately,
while ACC1 and ACC2 remain stable as $\lambda$ increases. This implies that
with the increase of interpretation pathway weights, the filters at the target
conv-layer become more class-specific and the discrimination power of the
discrimination pathway is preserved.

Figure 6:  Effect of regularization parameter $\lambda$ of PICNN using
ResNet-18 as the backbone on the CIFAR-10 dataset. PICNN works well across a
wide range of $\lambda$.  Figure 7:  The Grad-CAM visualizations of filters on
the PASCAL VOC Part dataset learned by STD (ResNet-18), CSG, and PICNN. CSG
and PICNN use the same STD (ResNet-18) as the backbone. In each method, the
(a) column represents CAMs using all filters, (b) column represents CAMs using
a cluster of class-specific filters, and (c) column represents CAMs using the
complementary clusters of class-specific filters. $N_{p}$ is the number of
pictures in a class and $N_{f}$ is the number of class-specific filters
assigned to a class.

#### Visualization

Figure 7 displays the Grad-CAM visualizations on the PASCAL VOC Part dataset
using three different sets of filters: (a) all filters, (b) a cluster of
class-specific filters, and (c) the complementary clusters of class-specific
filters. As above, we randomly initialize the correspondence matrix
$\mathbf{P}$ and use our Bernoulli sampling to generate the filter-cluster
assignment matrix for STD (ResNet-18). We observe few differences among these
three Grad-CAM visualizations in STD. This phenomenon proves that the filters
and classes are entangled in STD. The Grad-CAM visualizations using the second
set of filters (column (b)) of PICNN sometimes capture more class-specific
information than using the first set of filters (column (a)), such as
highlighting more parts of the cat’s body and the horse’s head. The Grad-CAM
visualizations using the third set of filters (column (c)) of PICNN mostly
highlight unimportant image regions, such as the meadow and road. We can see
that the Grad-CAM visualizations using the second and third sets of filters of
PICNN differ a lot. But we do not see a significant difference in STD and CSG.

## Conclusion

In this paper, we have proposed a novel pathway to transform a standard CNN
into an interpretable CNN without compromising its high discrimination power.
The proposed pathway uses the Bernoulli sampling to generate the filter-
cluster assignment matrix from a learnable filter-class correspondence matrix.
The filter-cluster assignment matrix groups the filters in the target late
conv-layer in CNN into class-specific clusters and thus mitigates the filter-
class entanglement problem. Because the Bernoulli sampling is non-
differentiable, we propose a reparameterization trick for end-to-end learning.
Experiments have shown that our method PICNN is superior to standard CNNs in
terms of both interpretability and discrimination power. Moreover, our pathway
has good versatility and can be combined with various network architectures.

## Acknowledgments

We thank the anonymous reviewers for their valuable and constructive comments.
This work was supported in part by the National Key Research and Development
Program of China under Grant 2020AAA0108100.

## References

  * Bau et al. (2017) Bau, D.; Zhou, B.; Khosla, A.; Oliva, A.; and Torralba, A. 2017.  Network dissection: Quantifying interpretability of deep visual representations.  In _CVPR_.
  * Chen et al. (2014) Chen, X.; Mottaghi, R.; Liu, X.; Fidler, S.; Urtasun, R.; and Yuille, A. 2014.  Detect what you can: Detecting and representing objects using holistic models and body parts.  In _CVPR_.
  * Coates, Ng, and Lee (2011) Coates, A.; Ng, A.; and Lee, H. 2011.  An analysis of single-layer networks in unsupervised feature learning.  In _AISTATS_.
  * D’Amour et al. (2022) D’Amour, A.; Heller, K.; Moldovan, D.; Adlam, B.; Alipanahi, B.; Beutel, A.; Chen, C.; Deaton, J.; Eisenstein, J.; Hoffman, M. D.; et al. 2022.  Underspecification presents challenges for credibility in modern machine learning.  _JMLR_.
  * Deng et al. (2009) Deng, J.; Dong, W.; Socher, R.; Li, L.-J.; Li, K.; and Fei-Fei, L. 2009.  Imagenet: A large-scale hierarchical image database.  In _CVPR_.
  * Dosovitskiy et al. (2020) Dosovitskiy, A.; Beyer, L.; Kolesnikov, A.; Weissenborn, D.; Zhai, X.; Unterthiner, T.; Dehghani, M.; Minderer, M.; Heigold, G.; Gelly, S.; et al. 2020\.  An image is worth 16x16 words: Transformers for image recognition at scale.
  * Dosovitskiy and Brox (2016) Dosovitskiy, A.; and Brox, T. 2016.  Inverting visual representations with convolutional networks.  In _CVPR_.
  * Fong and Vedaldi (2018) Fong, R.; and Vedaldi, A. 2018.  Net2vec: Quantifying and explaining how concepts are encoded by filters in deep neural networks.  In _CVPR_.
  * Goodfellow et al. (2014) Goodfellow, I.; Pouget-Abadie, J.; Mirza, M.; Xu, B.; Warde-Farley, D.; Ozair, S.; Courville, A.; and Bengio, Y. 2014.  Generative adversarial nets.  In _NeurIPS_.
  * Hasany, Petitjean, and Mériaudeau (2023) Hasany, S. N.; Petitjean, C.; and Mériaudeau, F. 2023.  Seg-XRes-CAM: Explaining Spatially Local Regions in Image Segmentation.  In _CVPR_.
  * He et al. (2016) He, K.; Zhang, X.; Ren, S.; and Sun, J. 2016.  Deep residual learning for image recognition.  In _CVPR_.
  * Hendricks et al. (2016) Hendricks, L. A.; Akata, Z.; Rohrbach, M.; Donahue, J.; Schiele, B.; and Darrell, T. 2016.  Generating visual explanations.  In _ECCV_.
  * Huang et al. (2017) Huang, G.; Liu, Z.; Van Der Maaten, L.; and Weinberger, K. Q. 2017.  Densely connected convolutional networks.  In _CVPR_.
  * Jang, Gu, and Poole (2016) Jang, E.; Gu, S.; and Poole, B. 2016.  Categorical Reparameterization with Gumbel-Softmax.  In _ICLR_.
  * Jiang et al. (2017) Jiang, Z.; Wang, Y.; Davis, L.; Andrews, W.; and Rozgic, V. 2017.  Learning discriminative features via label consistent neural network.  In _WACV_.
  * Kingma and Welling (2014) Kingma, D. P.; and Welling, M. 2014.  Auto-Encoding Variational Bayes.  In _ICLR_.
  * Krizhevsky (2014) Krizhevsky, A. 2014.  One weird trick for parallelizing convolutional neural networks.  _CoRR_.
  * Krizhevsky, Hinton et al. (2009) Krizhevsky, A.; Hinton, G.; et al. 2009.  Learning multiple layers of features from tiny images.
  * Kweon et al. (2021) Kweon, H.; Yoon, S.-H.; Kim, H.; Park, D.; and Yoon, K.-J. 2021.  Unlocking the potential of ordinary classifier: Class-specific adversarial erasing framework for weakly supervised semantic segmentation.  In _ICCV_.
  * Lee et al. (2021) Lee, J. R.; Kim, S.; Park, I.; Eo, T.; and Hwang, D. 2021.  Relevance-cam: Your model already knows where to look.  In _CVPR_.
  * Li et al. (2021) Li, Y.; Li, Y.; Lu, J.; Shechtman, E.; Lee, Y. J.; and Singh, K. K. 2021.  Collaging class-specific gans for semantic image synthesis.  In _ICCV_.
  * Liang et al. (2020) Liang, H.; Ouyang, Z.; Zeng, Y.; Su, H.; He, Z.; Xia, S.-T.; Zhu, J.; and Zhang, B. 2020.  Training interpretable convolutional neural networks by differentiating class-specific filters.  In _ECCV_.
  * Mahendran and Vedaldi (2015) Mahendran, A.; and Vedaldi, A. 2015.  Understanding deep image representations by inverting them.  In _CVPR_.
  * Martinez et al. (2019) Martinez, B.; Modolo, D.; Xiong, Y.; and Tighe, J. 2019.  Action recognition with spatial-temporal discriminative filter banks.  In _ICCV_.
  * Paszke et al. (2019) Paszke, A.; Gross, S.; Massa, F.; Lerer, A.; Bradbury, J.; Chanan, G.; Killeen, T.; Lin, Z.; Gimelshein, N.; Antiga, L.; et al. 2019.  Pytorch: An imperative style, high-performance deep learning library.  In _NeurIPS_.
  * Rudin (2019) Rudin, C. 2019.  Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead.  _Nat. Mac. Intell._
  * Sandler et al. (2018) Sandler, M.; Howard, A.; Zhu, M.; Zhmoginov, A.; and Chen, L.-C. 2018.  Mobilenetv2: Inverted residuals and linear bottlenecks.  In _CVPR_.
  * Sarkar et al. (2023) Sarkar, S.; Babu, A. R.; Gundecha, V.; Guillen, A.; Mousavi, S.; Luna, R.; Ghorbanpour, S.; and Naug, A. 2023.  RL-CAM: Visual Explanations for Convolutional Networks Using Reinforcement Learning.  In _CVPR_.
  * Selvaraju et al. (2017) Selvaraju, R. R.; Cogswell, M.; Das, A.; Vedantam, R.; Parikh, D.; and Batra, D. 2017.  Grad-cam: Visual explanations from deep networks via gradient-based localization.  In _ICCV_.
  * Shen et al. (2021) Shen, W.; Wei, Z.; Huang, S.; Zhang, B.; Fan, J.; Zhao, P.; and Zhang, Q. 2021.  Interpretable compositional convolutional neural networks.  In _IJCAI_.
  * Shi and Malik (2000) Shi, J.; and Malik, J. 2000.  Normalized cuts and image segmentation.  _TPAMI_.
  * Simonyan, Vedaldi, and Zisserman (2013) Simonyan, K.; Vedaldi, A.; and Zisserman, A. 2013.  Deep inside convolutional networks: Visualising image classification models and saliency maps.  In _ICLR_.
  * Simonyan and Zisserman (2014) Simonyan, K.; and Zisserman, A. 2014.  Very deep convolutional networks for large-scale image recognition.  In _ICLR_.
  * Springenberg et al. (2014) Springenberg, J. T.; Dosovitskiy, A.; Brox, T.; and Riedmiller, M. 2014.  Striving for simplicity: The all convolutional net.  In _ICLR_.
  * Tan and Le (2019) Tan, M.; and Le, Q. 2019.  Efficientnet: Rethinking model scaling for convolutional neural networks.  In _ICML_.
  * Tang et al. (2020) Tang, H.; Xu, D.; Yan, Y.; Torr, P. H.; and Sebe, N. 2020.  Local class-specific and global image-level generative adversarial networks for semantic-guided scene generation.  In _CVPR_.
  * Wang et al. (2021) Wang, J.; Liu, H.; Wang, X.; and Jing, L. 2021.  Interpretable image recognition by constructing transparent embedding space.  In _ICCV_.
  * Wang, Morariu, and Davis (2018) Wang, Y.; Morariu, V. I.; and Davis, L. S. 2018.  Learning a discriminative filter bank within a cnn for fine-grained recognition.  In _CVPR_.
  * Yang, Kim, and Joo (2022) Yang, Y.; Kim, S.; and Joo, J. 2022.  Explaining deep convolutional neural networks via latent visual-semantic filter attention.  In _CVPR_.
  * Yosinski et al. (2015) Yosinski, J.; Clune, J.; Nguyen, A.; Fuchs, T.; and Lipson, H. 2015.  Understanding neural networks through deep visualization.
  * Zablocki et al. (2022) Zablocki, É.; Ben-Younes, H.; Pérez, P.; and Cord, M. 2022.  Explainability of deep vision-based autonomous driving systems: Review and challenges.  _IJCV_.
  * Zeiler and Fergus (2014) Zeiler, M. D.; and Fergus, R. 2014.  Visualizing and understanding convolutional networks.  In _ECCV_.
  * Zerhouni et al. (2017) Zerhouni, E.; Lányi, D.; Viana, M.; and Gabrani, M. 2017.  Wide residual networks for mitosis detection.  In _ISBI_. IEEE.
  * Zhang et al. (2018) Zhang, Q.; Cao, R.; Shi, F.; Wu, Y. N.; and Zhu, S.-C. 2018.  Interpreting CNN knowledge via an explanatory graph.  In _AAAI_.
  * Zhang, Wu, and Zhu (2018) Zhang, Q.; Wu, Y. N.; and Zhu, S.-C. 2018.  Interpretable convolutional neural networks.  In _CVPR_.
  * Zhang et al. (2019) Zhang, Q.; Yang, Y.; Ma, H.; and Wu, Y. N. 2019.  Interpreting cnns via decision trees.  In _CVPR_.
  * Zhou et al. (2015) Zhou, B.; Khosla, A.; Lapedriza, À.; Oliva, A.; and Torralba, A. 2015.  Object Detectors Emerge in Deep Scene CNNs.  In _ICLR_.
  * Zhou et al. (2016) Zhou, B.; Khosla, A.; Lapedriza, A.; Oliva, A.; and Torralba, A. 2016.  Learning deep features for discriminative localization.  In _CVPR_.


# Binary black holes mergers from hierarchical triples in open clusters

Dylan Britt 1,2, Ben Johanson 1, Logan Wood 1,3, M. Coleman Miller 1,4 and
Erez Michaely 1
1Department of Astronomy, University of Maryland, College Park, MD 20742, USA
2Kavli Institute for Particle Astrophysics & Cosmology, P.O. Box 2450,
Stanford University, Stanford, CA 94305, USA
3Department of Physics and Astronomy, George Mason University, Fairfax,
Virginia 22030, USA
4Joint Space-Science Institute, University of Maryland, College Park, MD
20742, USA E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

A promising channel for producing binary black hole mergers is the Lidov-Kozai
orbital resonance in hierarchical triple systems. While this mechanism has
been studied in isolation, the distribution of such mergers in time and across
star-forming environments is not well characterized. In this work, we explore
Lidov-Kozai-induced black hole mergers in open clusters, combining semi-
analytic and Monte Carlo methods to calculate merger rates and delay times for
eight different population models. We predict a merger rate density of
$\sim$1–10 Gpc-3 yr-1 for the Lidov-Kozai channel in the local universe, and
all models yield delay-time distributions in which a significant fraction of
binary black hole mergers (e.g., $\sim$20%–50% in our baseline model) occur
during the open cluster phase. Our findings suggest that a substantial
fraction of mergers from hierarchical triples occur within star-forming
regions in spiral galaxies.

###### keywords:

black hole mergers – celestial mechanics – galaxies: clusters: general

††pubyear: 2021††pagerange: Binary black holes mergers from hierarchical
triples in open clusters–Binary black holes mergers from hierarchical triples
in open clusters

##  1 Introduction

The detection of binary black hole (BBH) mergers via gravitational wave (GW)
emission became routine by the O3 observational run of the LIGO and Virgo
collaborations. To date, tens of BBH mergers have been detected, with an
overall merger rate density of $\mathcal{R}_{\rm
BBH}=23.9^{+14.9}_{-8.6}$Gpc-3 yr-1 (Abbott et al., 2020). The identification
of relevant channels which lead to mergers via GW emission is an ongoing
endeavour which spans a number of subfields, including orbital dynamics,
stellar evolution, and dynamics on the scale of galaxies.

Channels for BBH mergers may be grouped into four broad categories. The first,
isolated binary stellar evolution of massive stars (e.g., Tutukov & Yungelson,
1973, 1993; Lipunov et al., 1997; Bethe & Brown, 1998; Portegies Zwart &
Yungelson, 1998; Kalogera, 2000; Mandel & O’Shaughnessy, 2010; Voss & Tauris,
2003; Kalogera et al., 2007; Belczynski et al., 2008; Dominik et al., 2012,
2013, 2015; de Mink & Belczynski, 2015; Belczynski et al., 2016; Eldridge et
al., 2017; Giacobbo et al., 2018; Olejak et al., 2020), proposes that some
massive stellar binaries evolve to short-period binaries prior to either star
forming a BH. One type of such evolution occurs during one or two common
envelope episodes, in which one star swells during the giant phase, imparting
drag on the other and shrinking their mutual orbit. The total orbital energy
loss is directly related to the amount of energy transferred to the envelope
of the giant star. If the energy transfer is too efficient then the binary
merges before the objects turn into BHs; if the transfer is too inefficient
then the binary does not lose enough orbital energy to merge via GW emission.
The result is a short-period stellar binary which then evolves to a BBH and
merges via GW emission within a Hubble time. Studies of this channel predict a
delay-time distribution $\propto t^{-1}$ that starts $10-100\rm{Myr}$ after
star formation. They also predict no measurable eccentricity in the LIGO
detection band (due to circularisation during the binary interaction phase and
the subsequent circularisation from gravitational radiation) and merger rate
densities of $\sim$10-2–103 Gpc-3 yr-1. Another isolated binary formation
scenario is the chemically homogeneous stellar evolution (Marchant et al.,
2016; de Mink & Mandel, 2016; Mandel & de Mink, 2016). In this scenario a
massive binary that is close to contact experiences intense internal mixing
that keeps the stars chemically homogeneous while the cores are burning
hydrogen. The hydrogen in the star is thus nearly exhausted and thus a common
envelope phase is avoided. The predicted BBH merger rate is up to $500~{}{\rm
Gpc}^{-3}~{}{\rm yr}^{-1}$ (de Mink & Mandel, 2016).

A second merger channel is dynamical in nature and proposes that existing BBHs
are induced to merge in dense environments such as galactic centers, AGN
accretion disks, or globular clusters. In these settings, BBHs experience
strong gravitational interactions with individual stars or high-multiplicity
systems, and these interactions tend to harden the target binaries and may
increase their eccentricities Sigurdsson & Phinney (1993); Kulkarni et al.
(1993); Portegies Zwart & McMillan (2000); Madau & Rees (2001); Miller &
Hamilton (2002a); Gültekin et al. (2004, 2006); Miller & Lauburg (2009);
McKernan et al. (2012); Samsing et al. (2014); Rodriguez et al. (2016); Stone
et al. (2017); Rodriguez et al. (2018); Fragione & Kocsis (2018); Banerjee
(2018); Hamers et al. (2018); Leigh et al. (2018); Rodriguez et al. (2021).
Models of these interactions predict merger rate densities of $\sim$2–25 Gpc-3
yr-1.

The third channel concerns mergers of initially wide, isolated systems, either
binaries or triples, in the field of the host galaxy (Michaely & Perets, 2019,
2020; Michaely, 2020). For wide systems, the field of the host galaxy is
considered a collisional environment due to frequent flyby interactions with
field stars. These interactions are capable of exciting the eccentricity (in
the case of binaries) or outer eccentricity (in case of triples), with the
result that mergers occur via increased GW emission (binaries) or three-body
instabilities (triples). Predicted BBH merger rate densities for this channel
are $\sim$1–100 Gpc-3 yr-1.

The fourth merger channel, and the focus of this paper, is secular evolution
in hierarchical triple systems. These systems reside either in the field of
the host galaxy (e.g., Antonini et al., 2016, 2017; Silsbee & Tremaine, 2017)
or in dense environments (Miller & Hamilton, 2002b; Antonini & Perets, 2012;
Antonini et al., 2014; Kimpson et al., 2016; Petrovich & Antonini, 2017;
Samsing & D’Orazio, 2018; Hoang et al., 2018; Fragione et al., 2019; Hamilton
& Rafikov, 2019; Martinez et al., 2020; Wang et al., 2020). In this channel, a
BBH experiences secular effects due to its tertiary companion in the form of
the Lidov-Kozai resonance (Lidov, 1962; Kozai, 1962a; Harrington, 1968; Lidov
& Ziglin, 1976; Innanen et al., 1997; Ford et al., 2000; Blaes et al., 2002);
for a recent review, see Naoz (2016). In this resonance, the eccentricity of
the BBH experiences cyclic changes which boost the GW emission rate of the
inner binary and lead to a merger. Predicted merger rate densities due to this
channel are $\sim$0.5–15 Gpc-3 yr-1.

Distinguishing the various channels for producing BBH mergers is important,
given that each may yield mergers with particular observational signatures and
with different spatial or temporal distributions. BBH mergers in open clusters
have been studied previously via $N$-body simulations (Banerjee, 2018;
Kumamoto et al., 2019; Di Carlo et al., 2019; Di Carlo et al., 2020; González
et al., 2020; Weatherford et al., 2021), which predict merger rate densities
of $\sim$0.3 Gpc-3 yr-1 in these environments. Michaely & Perets (2018) found
that a small fraction, up to fraction of a percent, of mergers are expected to
occur extremely close in time to the formation of the second BH, specifically
within years to decades following the supernova. Open clusters are loosely
bound groups of young stars with stellar number densities $n_{*}\sim 0.1$–$10$
pc-3 and typical velocity dispersion $\sigma\sim 1$–$5$ km s-1 (Moraux, 2016).
We assume that effectively all star formation occurs in these clusters (Lada &
Lada, 2003), which remain bound for lifetimes ranging from $\sim$100$\,\rm
Myr$ for the sparsest examples to a few Gyr for the densest clusters (Moraux,
2016).

In this work, we apply semi-analytic modeling and Monte Carlo simulations to
the hierarchical triple channel, studying a set of models describing different
initial triple system populations. For each model, we calculate the total BBH
merger rate density as well as the cumulative distribution of mergers as a
function of time since star formation; this is known as the delay-time
distribution (DTD). In particular, we calculate the fraction of mergers which
occur while a triple still resides in its birth cluster. The main focus of
this work is to estimate the fraction of mergers in the open cluster phase out
of the total mergers induced by the secular evolution.

We begin by describing our semi-analytic treatment of BBH mergers induced by
the secular Lidov-Kozai resonance in Section 2. In Section 3, we then
establish our numerical approach and the different population models
considered. Section 4 presents the simulation results, including the DTD and
merger rate for each model. Section 5 discusses our model assumptions and
limitations, and in Section 6, we summarise and offer broader context for our
results.

##  2 BBH mergers from hierarchical triples

In the following section, we briefly describe the secular evolution of triple
systems under the Lidov-Kozai resonance. For a more detailed description of
this mechanism, see Naoz (2016).

###  2.1 Newtonian treatment

A hierarchical triple system is composed of an inner binary with masses
denoted $m_{1},m_{2}$ and a distant tertiary of mass $m_{3}$. The inner binary
is characterised by its orbital semimajor axis (SMA) $a_{1}$ and eccentricity
$e_{1}$. The center of mass of the inner binary then hierarchically
constitutes an additional two-body system with the tertiary; this system is
referred to as the outer binary, with SMA $a_{2}$ and eccentricity $e_{2}$.
Each binary defines a unique orbital plane, and the angle between these two
planes is the inclination $I$ associated with the triple system. Within these
planes, the orientations of the inner and outer orbits are given by their
arguments of pericenter $\omega_{1}$ and $\omega_{2}$, respectively. See Fig.
1 for a diagram of a general hierarchical triple system.

Figure 1: Illustration of a hierarchical triple system. The inner binary
consists of two objects with masses $m_{1}$ and $m_{2}$ whose orbit is defined
by a SMA $a_{1}$ and eccentricity $e_{1}$. In this study we set $e_{1}=0$. The
outer binary consists of the tertiary of mass $m_{3}$ and the center of mass
of the inner binary. The orbit of the outer binary is defined by a SMA $a_{2}$
and eccentricity $e_{2}$. The angle between the planes of the inner and outer
binaries is the system inclination $I$.

A three-body system is chaotic when the system masses and separations are
similar. Such a system thus tends to break apart on dynamical timescales.
Hence, on grounds of system stability, most astrophysical triple systems are
hierarchical in scale; i.e., $a_{1}\ll a_{2}$. This hierarchy of spatial
scales sets a corresponding hierarchy of timescales for these systems: the
inner binary orbital period $P_{1}$ is much shorter than the outer binary
orbital period $P_{2}$, and any dynamical evolution of the system occurs on
timescales much longer than both.

When the secular approximation is applied to hierarchical triple systems, one
can show that the orbital energies of each binary are conserved quantities,
and therefore the SMAs $a_{1}$ and $a_{2}$ are constant in time. Long-term
changes to the system do occur, however, due to mutual torque and angular
momentum transfer between the inner and outer binaries. The result of this
secular evolution is simultaneous oscillations of the inner eccentricity
$e_{1}$ and system inclination $I$, such that the total angular momentum of
the triple system is conserved; at higher order, the outer orbit can evolve as
well. Peak eccentricity in the inner binary occurs at the time of minimum
inclination, and vice versa, and these oscillations are referred to as Lidov-
Kozai cycles (Lidov, 1962; Kozai, 1962b).

Following Miller & Hamilton (2002b) and VanLandingham et al. (2016), we define
a conserved quantity derived from the quadrupole-order Hamiltonian for a
hierarchical triple system:

$W_{\rm
N}=-2\epsilon+\epsilon\cos^{2}I+5(1-\epsilon)\sin^{2}\omega_{1}(\cos^{2}I-1)\,,$
(1)

where $\epsilon\equiv 1-e_{1}^{2}$. The minimum value of $\epsilon$, which
corresponds to the maximum value of $e_{1}$, occurs when $\omega_{1}=\pi/2$.
Hence, knowing the initial values $\omega_{1,0}$ and $e_{1,0}$, one can
exploit the conservation of $W$ to calculate the maximum value of the inner
binary eccentricity, denoted $e_{\rm max}$. We note here that the octupole-
order result is different (Harrington, 1968; Ford et al., 2000; Blaes et al.,
2002; Thompson, 2011; Naoz et al., 2013; Michaely & Perets, 2014; Naoz, 2016)
but is beyond the scope of this work.

Innanen et al. (1997) provide a concise and useful relation between the
initial inclination and the maximal eccentricity due to the Lidov-Kozai
resonance in the quadrupole approximation when the tertiary dominates the
system angular momentum:

$e_{\rm max}=\left(1-\frac{5}{3}\cos^{2}I_{0}\right)^{1/2}\,,$ (2)

which implies that for the restricted three-body problem, the inner binary
eccentricity tends to unity if $I_{0}=\pi/2$. The growth of the inner
eccentricity to its maximum value over long timescales is a consequence of
coherent perturbations by the potential of the tertiary, specifically inner
binary precession. If the inner eccentricity is sufficiently high, one might
expect the inner binary’s components to interact and thus to disrupt this
precession. In the following subsection, we consider such an effect in general
relativity (GR), namely GR pericenter precession in the inner binary.

###  2.2 Post-Newtonian treatment

In a triple system whose inner binary evolves to high sufficiently high
eccentricity, GR precession of the inner binary pericenter becomes
nonnegligible. This precession interferes with the coherent perturbations due
to the tertiary and suppresses the Lidov-Kozai resonance. Following Miller &
Hamilton (2002b), we account for this quenching effect of GR precession by
adding to equation (1) the following post-Newtonian term:

$W_{\rm
PN}=\frac{8}{\sqrt{\epsilon}}\frac{M_{1}}{m_{3}}\left(\frac{b_{2}}{a_{1}}\right)^{3}\frac{GM_{1}}{a_{1}c^{2}}\equiv\theta_{{\rm
PN}}\epsilon^{-1/2}\,.$ (3)

Here $M_{1}\equiv m_{1}+m_{2}$ is the total mass of the inner binary,
$b_{2}=a_{2}\left(1-e_{2}^{2}\right)^{1/2}$ is the semi-minor axis of the
outer binary, $G$ is the Newtonian gravitational constant, and $c$ is the
speed of light. Note that we include a term for GR pericenter precession but
continue to treat GW emission as negligible for the purposes of the Lidov-
Kozai resonance; as a result, the sum of equations (1) and (3),

$W=W_{N}+W_{\rm PN}\,,$ (4)

remains a conserved quantity. As before, the maximal eccentricity (minimal
$\epsilon$) is obtained when $\omega=\pi/2$, and the result in this post-
Newtonian treatment becomes

$\epsilon_{{\rm min}}^{1/2}\approx\frac{1}{6}\left(\theta_{{\rm
PN}}+\sqrt{\theta_{{\rm PN}}^{2}+60\cos^{2}I_{0}}\right)\,.$ (5)

This maximal eccentricity can be used to estimate the merger time of the inner
binary due to GW emission. The merger timescale for a binary of eccentricity
$e\approx 1$ is given by Peters (1964) as

$T_{\rm
GW}\approx\frac{768}{425}T_{c}\text{$\left(a_{1}\right)$}\left(1-e^{2}\right)^{7/2}\,,$
(6)

where $T_{c}\equiv a_{1}^{4}/\beta$ is the merger timescale for a circular
binary and $\beta\equiv 64G^{3}m_{1}m_{2}\left(m_{1}+m_{2}\right)/(5c^{5})$.
However, in the case of a triple system whose inner binary oscillates between
its initial eccentricity $e_{1,0}$ and maximal eccentricity $e_{\rm max}$, the
merger timescale due to GW emission is necessarily longer. Randall & Xianyu
(2018) analytically estimate the merger time in this case to be

$T_{\rm merger}=\frac{T_{\rm GW}}{\epsilon_{\rm min}^{1/2}}\,.$ (7)

In this work, we are interested in merger times $T_{\rm merger}<T_{\rm
Hubble}$ for the purpose of calculating the merger rate density and DTD of BBH
mergers originating from hierarchical triples. In the following section, we
describe a method for numerically selecting different triple populations in
order to calculate these statistics.

##  3 Numerical method

Our approach to calculating merger rate densities and DTDs is as follows. In
each of several population models, described in Sections 3.2 and 3.3, we
employ a Monte Carlo simulation to generate $10^{6}$ representative triple
systems. For each triple system in a given model, the model analytically
determines whether an inner BBH merger occurs within $T_{\rm Hubble}$ using
equation (7) and records the value of $T_{\rm merger}$ in order to calculate
the theoretical DTD.

As stated earlier, we assume that star formation occurs entirely within open
clusters (Lada & Lada, 2003) and that as a result, black hole progenitor stars
all form simultaneously; i.e., we do not calculate any detailed dynamical
effects during the main sequence (MS) phase.

###  3.1 Creating a population model

Each population model uses a Monte Carlo approach to generate a set of stable,
hierarchical triple systems whose inner binaries evolve to BBHs. Although
binary stellar evolution processes are beyond the scope of this study, we do
consider basic restrictions imposed on triple systems due to their passage
through the MS phase. Specifically, we exclude triples whose inner binary
components would have interacted as MS stars. We exclude any triple that is
considered dynamically unstable by the criterion of Mardling & Aarseth (2001).
Additionally, we work under the simplifying assumption that the probability
distributions of all system parameters are independent, meaning that an
individual triple system can be generated by drawing each of its parameters
independently.

Triple systems are produced in this model by drawing initial stellar masses
and orbital parameters, then mapping those stellar masses to final BH masses.
To generate the inner binary for a system, we draw the primary mass $m_{1}$
from the Kroupa initial mass function (IMF) (Kroupa, 2001), denoted $f_{\rm
IMF}\left(m\right)$, with a range $[m_{\rm min},m_{\rm max}]$. Because we are
interested in masses of BH progenitors, we concern ourselves only with the
upper end of the range of initial masses. The Kroupa and Salpeter IMFs
(Salpeter, 1955) are similar in the high-mass regime, and therefore we do not
expect that a different choice of IMF would affect the results presented here.
However, the specific choice of IMF is important for the normalisation of the
results; see Section 4.1.

With $m_{1}$ determined, the next parameter drawn is the inner binary SMA,
$a_{1}$. Motivated by observations (Duchêne & Kraus, 2013; Moe & Di Stefano,
2016), we draw the inner SMA from a log-uniform distribution (Öpik’s law) over
a range $[a_{1,\rm min},a_{1,\rm max}]$. The mass of the second inner binary
object is determined by

$m_{2}=m_{1}q_{1}\,,$ (8)

where $q_{1}$ is the inner binary mass ratio, drawn from a power law
distribution $f(q)\propto q^{\gamma}$. For high-mass stars ($M_{*}\gtrsim
16\,M_{\odot}$), this distribution covers the range $[0.1,1]$. The power law
index $\gamma$ is determined by the SMA of the inner binary, with $\gamma=0$
for $a_{1}<100\,\rm AU$ and $\gamma=-1/2$ for $a_{1}>100\,\rm AU$ (Duchêne &
Kraus, 2013).

The remaining parameters of the inner binary orbit are its eccentricity
$e_{1}$ and argument of pericenter $\omega_{1}$. In order to be conservative
with merger time we set the inner binary eccentricity to be zero,
$e_{1}\rightarrow 0$. Any other choise of inner eccentricity distribution
would shorten the merger timescale because of the increase of the maximal
eccentricity reached in the Lidov-Kozai resonance (Lidov & Ziglin, 1976; Naoz,
2016). Finally, $\omega_{1}$ is drawn from a uniform distribution on
$[0,2\pi]$. These five parameters define our inner binary progenitor star
system.

As mentioned previously, we discard any system whose inner binary would have
interacted during the MS phase. To check for such interactions, our method
calculates the radii of the progenitor stars and compares these to the stars’
respective Roche limits. The stellar radius-mass relation is given by
$r_{i}\propto m_{i}^{0.57}$ (Demircan & Kahraman, 1991) and the Roche limit by

$R_{1(2)}=a_{1}\times\frac{0.49q_{1(2)}^{2/3}}{0.6q_{1(2)}^{2/3}+\ln{\left(1+q_{1(2)}^{1/3}\right)}}\,,$
(9)

where $q_{1(2)}=m_{1(2)}/m_{2(1)}$ (Eggleton, 1983). A system is discarded if
$r_{i}>R_{i}$ for either progenitor star, reflecting the likelihood that such
a system would have interacted significantly during the MS phase and might
have failed to produce a BBH.

We now address the parameters characterising the outer binary. The tertiary
mass $m_{3}$ is set by drawing the outer mass ratio $q_{2}\equiv m_{3}/M_{1}$
from a power law distribution $q\propto M_{1}^{\gamma}$ with $\gamma=-2$ (Moe
& Di Stefano, 2016) over a range $[0.1,1]$. The outer eccentricity $e_{2}$ is
drawn from a thermal distribution $f(e)=2e$ and the outer SMA $a_{2}$ from a
log-uniform distribution over a range $[a_{2,\rm min},a_{2,\rm max}]$. The
final parameter needed to specify the triple system is the mutual orbital
inclination $I$; this value is drawn from a distribution function
$f\left(I\right)$ which varies by model and is discussed further in the
following sections.

With the system parameters fully determined, our method next checks that the
triple is indeed dynamically stable. The outer pericenter distance is given by
$R^{\rm out}_{P}=a_{2}(1-e_{2})$, and Mardling & Aarseth (2001) define the
stability threshold

$\kappa=2.8\left[(1+q_{2})\frac{(1+e_{2})}{(1-e_{2})}\right]^{2/5}a_{1}\,,$
(10)

which specifies the smallest outer pericenter value for which the system
remains stable. Accordingly, a system is discarded by our model if $R_{P}^{\rm
out}<\kappa$.

The steps described to this point are sufficient to generate a stable,
hierarchical, stellar triple. The initial masses of the three system
components must now be mapped to the final masses of the BHs or other objects
to which they evolve. When the simulation generates a star of sufficient mass,
it converts it into a BH by establishing two mass regimes. For a star whose
initial mass $m_{i}$ falls in the range $20\,M_{\odot}\leq m_{i}\leq
60\,M_{\odot}$, the resulting BH is assigned a final mass $m_{i}/2$, in
keeping with the approximate relation between progenitor mass and final BH
mass for stars in this range. A star with initial mass $m_{i}>60\,M_{\odot}$
is converted to a BH with a final mass of $30\,M_{\odot}$, reflecting the
significant mass loss experienced by very massive MS stars.

The tertiary is treated differently from the initial binary, as it does not
necessarily evolve to a BH. For $m_{3}\leq 8\,M_{\odot}$, the simulation
checks the MS lifetime for that mass; if it is less than $T_{\rm merger}$,
then $m_{3}$ is converted to a $1$-$M_{\odot}$ white dwarf. In this case, we
ignore any expansion of the outer SMA $a_{2}$, given that the expected mass
loss of the tertiary stellar companion in this case is negligible relative to
the total mass of the triple. To obtain the total merger time, the original MS
lifetime is then added to the merger time for the white dwarf system.

For a tertiary in the range $8\,M_{\odot}<m_{3}<20\,M_{\odot}$, i.e., the mass
range for forming a neutron star (NS), the calculation is stopped and the
system discarded. In this case, it is expected that the triple system will be
disrupted by the natal kick of the NS (Hobbs et al., 2005), precluding any
secular evolution.

The following section introduces the baseline (“standard”) population model,
which adopts the most plausible assumptions for the various triple system
parameter distributions. A set of additional models then extends the standard
model by modifying a single assumption at a time.

###  3.2 Standard model

The baseline model assumes that BHs are formed with no natal kicks, either
because of a failed supernova or massive fallback. The limits on the primary
mass are set to $m_{1,\rm min}=30\,M_{\odot}$ and $m_{1,\rm
max}=100\,M_{\odot}$. While $20$-$M_{\odot}$ O-type stars might produce BHs,
there is considerable speculation regarding which mass ranges will yield a
natal kick when forming compact objects. By raising the minimum mass for BH
formation in our simulation, we impose a conservative buffer which makes it
more likely that natal kicks can be neglected.

This model sets the bounds on the inner binary SMA to $a_{1,\rm min}=0.1\,\rm
AU$ and $a_{1,\rm max}=100\,\rm AU$. In keeping with the focus on hierarchical
triples, the outer binary SMA is assigned a lower bound of $a_{2,\rm
min}=5a_{1}$ and an upper bound of $a_{2,\rm max}=1000\,\rm AU$. This upper
bound is determined by the environment: we do not expect open clusters to
contain ultra-wide systems, as these would be ionized due to the relatively
high stellar density in the cluster. For a more detailed treatment of the open
cluster environment, see Section 5. As mentioned in section 3.1 both the inner
SMA, $a_{1}$, and the outer SMA, $a_{2}$ distributions are equal in log
intervals of $a_{1}$ and $a_{2}$ respectively.

The inclination $I$ of each system is of particular interest when studying the
Lidov-Kozai resonance. Given the dearth of observational constraints on the
inclinations of high-multiplicity systems within open clusters, we make the
reasonable assumption that open clusters and their constituents exhibit a bias
toward aligned angular momenta. For triple systems, such a bias favors
coplanar orbits. To account for this preference, the standard model draws
inclinations from a distribution which increases linearly with $\cos{I}$ in
the range $\cos{I}\in\left[-1,1\right]$; see Fig. 2.

###  3.3 Additional models

To probe the sensitivity of merger rates and the DTD to the assumptions used
in the standard model, we present several additional models. Each isolates and
modifies a single assumption in order to test the robustness of our results.

#### No Natal Kicks

The standard model excludes primary object masses below $m_{1,\rm
min}=30\,M_{\odot}$ due to uncertainty regarding BH natal kicks below this
mass. In this No Natal Kicks model, it is assumed that _all_ BHs are born with
no natal kick, and thus $m_{1,\rm min}$ is lowered to the traditionally
accepted lower limit of $20\,M_{\odot}$ for BH progenitors.

#### Isotropic Distribution

In order to test the sensitivity of our results to the initial distribution of
mutual inclinations, this model implements an isotropic (rather than prograde-
biased) distribution for $I$. Inclinations are drawn from a uniform
distribution of $\cos I\in\left[-1,1\right]$, i.e., from prograde to
retrograde mutual inclinations. See Fig. 2 for the distribution of
inclinations generated by this model.

#### Prograde-Only

This model restricts the mutual inclination $I$ to prograde values by drawing
from a linear distribution of $\cos I\in\left[0,1\right]$. See again Fig. 2
for the initial distribution of inclinations.

#### BH Tertiary

This and the following model concern modifications to the tertiary object in
the triple system. In the standard model, the tertiary star is either massive
enough to become a BH, forming a hierarchical triple BH, or has a mass low
enough to evolve to a white dwarf. Recall that if the tertiary mass falls in
the intermediate regime $8\,M_{\odot}<m_{3}<20\,M_{\odot}$, it is assumed to
form a NS and disrupt the triple via a high natal kick velocity. In this BH
Tertiary population model, only tertiary companions which form black holes are
included, and so only systems with tertiary masses $m_{3}>20\,M_{\odot}$ are
considered.

#### Stellar Tertiary

Complementary to the previous model, here only lower-mass tertiary objects are
allowed. The evolution of these stars is modeled in two phases, as previously
described in Section 3.1. In the first, a star retains its zero-age MS mass
$m_{3}$. In the second phase, the mass of the tertiary star is set to
$1\,M_{\odot}$ to account for mass loss during the giant phases and final
evolution to a white dwarf. As before, these low-mass tertiaries are
restricted to $m_{3}\leq 8\,M_{\odot}$.

#### SMA Boundaries Model a

In all previous models, the inner binary SMA $a_{1}$ is drawn from the range
$\left[0.1\,\rm AU,100\,AU\right]$. This model considers only larger inner
binaries by increasing the lower bound of the inner binary SMA by an order of
magnitude, drawing $a_{1}\in\left[1\,\rm AU,100\,\rm AU\right]$.

#### SMA Boundaries Model b

This model complements the previous model by doubling the upper bound on the
inner binary SMA, drawing $a_{1}\in\left[0.1\,\rm AU,200\,\rm AU\right]$.

Figure 2: Top row: Initial parameter distributions produced by the Monte Carlo simulation for the various population models. Left panel: inner binary SMA; middle panel: outer binary SMA; right panel: system inclination plotted as $\cos I_{0}$. Note that in the distribution of inclinations, only the standard, Isotropic, and Prograde models are shown; all others exhibit no significant differences from the standard model in their distributions of inclinations. Bottom row: the same three parameter distributions shown in the top row, but restricted to the subset of triple systems which merge within a Hubble time in our model. All plots are normalized to unity. Table 1: Parameter distributions and merger results for the standard and additional models. For comparison with the lifetimes of open clusters, the percentages of systems which have merged at $10^{8}$ and $10^{9}$ yr are reported. Model |  $\boldsymbol{m_{1}\,(\rm M_{\odot})}$ |  $\boldsymbol{m_{3}\,(\rm M_{\odot})}$ |  Inclination ($\boldsymbol{f(I)}$) |  $\boldsymbol{a_{1}}$ (AU) |  Local Rate |  Merger Time
---|---|---|---|---|---|---
|  |  |  |  |  (Gpc-3 yr-1) |  $\boldsymbol{\leq$$10^{9}}$ yr ($\boldsymbol{10^{8}}$ yr)
Standard |  30–100 |  $m\leq 8$ |  linear in $\cos I$ |  0.1–100 |  6.2 |  49.9 (18.9) %
|  |  $m\geq 30$ |  |  |  |
No Natal Kicks |  20–100 |  $m\leq 8$ |  linear in $\cos I$ |  0.1–100 |  4.5 |  48.4 (18.1)%
|  |  $m\geq 20$ |  |  |  |
Isotropic Distribution |  30–100 |  $m\leq 8$ |  uniform in $\cos I$ |  0.1–100 |  6.6 |  51.0 (19.5) %
|  |  $m\geq 30$ |  |  |  |
Prograde-Only |  30–100 |  $m\leq 8$ |  linear in $\cos I$, |  0.1–100 |  2.1 |  33.3 (7.9)%
|  |  $m\geq 30$ |  0 $\leq I\leq$ 1 |  |  |
BH Tertiary |  30–100 |  $m\geq 30$ |  linear in $\cos I$ |  0.1–100 |  6.7 |  50.5 (19.3)%
Stellar Tertiary |  30–100 |  $m\leq 8$ |  linear in $\cos I$ |  0.1–100 |  0.9 |  15.8 (0.6)%
SMA Boundaries a |  30–100 |  $m\leq 8$ |  linear in $\cos I$ |  1–100 |  3.7 |  51.5 (20.0)%
|  |  $m\geq 30$ |  |  |  |
SMA Boundaries b |  30–100 |  $m\leq 8$ |  linear in $\cos I$ |  0.1–200 |  6.2 |  33.3 (7.9)%
|  |  $m\geq 30$ |  |  |  |

##  4 Results

###  4.1 Normalisation and rates

We calculate the merger rate density for Lidov-Kozai-assisted BBHs under the
assumption that the Milky Way is the prototypical spiral galaxy with a
population of $N\approx 10^{10}$ stars. The fraction of primary objects in our
triple systems which will form BHs is given by

$f_{\rm
p}=\frac{\int_{30\,M_{\odot}}^{100\,M_{\odot}}m^{-2.3}\,dm}{\int_{0.08\,M_{\odot}}^{100\,M_{\odot}}f_{\rm
IMF}(m)\,dm}\ ,$ (11)

We continue to treat BHs as forming in high-multiplicity systems (Duchêne &
Kraus, 2013) and without natal kicks. Therefore, taking a uniform distribution
of mass ratios $q_{1}\in\left[0.1,1\right]$ for the inner binary, the fraction
of secondary stars forming BHs is $f_{\rm s}\approx 0.4$. Drawing from a mass
ratio distribution $q_{2}\propto M_{1}^{-2}$ for a tertiary at large distances
to the inner binary, the fraction of tertiary objects which remain in triple
systems is $f_{\rm t}\approx 0.25$. Recall that all tertiary masses in the
range $\left[8\,M_{\odot},30\,M_{\odot}\right]$ are rejected, as these are
expected to disrupt the triple system due to large natal kicks during NS
formation (Hobbs et al., 2005). The fraction of the total stellar population
which resides in triple systems is taken to be $f_{\rm triple}\approx 0.1$
(Tokovinin, 2004). Recall that because this work concerns triples within open
clusters, we consider only those triples with an maximum outer binary SMA of
$1000\,\rm AU$ and maximum inner binary SMA of $100\,\rm AU$; see Section 5
for a discussion of this choice of values. The fraction of stars which form
triple systems with inner binary BBHs that merge via the Lidov-Kozai resonance
is then given by

$F_{\rm model}=10^{-5}\left(\frac{f_{\rm p}}{10^{-3}}\right)\left(\frac{f_{\rm
s}}{0.4}\right)\left(\frac{f_{\rm t}}{0.25}\right)\left(\frac{f_{\rm
triple}}{0.1}\right)f_{\rm merger}\,,$ (12)

where $f_{\rm merger}$ is the merger fraction for hierarchical triples
calculated by our numerical model. Recall that a triple system is considered
to have merged if $T_{\rm merger}<T_{\rm Hubble}$. The average merger rate for
a single Milky Way-like galaxy over a Hubble time is therefore

$\Gamma_{\rm MW}=N\times\frac{F_{\rm model}}{T_{\rm Hubble}}\approx
0.53\,\rm{Myr^{-1}}\,.$ (13)

Following Belczynski et al. (2016), the merger rate density in the local
universe is given by

$\mathcal{R}=\rho_{\rm gal}\times N\times\frac{F_{\rm model}}{T_{\rm
Hubble}}\approx 6\left(\frac{F_{\rm
model}}{10^{-5}}\right)\rm{Gpc^{-3}\,yr^{-1}}\,,$ (14)

where $\rho_{\rm gal}\approx 0.0116\,\rm{Mpc^{-3}}$ is the Milky Way-like
galaxy density in the local universe (Belczynski et al., 2016). Depending on
the values of the factors that determine $F_{\rm model}$ (see Equation (12)
above), this rate is plausibly comparable to the observed LIGO rate of
$\sim$20$\,\rm{Gpc^{-3}\,yr^{-1}}$.

###  4.2 Delay-time Distribution

Having recorded $T_{\rm merger}$ for each triple system, we can calculate the
fraction of systems which merge within a given time after star formation. Fig.
3 shows the standard model DTD, i.e., the cumulative merger fraction as a
function of time. We find that approximately half of mergers in the standard
model occur within the lifetime of open clusters, suggesting that a
significant fraction of Lidov-Kozai-induced mergers may occur in these
clusters before their dissolution.

Fig. 4 compares the DTD for the standard model to those for the additional
models. Accounting for white dwarf formation in low-mass tertiaries and
allowing all viable systems to evolve in time, we find that $\sim$20%–50% of
Lidov-Kozai-assisted BBH mergers occur within the lifetime of open clusters.
We find that the DTD is not particularly sensitive to model assumptions, with
the exception of the Stellar Tertiary model, which is skewed toward later
merger times and yields a smaller merger fraction within the lifetime of open
clusters. This difference can be understood as the result of lower-mass
tertiary objects, which have weaker effects on the secular evolution of triple
systems.

Figure 3: Delay-time distribution for simulated mergers in the standard model.
The blue curve shows the distribution function for systems which merged within
a Hubble time. The gray box indicates the fraction of mergers occurring during
the open cluster phase, using $10^{9}\,\rm yr$ as the upper limit of an open
cluster lifetime. For comparison, the orange curve shows the DTD which would
result if all BBHs merged in isolation via GW emission alone. Figure 4:
Comparison of delay-time distributions for simulated mergers across all
population models. Curves show the distribution function for systems which
merged within a Hubble time in all models. The gray box indicates the fraction
of mergers which occurred during the open cluster phase in the standard model.
The distribution shows little sensitivity to initial assumptions, with the
exception of low tertiary masses in the Stellar Tertiary model.

##  5 Discussion

###  5.1 Assumptions

Each of the population models developed in this work rests on a set of
underlying assumptions regarding the parameter distributions of its triple
systems. In what follows, we discuss the justification for and implications of
several key model assumptions.

BH natal kicks. The first and most important of these assumptions is that BHs
are born with little or no natal kick. While it remains unclear whether such
kicks are significant (Nelemans et al., 1999; Willems et al., 2005; Wong et
al., 2012; Repetto et al., 2012; Wong et al., 2014; Mandel, 2016; Repetto et
al., 2017), observational evidence supports BH formation via failed SN or
direct collapse (Fryer et al., 1999; Ertl et al., 2016; Adams et al., 2017).
Both mechanisms imply small natal kicks or none at all, supporting the use of
our simplifying assumption. In future work, however, we aim to test the
importance and sensitivity of this assumption by implementing a more
sophisticated population synthesis method.

Triple formation. Throughout this study, we assume that all star formation
occurs in open clusters or associations (Lada & Lada, 2003). The issue of
triple formation is not explicitly addressed; our standard model effectively
treats each hierarchical triple as primordial. This assumption of primordial
system formation is reflected in the non-isotropic distribution of
inclinations used in the standard model. We explore a deviation from this
assumption by including the Isotropic Distribution model, which draws from a
uniform distribution of $\cos I$ and thus simulates triples formed by
dynamical processes. In our results, neither the merger rate density nor the
DTD depends sensitively on the initial distribution of inclinations.

SMA bounds. In the standard model, the lower bound on the inner binary SMA is
set to $a_{1,\rm min}=0.1$ AU. For an isolated BBH with
$m_{1}=m_{2}=20\,M_{\odot}$ and $a_{1}=0.1$ AU in a circular orbit, equation
(6) gives an inspiral time (via GW emission only) of $T_{\rm GW}\approx
10^{10}$ yr, which is on the order of a Hubble time. Therefore, for smaller
values of $a_{1,\rm min}$, we would not expect our Lidov-Kozai channel to
increase the overall rate of BBH mergers. The upper bound $a_{2,\rm max}$ on
the outer binary SMA is set by environmental constraints, specifically the
lifetime of a wide orbit in a collisional environment. Following Bahcall et
al. (1985) one can calculate the half-life of a wide system of SMA $a_{2}$ in
a collisional environment according to

$t_{1/2}=0.00233\,\frac{v_{\text{enc}}}{Gm_{\rm p}n_{*}a_{2}}\ ,$ (15)

where $v_{\rm enc}$ is the typical encounter velocity at infinity, $m_{\rm p}$
the mass of the perturbing body, and $n_{*}$ the local stellar number density.
For an open cluster, we take $v_{\text{enc}}$ to be a typical velocity
dispersion $\sigma\approx 5$ km sec-1 and assume a stellar number density
$n_{*}\approx 0.5$ pc-3 and a perturber mass $m_{\rm p}=1\,M_{\odot}$. Taking
$10^{9}$ yr to be a typical open cluster lifetime, the outer binary SMA of a
system whose half-life is equal to the lifetime of the cluster is
$a_{2}\approx 1000\,\rm AU$; this serves as the upper limit for the size of
the outer binary.

###  5.2 Mergers in the open cluster phase

As summarized in Fig. 4 and in Table 1 for all models considered, the fraction
of mergers occurring during the lifetime of open clusters is significant. In
the standard model, assuming an open cluster lifetime of $10^{9}$ Myr
($10^{8}$ Myr), we find that $49.9\%$ ($18.9\%$) of BBH mergers induced by the
Lidov-Kozai resonance occur in open clusters. This result implies that at
least this fraction of mergers from the secular triple channel occur in young
environments within star-forming galaxies.

##  6 Conclusions

In this work, we calculate the merger rates and DTD of BBH mergers occuring in
hierarchical triple systems within open clusters via the Lidov-Kozai
resonance. This resonance increases the inner binary eccentricity in cycles,
allowing the binary to dissipate orbital energy and inspiral via GW emission.
Given the sensitive dependence of merger time on orbital eccentricity, BBH
mergers in triple systems experiencing the Lidov-Kozai resonance are expected
to occur on much shorter timescales than those in isolated binaries.
Calculating the DTD for hierarchical triples in open clusters, we find that a
significant fraction of mergers ($18\%$–$50\%$ in our baseline model) occur
before the open cluster has dissolved. This result suggests that many mergers
in hierarchical triples occur in star-forming regions and hence in spiral
galaxies.

## Acknowledgements

E.M. thanks the University of Maryland CTC prize fellowship for supporting
this research. The authors thank Selma de Mink, Chris Belczynski and Ilya
Mandel for their useful comments on this manuscript.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Abbott et al. (2020) Abbott R., et al., 2020, arXiv e-prints, p. arXiv:2010.14527
  * Adams et al. (2017) Adams S. M., Kochanek C. S., Gerke J. R., Stanek K. Z., Dai X., 2017, MNRAS, 468, 4968
  * Antonini & Perets (2012) Antonini F., Perets H. B., 2012, ApJ, 757, 27
  * Antonini et al. (2014) Antonini F., Murray N., Mikkola S., 2014, ApJ, 781, 45
  * Antonini et al. (2016) Antonini F., Chatterjee S., Rodriguez C. L., Morscher M., Pattabiraman B., Kalogera V., Rasio F. A., 2016, ApJ, 816, 65
  * Antonini et al. (2017) Antonini F., Toonen S., Hamers A. S., 2017, ApJ, 841, 77
  * Bahcall et al. (1985) Bahcall J. N., Hut P., Tremaine S., 1985, ApJ, 290, 15
  * Banerjee (2018) Banerjee S., 2018, MNRAS, 473, 909
  * Belczynski et al. (2008) Belczynski K., Kalogera V., Rasio F. A., Taam R. E., Zezas A., Bulik T., Maccarone T. J., Ivanova N., 2008, ApJS, 174, 223
  * Belczynski et al. (2016) Belczynski K., Repetto S., Holz D. E., O’Shaughnessy R., Bulik T., Berti E., Fryer C., Dominik M., 2016, ApJ, 819, 108
  * Bethe & Brown (1998) Bethe H. A., Brown G. E., 1998, ApJ, 506, 780
  * Blaes et al. (2002) Blaes O., Lee M. H., Socrates A., 2002, ApJ, 578, 775
  * Demircan & Kahraman (1991) Demircan O., Kahraman G., 1991, Ap&SS, 181, 313
  * Di Carlo et al. (2019) Di Carlo U. N., Giacobbo N., Mapelli M., Pasquato M., Spera M., Wang L., Haardt F., 2019, MNRAS, 487, 2947
  * Di Carlo et al. (2020) Di Carlo U. N., et al., 2020, MNRAS, 498, 495
  * Dominik et al. (2012) Dominik M., Belczynski K., Fryer C., Holz D. E., Berti E., Bulik T., Mandel I., O’Shaughnessy R., 2012, ApJ, 759, 52
  * Dominik et al. (2013) Dominik M., Belczynski K., Fryer C., Holz D. E., Berti E., Bulik T., Mandel I., O’Shaughnessy R., 2013, ApJ, 779, 72
  * Dominik et al. (2015) Dominik M., et al., 2015, ApJ, 806, 263
  * Duchêne & Kraus (2013) Duchêne G., Kraus A., 2013, ARA&A, 51, 269
  * Eggleton (1983) Eggleton P. P., 1983, ApJ, 268, 368
  * Eldridge et al. (2017) Eldridge J. J., Stanway E. R., Xiao L., McClelland L. A. S., Taylor G., Ng M., Greis S. M. L., Bray J. C., 2017, Publ. Astron. Soc. Australia, 34, e058
  * Ertl et al. (2016) Ertl T., Janka H. T., Woosley S. E., Sukhbold T., Ugliano M., 2016, ApJ, 818, 124
  * Ford et al. (2000) Ford E. B., Kozinsky B., Rasio F. A., 2000, ApJ, 535, 385
  * Fragione & Kocsis (2018) Fragione G., Kocsis B., 2018, Phys. Rev. Lett., 121, 161103
  * Fragione et al. (2019) Fragione G., Grishin E., Leigh N. W. C., Perets H. B., Perna R., 2019, MNRAS, 488, 47
  * Fryer et al. (1999) Fryer C. L., Woosley S. E., Hartmann D. H., 1999, ApJ, 526, 152
  * Giacobbo et al. (2018) Giacobbo N., Mapelli M., Spera M., 2018, MNRAS, 474, 2959
  * González et al. (2020) González E., Kremer K., Chatterjee S., Fragione G., Rodriguez C. L., Weatherford N. C., Ye C. S., Rasio F. A., 2020, arXiv e-prints, p. arXiv:2012.10497
  * Gültekin et al. (2004) Gültekin K., Miller M. C., Hamilton D. P., 2004, ApJ, 616, 221
  * Gültekin et al. (2006) Gültekin K., Miller M. C., Hamilton D. P., 2006, ApJ, 640, 156
  * Hamers et al. (2018) Hamers A. S., Bar-Or B., Petrovich C., Antonini F., 2018, ApJ, 865, 2
  * Hamilton & Rafikov (2019) Hamilton C., Rafikov R. R., 2019, MNRAS, 488, 5512
  * Harrington (1968) Harrington R. S., 1968, AJ, 73, 190
  * Hoang et al. (2018) Hoang B.-M., Naoz S., Kocsis B., Rasio F. A., Dosopoulou F., 2018, ApJ, 856, 140
  * Hobbs et al. (2005) Hobbs G., Lorimer D. R., Lyne A. G., Kramer M., 2005, MNRAS, 360, 974
  * Innanen et al. (1997) Innanen K. A., Zheng J. Q., Mikkola S., Valtonen M. J., 1997, AJ, 113, 1915
  * Kalogera (2000) Kalogera V., 2000, ApJ, 541, 319
  * Kalogera et al. (2007) Kalogera V., Belczynski K., Kim C., O’Shaughnessy R., Willems B., 2007, Phys. Rep., 442, 75
  * Kimpson et al. (2016) Kimpson T. O., Spera M., Mapelli M., Ziosi B. M., 2016, MNRAS, 463, 2443
  * Kozai (1962a) Kozai Y., 1962a, AJ, 67, 591
  * Kozai (1962b) Kozai Y., 1962b, AJ, 67, 591
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Kulkarni et al. (1993) Kulkarni S. R., Hut P., McMillan S., 1993, Nature, 364, 421
  * Kumamoto et al. (2019) Kumamoto J., Fujii M. S., Tanikawa A., 2019, MNRAS, 486, 3942
  * Lada & Lada (2003) Lada C. J., Lada E. A., 2003, ARA&A, 41, 57
  * Leigh et al. (2018) Leigh N. W. C., et al., 2018, MNRAS, 474, 5672
  * Lidov (1962) Lidov M. L., 1962, Planet. Space Sci., 9, 719
  * Lidov & Ziglin (1976) Lidov M. L., Ziglin S. L., 1976, Celestial Mechanics, 13, 471
  * Lipunov et al. (1997) Lipunov V. M., Postnov K. A., Prokhorov M. E., 1997, MNRAS, 288, 245
  * Madau & Rees (2001) Madau P., Rees M. J., 2001, ApJ, 551, L27
  * Mandel (2016) Mandel I., 2016, MNRAS, 456, 578
  * Mandel & O’Shaughnessy (2010) Mandel I., O’Shaughnessy R., 2010, Classical and Quantum Gravity, 27, 114007
  * Mandel & de Mink (2016) Mandel I., de Mink S. E., 2016, MNRAS, 458, 2634
  * Marchant et al. (2016) Marchant P., Langer N., Podsiadlowski P., Tauris T. M., Moriya T. J., 2016, A&A, 588, A50
  * Mardling & Aarseth (2001) Mardling R. A., Aarseth S. J., 2001, MNRAS, 321, 398
  * Martinez et al. (2020) Martinez M. A. S., et al., 2020, ApJ, 903, 67
  * McKernan et al. (2012) McKernan B., Ford K. E. S., Lyra W., Perets H. B., 2012, MNRAS, 425, 460
  * Michaely (2020) Michaely E., 2020, MNRAS,
  * Michaely & Perets (2014) Michaely E., Perets H. B., 2014, ApJ, 794, 122
  * Michaely & Perets (2018) Michaely E., Perets H. B., 2018, ApJ, 855, L12
  * Michaely & Perets (2019) Michaely E., Perets H. B., 2019, ApJ, 887, L36
  * Michaely & Perets (2020) Michaely E., Perets H. B., 2020, Monthly Notices of the Royal Astronomical Society
  * Miller & Hamilton (2002a) Miller M. C., Hamilton D. P., 2002a, MNRAS, 330, 232
  * Miller & Hamilton (2002b) Miller M. C., Hamilton D. P., 2002b, ApJ, 576, 894
  * Miller & Lauburg (2009) Miller M. C., Lauburg V. M., 2009, ApJ, 692, 917
  * Moe & Di Stefano (2016) Moe M., Di Stefano R., 2016, preprint,  (arXiv:1606.05347)
  * Moraux (2016) Moraux E., 2016, in EAS Publications Series. pp 73–114 (arXiv:1607.00027), doi:10.1051/eas/1680004
  * Naoz (2016) Naoz S., 2016, ARA&A, 54, 441
  * Naoz et al. (2013) Naoz S., Farr W. M., Lithwick Y., Rasio F. A., Teyssandier J., 2013, MNRAS, 431, 2155
  * Nelemans et al. (1999) Nelemans G., Tauris T. M., van den Heuvel E. P. J., 1999, A&A, 352, L87
  * Olejak et al. (2020) Olejak A., Fishbach M., Belczynski K., Holz D. E., Lasota J. P., Miller M. C., Bulik T., 2020, ApJ, 901, L39
  * Peters (1964) Peters P. C., 1964, Physical Review, 136, 1224
  * Petrovich & Antonini (2017) Petrovich C., Antonini F., 2017, ApJ, 846, 146
  * Portegies Zwart & McMillan (2000) Portegies Zwart S. F., McMillan S. L. W., 2000, ApJ, 528, L17
  * Portegies Zwart & Yungelson (1998) Portegies Zwart S. F., Yungelson L. R., 1998, A&A, 332, 173
  * Randall & Xianyu (2018) Randall L., Xianyu Z.-Z., 2018, ApJ, 864, 134
  * Repetto et al. (2012) Repetto S., Davies M. B., Sigurdsson S., 2012, MNRAS, 425, 2799
  * Repetto et al. (2017) Repetto S., Igoshev A. P., Nelemans G., 2017, MNRAS, 467, 298
  * Rodriguez et al. (2016) Rodriguez C. L., Chatterjee S., Rasio F. A., 2016, Phys. Rev. D, 93, 084029
  * Rodriguez et al. (2018) Rodriguez C. L., Amaro-Seoane P., Chatterjee S., Kremer K., Rasio F. A., Samsing J., Ye C. S., Zevin M., 2018, Phys. Rev. D, 98, 123005
  * Rodriguez et al. (2021) Rodriguez C. L., Kremer K., Chatterjee S., Fragione G., Loeb A., Rasio F. A., Weatherford N. C., Ye C. S., 2021, Research Notes of the American Astronomical Society, 5, 19
  * Salpeter (1955) Salpeter E. E., 1955, ApJ, 121, 161
  * Samsing & D’Orazio (2018) Samsing J., D’Orazio D. J., 2018, MNRAS, 481, 5445
  * Samsing et al. (2014) Samsing J., MacLeod M., Ramirez-Ruiz E., 2014, ApJ, 784, 71
  * Sigurdsson & Phinney (1993) Sigurdsson S., Phinney E. S., 1993, ApJ, 415, 631
  * Silsbee & Tremaine (2017) Silsbee K., Tremaine S., 2017, ApJ, 836, 39
  * Stone et al. (2017) Stone N. C., Metzger B. D., Haiman Z., 2017, MNRAS, 464, 946
  * Thompson (2011) Thompson T. A., 2011, ApJ, 741, 82
  * Tokovinin (2004) Tokovinin A., 2004, in Allen C., Scarfe C., eds, Rev. Mex. AA Conf. Ser. 27, Environment and Evolution of Double and Multiple Stars. p. 7
  * Tutukov & Yungelson (1973) Tutukov A., Yungelson L., 1973, Nauchnye Informatsii, 27, 70
  * Tutukov & Yungelson (1993) Tutukov A. V., Yungelson L. R., 1993, MNRAS, 260, 675
  * VanLandingham et al. (2016) VanLandingham J. H., Miller M. C., Hamilton D. P., Richardson D. C., 2016, ApJ, 828, 77
  * Voss & Tauris (2003) Voss R., Tauris T. M., 2003, MNRAS, 342, 1169
  * Wang et al. (2020) Wang H., Stephan A. P., Naoz S., Hoang B.-M., Breivik K., 2020, arXiv e-prints, p. arXiv:2010.15841
  * Weatherford et al. (2021) Weatherford N. C., Fragione G., Kremer K., Chatterjee S., Ye C. S., Rodriguez C. L., Rasio F. A., 2021, ApJ, 907, L25
  * Willems et al. (2005) Willems B., Henninger M., Levin T., Ivanova N., Kalogera V., McGhee K., Timmes F. X., Fryer C. L., 2005, ApJ, 625, 324
  * Wong et al. (2012) Wong T.-W., Valsecchi F., Fragos T., Kalogera V., 2012, ApJ, 747, 111
  * Wong et al. (2014) Wong T.-W., Valsecchi F., Ansari A., Fragos T., Glebbeek E., Kalogera V., McClintock J., 2014, ApJ, 790, 119
  * de Mink & Belczynski (2015) de Mink S. E., Belczynski K., 2015, ApJ, 814, 58
  * de Mink & Mandel (2016) de Mink S. E., Mandel I., 2016, MNRAS, 460, 3545