# Fast, Accurate Barcode Detection in Ultra High-Resolution Images

###### Abstract

Object detection in Ultra High-Resolution (UHR) images has long been a
challenging problem in computer vision due to the varying scales of the
targeted objects. When it comes to barcode detection, resizing UHR input
images to smaller sizes often leads to the loss of pertinent information,
while processing them directly is highly inefficient and computationally
expensive. In this paper, we propose using semantic segmentation to achieve a
fast and accurate detection of barcodes of various scales in UHR images. Our
pipeline involves a modified Region Proposal Network (RPN) on images of size
greater than 10k$\times$10k and a newly proposed Y-Net segmentation network,
followed by a post-processing workflow for fitting a bounding box around each
segmented barcode mask. The end-to-end system has a latency of 16
milliseconds, which is $2.5\times$ faster than YOLOv4 and $5.9\times$ faster
than Mask R-CNN. In terms of accuracy, our method outperforms YOLOv4 and Mask
R-CNN by a $mAP$ of 5.5% and 47.1% respectively, on a synthetic dataset. We
have made available the generated synthetic barcode dataset and its code at
http://www.github.com/viplabB/SBD/.

Index Terms—  Barcode detection with deep neural networks, barcode
segmentation, Ultra High-Resolution images.

##  1 Introduction

Barcodes are digital signs often made of adjacent and alternating black and
white smaller rectangles that have become an intrinsic part of human society.
In administration, for example, they are used to encode, save, and retrieve
various users’ information. At grocery stores, they are used to track sales
and inventories. More interestingly in e-commerce, they are used to track and
speed up processing time in warehouses and fulfillment centers.

In classical signal processing, filters used for detection are image-specific
since input images are not all necessarily acquired with the same
illumination, brightness, angle, or camera. Consequently, adaptive image
processing algorithms are required, which can impact detection accuracy [1].
In addition, because classical signal processing methods often run on Central
Processing Units, they tend to be much slower compared with deep learning
implementations that are easily optimized on Graphics Processing Units (GPUs).

Fig. 1: Proposed Approach, the modified RPN is followed by Y-Net and the
bounding box extractor.

Over the years, a number of methods have been proposed to detect barcodes
using classical signal processing [1, 2, 3, 4, 5], but nearly all of them take
too long to process Ultra High-Resolution (UHR) images. More specifically, [5]
used parallel segment detectors which improved on their previous work [6] of
finding imaginary perpendicular lines in Hough space with maximal stable
extremal regions to detect barcodes. Katona et al. [3] used morphological
manipulation for barcode detection, but this method did not generalize well as
different barcode types have varying detection performances. Similarly, [7]
proposed using x and y derivative differences, but varying input images yields
different outputs, and using such operation on UHR images often become highly
inefficient.

With neural networks, though there has been much improvement in barcode
detection tasks, few of them have addressed the fast and accurate detection
problem in UHR images. Zamberletti et al. [8] paved the way for using neural
networks to detect barcodes by investigating Hough spaces. This was followed
by [9] which adapted the You Only Looked Once (YOLO) detector to find barcodes
in Low Resolution (LR) images, but the YOLO algorithm is known to perform
poorly with long shaped objects such as code 39 barcodes. Instance
segmentation methods such as Mask R-CNN [10] perform better on images of size
$1024\times 1024$ pixels but on smaller size images, the outputted Region of
Interests (RoI) do not align well with long, 1D-barcode structures. This is
because it typically predicts masks on $28\times 28$ pixels irrespective of
object size, and thereby generates ”wiggly” artifacts on some barcode
predictions, losing spatial resolution. In the same way, dedicated object
detection pipelines, such as YOLOv4 [11], though they perform well on lower
Intersection over Union (IoU) thresholds, suffer accuracy at higher IoU
thresholds. Among those using segmentation on LR images as a means for
detection, [12] also tends to not perform well at higher IoU thresholds.

In this paper, we propose a pipeline for detecting barcodes using deep neural
networks, shown in Fig. 1, which consists of two stages trained separately.
When compared with classical signal processing methods, neural networks not
only provide a faster inference time, but also yield higher accuracy because
they learn meaningful filters for optimal feature extraction. As seen in Fig.
1, in the first stage, we expand on the Region Proposal Network (RPN)
introduced in Faster R-CNN [13] to extract high definition regions of
potential locations where barcodes might be. This stage allows us to
significantly reduce inference computation time that would have been required
otherwise in the second stage. In the second stage, we introduce Y-Net, a
semantic segmentation network that detects all instances of barcodes in a
given outputted RoI image ($400\times 400$). We then apply morphological
operations on the predicted masks to separate and extract the corresponding
bounding boxes as shown in Fig. 2.

One of the limitations of existing work on barcode detection is the
insufficient number of training examples. ArTe-Lab 1D Medium Barcode Dataset
[8] and the WWU Muenster Barcode Database [14] are two examples of existing
available datasets. They contain 365 and 595 images respectively, with ground
truth masks at a resolution of $640\times 480$. Most of the samples in the
ArTe-Lab dataset have only one EAN13 barcode per sample image, and few of them
in the Muenster database have more than one barcode instance on a given image.
To address this dataset availability problem, we have released 100,000 UHR and
100,000 LR synthetic barcode datasets along with their corresponding bounding
boxes ground truths, and their ground truth masks to facilitate further
studies. The outline of this paper is as follows: in Section 2, we describe
details of our approach; in Section 3, we summarize our experimental results
and in Section 4, we conclude and expend on our future work.

##  2 Proposed Approach

As seen in Fig. 1, our proposed method consists of three stages: the modified
Region Proposal Network stage, our Y-Net 111Our Y-Net architecture resembles
the English alphabet letter “Y” and differs from [15] which used a pre-trained
encoder network that is augmented with an untrained mirrored network and a
decoder network. segmentation network stage, and the bounding box extraction
stage.

###  2.1 Modified Region Proposal Network

Region proposals have been influential in computer vision and more so when it
comes to object detection in UHR images. It is common in UHR images that
barcodes are clustered in a small region of the image. To filter out most of
the non-barcode backgrounds, we modified the RPN introduced in Faster R-CNN
[13] to propose regions of barcodes for our next stages. By first transforming
the UHR input image to an LR input image of size $256\times 256$, the RPN was
trained to identify blobs in LR images. Once a bounding box is placed around
the identified blobs, the resulting proposed bounding box is remapped to the
input UHR image by a perspective transformation, and the resulting regions are
cropped out. The LR input to the RPN is chosen to be of size $256\times 256$
as a lower resolution results in the loss of pertinent information. Non-Max
Suppression (NMS) is used on the predictions to select the most probable
regions.

Fig. 2:  Sample outputs of our pipeline; yellow - segmented barcode pixels;
purple - segmented background pixels; boxes - bounding box extracted; (a)
synthetic barcode image; (b) real barcode image; (c) prediction results on
(a); (d) prediction results on (b).

###  2.2 Y-Net Segmentation Network

Fig. 3:  Y-Net Architecture.

As depicted in Fig. 3, Y-Net is made out of 3 main modules distributed in 2
branches: a Regular Convolutional Module shown in blue which constitutes the
left branch, and a Pyramid Pooling Module shown in brown, along with a Dilated
Convolution Module shown in orange which after concatenation and convolution
constitute the right branch.

The Regular Convolution Module takes in $400\times 400$ output images of the
RPN and consists of convolutional and pooling layers. It starts with 64 -
channel $3\times 3$ kernels and doubles the number at each layer. We alternate
between convolution and max-pooling until we reach a feature map size of
$25\times 25$ pixels. This module allows the model to learn general pixel-wise
information anywhere in the input image.

The Dilated Convolution Module takes advantage of the fact that barcodes have
alternating black and white rectangles to learn sparse features in their
structure. The motivation for this module comes from the fact that dilated
convolution operators play a significant role in the ”algorithme a trous” for
biorthogonal wavelet decomposition [16]. Therefore, the discontinuities in
alternating patterns and sharp edges in barcodes are more accurately learned
by such filters. In addition, they leverage a multiresolution and multiscale
decomposition as they allow the kernels to widen their receptive fields with
dilation rates from 1 up to 16. Here too a $400\times 400$ input image is used
and we maintain 32 – channel $3\times 3$ kernels throughout the module while
the dimensions of the layers are gradually reduced using a stride of 2 until a
feature map of $25\times 25$ pixels is obtained.

The Pyramid Pooling Module allows the model to learn global information about
potential locations of the barcodes at different scales and its layers are
concatenated with the layers on the dilated convolution module in order to
preserve the features extracted from both modules.

The resulting feature maps from the right branch are then added to the output
of the Regular Convolution Module, which allows for the correction of features
that would have been missed by either branch. In other words, the output of
each branch constitutes a residual correction for the other thereby refining
the result at each node as shown in white. The nodes are then up-sampled and
concatenated with transposed convolution feature maps shown in red and yellow
of the corresponding dimension. Throughout the network, we use ReLU as a non-
linearity after each layer and add $L_{2}$ regularization to account for
possible over-fitting scenarios that could have occurred during training. On
all datasets, we use 80% for the training set, 10% for the validation set, and
the remaining 10 % for the testing set. We use one NVIDIA Tesla V100 GPU for
the training process. Since this is a segmentation network and we are
interested in classifying background and barcodes, we use binary cross-entropy
as loss function.

###  2.3 Bounding Box Extraction

Since some images contain barcodes that are really close to each other, their
Y-Net outputs reflect the same configuration which makes the extraction of
individual barcode bounding boxes complex as shown in Fig. 4(a). To separate
them effectively, we perform an erosion, contour extraction, and bounding box
expansion with a pixel correction margin. As shown in Fig. 4(b), the erosion
stage allows the algorithm to widen gaps between segmented barcodes that may
be separated by 1 or more pixels. The resulting mask is then used to infer
individual barcode bounding boxes in the contour extraction stage in Fig. 4(c)
through border following. A pixel correction margin is used to recover the
original bounding boxes’ dimensions during the expansion stage as shown in
Fig. 4(d). This post-processing stage of our pipeline has an average
processing time of 1.5 milliseconds (ms) because it is made of a set of Python
matrix operations to efficiently extract bounding boxes from predicted masks.

|  $mAP$ (all) |  $AP_{50}$ (all) |  $AP_{75}$ (all) |  $mAP$ (small) |  $mAP$ (medium) |  $AR_{50}$ (all) |  $AR_{70}$ (all) |  $AR_{80}$ (all) |  $AR_{90}$ (all) |  Latency (ms) |  Resolution (px)
---|---|---|---|---|---|---|---|---|---|---|---
Mask R-CNN [10] |  .466 |  .985 |  .317 |  .340 |  .489 |  .990 |  .740 |  .279 |  .023 |  94.8 |  $448\times 448$
YOLOv4 [11] |  .882 |  .990 |  .989 |  .815 |  .897 |  1. |  1. |  .995 |  .873 |  40.5 |  $320\times 320$
Ours |  .937 |  .990 |  .990 |  .903 |  .945 |  1. |  1. |  1. |  .972 |  16.0 |  $400\times 400$
|  |  |  |  |  |  |  |  |  |  |
Table 1: Average Precision for Max Detection of 100 and Average Recall for Max Detection of 10 computed using MS COCO API. | Muenster Dataset | ArTe Lab Dataset
---|---|---
| DR | Precision | Recall | mIoU | DR | Precision | Recall | mIoU
Creusot et al. [5] | .982 | - | - | - | .989 | - | - | -
Hansen et al. [9] | .991 | - | - | .873 | .926 | - | - | .816
Namane et al. [1] | .966 | - | - | .882 | .930 | - | - | .860
Zharkov et al. [12] | .980 | .777 | .990 | .842 | .989 | .814 | .995 | .819
ours | 1. | .984 | 1. | .921 | 1. | .974 | 1. | .934
|  |  |  |  |  |  |  |
Table 2: Mean IoU (mIoU), Precison and Recall and Detection Rate (DR) at IoU threshold of 0.5 (Muenster and ArTe-Lab Dataset). | Px Acc | Px mIoU | Px Prec | Px Rec
---|---|---|---|---
Mask R-CNN [10] | .993 | .990 | .989 | .890
Ours | 1. | 1. | .999 | .999
|  |  |  |
Table 3: Pixel-wise Metrics

##  3 Datasets and Results

For the synthetic dataset, we use treepoem
222https://github.com/adamchainz/treepoem and random-word
333https://github.com/vaibhavsingh97/random-word to generate UHR and LR
barcode images. We use Code 39, Code 93, Code 128, UPC, EAN, PD417, ITF, Data
Matrix, AZTEC, and QR among others. We model the number of barcodes in a given
image using a Poisson process and a combination of perspective transforms is
used to make the barcodes vary in shape and position from one image to the
other. We have also added random black blobs at random locations on the
original UHR and LR canvases. The real UHR barcodes dataset obtained from
Amazon.com, Inc is made of 3.8 million UHR images of resolution up to
$30k\times 30k$ grayscale images and could not be released due to
confidentiality reasons. Additionally, the Muenster and Artelab datasets are
used with some data augmentation schemes for more samples.

For the RPN, we accumulated the number of bounding boxes inside the proposed
regions and divided it by the total number of ground truth bounding boxes. Our
implementation yields an accuracy of 98.03% on the synthetic dataset at 10 ms
per image and 96.8% on the real dataset at 13 ms per image while the baseline
[13] yields the same accuracies and an average latency over 2.5 seconds (s)
per image for both datasets.

For Y-Net, we use the Microsoft (MS) COCO API, and Pixel-wise metrics to
evaluate against [10, 11]. By default, the MS COCO API configuration evaluates
on small, medium and large areas objects but in our application, the largest
detected barcode area is medium. Since Y-Net is a segmentation network and
does not output confidence scores for each segmented barcode, we propose using
pseudo scores, the ratio of the total number of nonzero pixels in a predicted
mask to the total number of nonzero pixels in the corresponding ground truth
mask at the location of a given object.

Table 1 shows $mAP$ and $mAR$ values of the models on the synthetic dataset.
As seen, our pipeline outperforms [10], and [11] by a $mAP$ of 47.1% and 5.5%
and $AP_{75}$ of 67.3% and 0.1% respectively. Also shown in Table 1, is a
$mAR_{90}$ improvement of 94.9% and 9.9% on [10] and, [11] respectively which
highlights that Y-Net continues to yield better $mAR$ results even at higher
IoU thresholds. Both our approach and [11] achieve an $AR_{50}$ of 100% and
outperform [10] by 1%. For small area barcodes, Y-Net outperforms [10] and
[11] by a $mAP$ of 56.3% and 8.8% and for medium area barcodes, Y-Net displays
a $mAP$ increase of 45.6% and 4.8% on [10] and [11] respectively. In addition,
Table 3 reveals that Y-Net a has much better semantic segmentation performance
than [10]. Table 1 displays that Y-Net performs at least $2.5\times$ faster
than the fastest of models [10] and, [11] on LR images.

Similarly, we have used the Detection Rate (DR), mIoU, Precision, and Recall,
as described in [1, 5, 9, 12] on the Arte-Lab and Muenster datasets and as can
be seen in Table 2, our method outperforms previous works on all of the
mentioned metrics. This indicates that our bounding box extraction algorithm
is working as expected to detect accurate bounding boxes. However, while it is
successful in separating barcodes that are relatively close to each other, it
has limitations when barcodes are overlapping as shown in Fig. 4(e). For those
occlusion scenarios, the algorithm tends to group the overlapping barcodes
into one bounding box instead of separate bounding boxes as shown in Fig.
4(f).

Fig. 4:  (a) Y-Net output; (b) Y-Net output after erosion; (c) extracted
bounding boxes –red, ground truth bounding boxes –green on eroded output; (d)
final bounding boxes after pixel correction margin on Y-Net output; (e) Y-Net
output of occluded barcodes scenarios; (f) final extracted bounding boxes are
grouped after pixel correction margin due to overlaping barcodes in input
image.

##  4 Conclusion

In this paper, we showed that barcodes can be efficiently, accurately, and
speedily detected using Y-Net on UHR images. With pseudo scores as confidence
scores, our approach outperforms existing detection pipelines with a much
better latency. In future work, we aim to extend this method to the multi-
class detection task for small objects in UHR images and videos in a weakly
supervised fashion.

## References

  * [1] A. Namane and M. Arezki,  “Fast real time 1d barcode detection from webcam images using the bars detection method,”  in Proceedings of the World Congress on Engineering (WCE), 2017, vol. 1.
  * [2] L. Hock, H. Hanaizumi, and E. Ohbuchi,  “Barcode readers using the camera device in mobile phones,”  in 2013 International Conference on Cyberworlds. 2004, pp. 260–265, IEEE Computer Society.
  * [3] M. Katona and L. G Nyúl,  “Efficient 1d and 2d barcode detection using mathematical morphology,”  in Mathematical Morphology and Its Applications to Signal and Image Processing. 2013, pp. 464–475, Springer Berlin Heidelberg.
  * [4] G. Sörös and C. Flörkemeier,  “Blur-resistant joint 1d and 2d barcode localization for smartphones,”  in Proceedings of the 12th International Conference on Mobile and Ubiquitous Multimedia. 2013, MUM ’13, Association for Computing Machinery.
  * [5] C. Creusot and A. Munawar,  “Low-computation egocentric barcode detector for the blind,”  in 2016 IEEE International Conference on Image Processing (ICIP), 2016, pp. 2856–2860.
  * [6] C. Creusot and A. Munawar,  “Real-time barcode detection in the wild.,”  IEEE Winter Conference on Applications of Computer Vision, p. 239–245, 2015.
  * [7] O. Gallo and R. Manduchi,  “Reading 1d barcodes with mobile phones using deformable templates.,”  IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 33, pp. 1834– 1843, 2011.
  * [8] A. Zamberletti, I. Gallo, M. Carullo, and E. Binaghi,  “Neural image restoration for decoding 1-d barcodes using common camera phones.,”  01 2010, vol. 1, pp. 5–11.
  * [9] D. K. Hansen, K. Nasrollahi, C. B. Rasmussen, , and T. B. Moeslund,  “Real-time barcode detection and classification using deep learning.,”  IJCCI, vol. 1, pp. 321–327, 2017.
  * [10] P. Dollár K. He, G. Gkioxari and R. Girshick,  “Mask R-CNN,”  in 2017 IEEE International Conference on Computer Vision (ICCV), 2017, pp. 2980–2988.
  * [11] A. Bochkovskiy, C. Wang, and H. M. Liao,  “YOLOv4: Optimal speed and accuracy of object detection,” 2020.
  * [12] A. Zharkov and I. Zagaynov,  “Universal barcode detector via semantic segmentation,”  2019 International Conference on Document Analysis and Recognition (ICDAR), pp. 837–843, 2019.
  * [13] S. Ren, K. He, R. Girshick, and J. Sun,  “Faster R-CNN: Towards real-time object detection with region proposal networks,”  IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 39, no. 6, pp. 1137–1149, 2017.
  * [14] X. Jiang S. Wachenfeld, S. Terlunen,  Robust 1-D Barcode Recognition on Camera Phones and Mobile Product Information Display, pp. 53–69,  Springer Berlin Heidelberg, Berlin, Heidelberg, 2010.
  * [15] Ahmed Mohammed, Sule Yildirim, Ivar Farup, Marius Pedersen, and Øistein Hovde,  “Y-net: A deep convolutional neural network for polyp detection,” 2018\.
  * [16] M. Holschneider, R. Kronland-Martinet, J. Morlet, and Ph. Tchamitchian,  “A real-time algorithm for signal analysis with the help of the wavelet transform,”  In Wavelets: Time-Frequency Methods and Phase Space. Proceedings of the International Conference, 1987.


11institutetext: Institute of Astronomy, KU Leuven, Celestijnenlaan 200D,
B-3001 Leuven, Belgium
11email<EMAIL_ADDRESS>22institutetext: University of Oxford, St
Edmund Hall, Oxford OX1 4AR, UK 33institutetext: School of Physics and
Astronomy, Monash University, Clayton, VIC 3800, Australia 44institutetext:
The ARC Centre of Excellence for Gravitational Wave Discovery – OzGrav,
Australia

# An upper limit on the spins of merging binary black holes formed through
binary evolution

Pablo Marchant 11 Philipp Podsiadlowski 22 Ilya Mandel 3344

###### Abstract

Context. As ground-based gravitational wave detectors progressively increase
their sensitivity, observations of black hole mergers will provide us with the
joint distribution of their masses and spins. This will be a critical
benchmark to validate different formation scenarios.

Aims. Merging binary black holes formed through the evolution of isolated
binary systems require both components to be stripped of their hydrogen
envelopes before core-collapse. The rotation rates of such stripped stars are
constrained by the critical rotation limit at their surface, setting a
restriction on their angular momentum content at core-collapse. We aim to
determine if this restriction plays a role in the spins of binary black hole
mergers.

Methods. We use detailed calculations of stripped stars with the MESA code at
low metallicities ($Z=Z_{\odot}/10$, $Z_{\odot}/50$ and $Z_{\odot}/250$) to
determine the dimensionless spins of black holes produced by critically
rotating stellar progenitors. To study how such progenitors can arise, we
consider their formation through chemically homogeneous evolution (CHE) in
binary stars. We use a semi-analytical model to study the physical processes
that determine the final angular momentum of CHE binaries, and compare our
results against available population synthesis models that rely on detailed
binary evolution calculations.

Results. We find that above black hole masses of $\simeq 25M_{\odot}$, the
dimensionless spin parameter of critically rotating stripped stars
($a=Jc/(GM^{2})$) is below unity. This results in an exclusion region at high
chirp masses and effective spins that cannot be populated by isolated binary
evolution. CHE can produce binaries where both black holes hit this limit,
producing a pile-up at the boundary of the excluded region. Highly spinning
black holes arise from very low-metallicity CHE systems with short delay
times, which merge at higher redshifts. On the other hand, the contribution of
CHE to merging binary black holes detected in the third observing run of the
LVK collaboration is expected to be dominated by systems with low spins
($\chi_{\mathrm{eff}}<0.5$) which merge near redshift zero. Owing to its
higher projected sensitivity and runtime, the fourth observing run of the LVK
collaboration can potentially place constraints on the high spin population
and the existence of a limit set by critical rotation.

###### Key Words.:

(Stars:) binaries (including multiple): close, Stars: black holes, Stars:
rotation, Gravitational waves

##  1 Introduction

Ever since the first observation of a merging binary black hole (BBH) (Abbott
et al., 2016), continuous improvements to ground-based interferometers have
resulted in weekly discoveries when current generation detectors are online.
As detector biases are well understood, with close to 100 detected sources
(Abbott et al., 2019, 2020b, 2021; Nitz et al., 2023) it has become possible
to place constraints on the intrinsic distribution of merging BBH masses and
spins at redshifts near zero (Abbott et al., 2023). Since the number of
sources is still limited the inference of intrinsic distributions mostly
relies on parameterized models, including features such as a bump or cutoff in
the mass distribution which is expected from pulsational pair-instability
supernovae (PPISNe, Belczynski et al. 2014; Marchant et al. 2016; Belczynski
et al. 2016; Woosley 2017; Fishbach & Holz 2017; Talbot & Thrane 2018;
Marchant et al. 2019; Farmer et al. 2019). Different parameterized
distributions can, however, lead to physically different interpretations. The
model used by the LIGO/Virgo/KAGRA collaboration for the distribution of black
hole spins indicated counter-aligned black hole spins (Abbott et al., 2023),
while Roulet et al. (2021); Galaudage et al. (2021); Callister et al. (2021);
Hoy et al. (2022) argued that this conclusion is dependent on the model
distribution used.

With an increasing number of measurements it is also becoming possible to
determine correlations between different BBH properties and to place
constraints on their joint distributions. For low-mass binaries the best
constrained intrinsic property is the chirp mass,

$\displaystyle
M_{\mathrm{chirp}}=\frac{(m_{1}m_{2})^{3/5}}{(m_{1}+m_{2})^{1/5}},$ (1)

where $m_{1}$ and $m_{2}$ are the individual component masses. Useful
constraints are also placed on both the effective spin $\chi_{\mathrm{eff}}$
and mass ratio $q$, defined as

$\displaystyle\chi_{\text{eff}}\equiv\frac{m_{1}a_{1}\cos\theta_{1}+m_{2}a_{2}\cos\theta_{2}}{m_{1}+m_{2}},\quad
q=\frac{m_{2}}{m_{1}}$ (2)

where $a_{1}$ and $a_{2}$ are the dimensionless spin parameters of each BH
while $\theta_{1}$ and $\theta_{2}$ are the angles between the individual spin
vectors and the orbital angular momentum. Despite there being a significant
correlated uncertainty between $\chi_{\mathrm{eff}}$ and $q$ for individual
detections (e.g., Hannam et al., 2013; Mandel & Smith, 2021), current
observations already suggest the presence of an intrinsic (anti)correlation
between these two properties, as well as a non trivial dependence of the
$\chi_{\mathrm{eff}}$ distribution with respect to $M_{\mathrm{chirp}}$
(Safarzadeh et al., 2020; Callister et al., 2021; Abbott et al., 2023;
Franciolini & Pani, 2022; Adamcewicz et al., 2023). Our understanding of these
distributions will rapidly evolve, as the currently ongoing O4 observing run
is expected to significantly increase the number of observed BBH mergers
(Abbott et al., 2020a).

The precise distribution of spins is critical to differentiate between various
formation scenarios of merging BBHs. Dynamical processes in dense environments
can result in an isotropic distribution of angles between the orbital plane
and the spins of each BH in contrast to near-aligned spins in binary formation
channels (Mandel & O’Shaughnessy, 2010; Rodriguez et al., 2016). High black
hole spins can also be indicative of hierarchichal BH mergers, as a generic
outcome from the merger of two similar-mass slowly-spinning BHs is the
formation of a BH with a dimensionless spin $\sim 0.7$ (Pretorius, 2005; Berti
& Volonteri, 2008; Gerosa & Berti, 2017; Rodriguez et al., 2019; Tagawa et
al., 2021; Gerosa & Fishbach, 2021). In isolated binary systems it has been
predicted that envelope stripping coupled with efficient angular momentum
transport between the stellar core and the outer layers should lead to small
aligned BH spins (Belczynski et al., 2020).

Merging BBHs with high spins ($\chi_{\mathrm{eff}}>0.5$) could also arise in
binary systems, either through tidal coupling in a short period post-common-
envelope binary (Qin et al., 2018; Zaldarriaga et al., 2018; Bavera et al.,
2020, 2021; Fuller & Lu, 2022; Ma & Fuller, 2023) or through chemically
homogeneous evolution (CHE). Large- scale flows are predicted to mix the
hydrogen burning core with the envelope in rapidly rotating massive stars,
preventing radial expansion in the CHE channel (Maeder, 1987). de Mink et al.
(2009) argued that this process operates in short period binaries (periods of
the order of one day) where rapid rotation arises from tidal synchronization.
If both components of a massive binary undergo CHE, it is then possible to
form a merging BBH without the need of an additional process to harden the
binary (Mandel & de Mink, 2016; Marchant et al., 2016; de Mink & Mandel, 2016;
du Buisson et al., 2020; Riley et al., 2021; Stevenson & Clarke, 2022). BBHs
formed through CHE can then produce pairs of high spin BHs (Marchant et al.,
2016; Fuller & Ma, 2019), although this is sensitive to mass loss after the
main sequence.

Various attempts have been made to constraint the individual contributions of
different formation channels to the observed sample of merging compact
objects. One possible approach to this is the use of individual predicted
distributions from different channels, which can then be scaled to attempt to
reproduce the observed distribution coupled with its observation biases (Zevin
et al., 2017; Bouffanais et al., 2019; Safarzadeh et al., 2020; Wong et al.,
2021; Zevin et al., 2021; Périgois et al., 2023). Alternatively, predicted
distributions from each channel can be used to infer individually what is the
likelihood that a single source is produced by a specific formation scenario
(e.g. Qin et al. 2022; Antonelli et al. 2023). However, conclusions from model
comparison can be severely affected if one lacks knowledge of some formation
channels (Qiu Cheng et al., 2023), and predicted distributions have large
uncertainties with different codes predicting drastically different outcomes
even for individual binary systems (e.g. Belczynski et al. 2022). It is then
critical to identify the qualitative features that arise from specific
formation channels both in masses and spins, independent of large
uncertainties in population modelling (e.g. van Son et al. 2022).

In this work we study the limit set by the critical rotation rate of a BH
progenitor on its resulting spin. In Section 2 we discuss the physics of this
limit and how it is expected to scale with stellar mass. In Section 3 we
connect this limit with evolutionary processes to form merging BBHs, by
developing a simple semi-analytical model that describes the evolution of CHE
binaries. We then compare our results to the detailed population synthesis
calculations of du Buisson et al. (2020) (hereafter dB20), which we have re-
analyzed to include predictions on the spins of BBHs produced via CHE. We
conclude our work and discuss its broader implications in Section 5.

##  2 Stellar critical rotation and black hole spins

In the past two decades, various results have pointed to the need of efficient
angular momentum transport in stellar interiors. This includes asteroseismic
measurements of the rotation rates of the cores of evolved stars (see Aerts et
al. 2019 for a recent review) as well as the angular momentum content of white
dwarfs and young neutron stars (Suijs et al., 2008). One potential mechanism
that allows for efficient angular momentum transport throughout the life of a
star is the Tayler-Spruit (TS) dynamo (Spruit, 2002). This dynamo process
amplifies seed magnetic fields through differential rotation in radiative
layers, leading to a stronger coupling of the stellar core and its outer
layers. Fuller et al. (2019) further developed the model for the TS dynamo,
and argued this process is more efficient than previously thought (we refer to
this model as the TSF dynamo). The model by Fuller et al. (2019) provides a
better fit to asteroseismic measurements of red giant and red clump stars, and
to the rotation rates of white dwarfs. It also leads to near solid-body
rotation of the hydrogen depleted cores.

The TSF dynamo reduces significantly the final rotation rates of massive
stripped stars, as it couples the outer mass-losing layers to the core that
remains until core-collapse (Fuller & Ma, 2019). If angular momentum transport
is efficient enough to make a stripped star remain near solid-body rotation
until late burning phases, the critical rotational frequency
$\Omega_{\mathrm{crit}}$ sets a limit on the maximum angular momentum of the
star,

$\displaystyle J_{\mathrm{max}}=I\Omega_{\mathrm{crit}},$ (3)

where $I$ is the moment of inertia. The critical rotation rate
$\Omega_{\mathrm{crit}}$ represents the limit at which centrifugal, radiation
and gravitational forces are perfectly balanced at the stellar equator. If the
star collapses directly into a BH, then its dimensionless BH spin is limited
to

$\displaystyle
a_{\mathrm{BH}}=\min\left(1,\frac{J_{\mathrm{max}}c}{M^{2}G}\right),$ (4)

where $M$ is the stellar mass. In practice, as the BH is assembled, if some
layers have enough angular momentum to support a disk at the innermost stable
circular orbit then outflows can remove the surplus angular momentum and form
a sub-critical BH (e.g. Batta & Ramirez-Ruiz 2019; Murguia-Berthier et al.
2020).

Owing to the large luminosities of massive stars, one needs to account for the
impact of radiation on the critical rotation rate. Ignoring gravity darkening,
the critical rate depends on the Eddington factor $\Gamma$ and the equatorial
radius $R_{\mathrm{eq}}$ as (Langer, 1997)

$\displaystyle\Omega_{\text{crit}}=\sqrt{\frac{GM(\Gamma-1)}{R_{\mathrm{eq}}^{3}}},\quad\Gamma=\frac{\kappa
L}{4\pi cGM},$ (5)

with $\kappa$ and $L$ being the surface opacity and luminosity, respectively.
The critical rotation rate including radiation is commonly referred to as the
$\Gamma-\Omega$ limit. In this section we aim to provide a qualitative
estimate of the impact of the $\Gamma-\Omega$ limit on BH spins. We will use
non-rotating stellar models to estimate $a_{\mathrm{BH}}$. Even if we were to
use rotationally deformed stellar models, the expression for the critical
rotation rate of Langer (1997) is only approximate as it ignores the impact of
gravity darkening (Maeder & Meynet, 2000; Sanyal et al., 2015; Fabry et al.,
2022).

The dependence of the limiting $a_{\mathrm{BH}}$ on mass can be understood by
expressing the moment of inertia as

$\displaystyle I=kMR^{2},$ (6)

where $k$ is the square of the ratio between the gyration radius and the
stellar radius. Here we ignore rotational deformation and use
$R_{\mathrm{eq}}=R$. If we take $k$ to be independent of mass at a given
evolutionary stage (which is the case for polytropic models with the same
polytropic index), and consider mass-radius and mass-luminosity relationships
of the form $R\propto M^{\beta}$ and $L\propto M^{\alpha}$ respectively, one
has that:

$\displaystyle\frac{J_{\mathrm{max}}c}{M^{2}G}\propto
M^{(\beta-1)/2}\sqrt{\Gamma-1},\quad\Gamma\propto M^{\alpha-1}.$ (7)

Power law relationships are not valid across all stellar masses, so the power
law exponents are normally also taken to be functions of mass. In particular,
the mass-luminosity exponent for main sequence stars is $\alpha>1$, while it
approaches unity at very high masses (e.g. Köhler et al. 2015). This leads to
$\Gamma\rightarrow 1$ as mass increases, which pushes down the spin limit in
Equation (7). Similarly, both stellar models and homology relationships show
that main sequence massive stars have $\beta<1$ (Kippenhahn et al., 2013),
which is also found from observations of eclipsing binaries (e.g. Torres et
al. 2010, Eker et al. 2018). Having $\beta<1$ indicates that the spin limit in
Equation (7) is reduced with mass. As more massive stars approach the
Eddington limit near their surface they might experience inflation, where a
small amount of mass expands to large radii (e.g. Sanyal et al. 2015). In this
case $\beta>1$ but also the approximation that $k$ is independent of mass is
not valid . If we want to understand in detail how the critical rotation limit
on spins works at various evolutionary stages and masses, we need to rely on
numerical calculations.

To make quantitative estimates we compute mock models of stripped stars using
version 23.05.1 of the MESA stellar evolution code111Input files to reproduce
our simulations together with our data are available at
https://zenodo.org/doi/10.5281/zenodo.10199794. (Paxton et al., 2011, 2013,
2015, 2018, 2019; Jermyn et al., 2023). Details about input microphysics are
described in Appendix A. We will also use these models to discuss the spins
produced by CHE stars, so evolution is initiated at the zero-age main-sequence
(ZAMS). We consider non-rotating models without mass loss which are forced to
evolve homogeneously until core helium depletion by enforcing a large
diffusion mixing coefficient throughout the entire star. Initial masses are
sampled logarithmically between $\log_{10}M/M_{\odot}=0.8$ ($6.3M_{\odot}$)
and $\log_{10}M/M_{\odot}=2.4$ ($251M_{\odot}$) and kept fixed during the
evolution. Three metallicites are sampled, $Z_{\odot}/10$, $Z_{\odot}/50$ and
$Z_{\odot}/250$ where the solar metallicity $Z_{\odot}=0.0142$ and the
relative metal fractions are adopted from Asplund et al. (2009). All
simulations are run until core carbon depletion.

As we do not include mass loss, the resulting models are not meant to
represent the evolution of a single star, but to be indicative of the
structure of a star with a given mass at a specific evolutionary stage. Since
CHE stars in binaries contract and tidally decouple after the main sequence,
stellar winds can efficiently spin them down and we do not expect them to
remain homogeneous. However, a stellar model evolved without mass loss will
produce a composition discontinuity between the CO core and the remaining
helium envelope, which is not present in mass-losing models. It is for this
reason that we keep our models homogeneous through helium burning, to better
represent the structure of a model at each evolutionary stage and metallicity.

Figure 1: The mock CHE models described in section 2 on the Hertzsprung-
Russell diagram. Colors indicate different evolutionary stages including the
zero-age main sequence (ZAMS), the terminal age main sequence (TAMS), the
helium ZAMS and TAMS, and core carbon depletion (C dep.). Different linestyles
indicate different values for the metallicity. Two evolutionary tracks with a
metallicity of $Z_{\odot}/50$ and masses of $10M_{\odot}$ and $100M_{\odot}$
are shown in gray for reference.

The position in the Hertzsprung-Russell diagram for the mock CHE models at
different evolutionary stages is illustrated in Figure 1. As the stars evolve,
they continuously increase their luminosity and (except for some of the most
massive models) continuously contract. At each of these stages we can consider
the maximum angular momentum that the star can contain under solid-body
rotation, $J_{\mathrm{max}}$, and compute the resulting maximum spin
$a_{\text{max}}=J_{\text{max}}c/GM^{2}$ a black hole would have were the star
to collapse at that stage. As shown in Figure 2 we find that for all
metallicities and masses $a_{\text{max}}$ decreases as the star evolves. In
particular, at core-carbon depletion, stars with masses above $25M_{\odot}$
would not be able to produce a maximally spinning BH. Barring the assumption
of solid-body rotation, which we discuss later, this is a feature that should
be common not only to CHE but to any formation scenario that involves the
formation of BHs through stripped envelope stars. This also includes stripped
stars that are tidally spun up by their companion (e.g. Qin et al. 2018). A
fit to $a_{\mathrm{max}}$ at core carbon depletion is provided in Appendix B.

The limit we find on BH spins is consistent with currently measured values in
X-ray binaries. The most massive BH with a dynamical mass estimate and a
constrained spin is Cyg X-1, with a mass of $21.2\pm 2.2M_{\odot}$ and
$a>0.9985$ (Miller-Jones et al., 2021). It might be possible in the near
future to identify more massive BHs in X-ray binaries, which would allow to
test our predicted spin limit with electromagnetic observations. On the other
hand, merging binary BHs with masses $>20M_{\odot}$ are now frequently
detected, possibly allowing a direct test. Even if critical rotation sets a
limit on the spins of BHs formed from stripped stars, for this spin limit to
appear in GW observations we require that the evolutionary processes which
form merging binary BHs also lead to critically rotating stars close to core-
collapse. In the context of long gamma-ray burst (lGRB) formation via the
collapsar model (Woosley, 1993), tidal spin-up of stripped stars by degenerate
companions has been suggested as a mechanism to form a rapidly spinning BH
(Podsiadlowski et al., 2004; Izzard et al., 2004). If such systems also result
in a compact object merger, then the second formed compact object would
dominate the effective spin of the system (Qin et al., 2018; Bavera et al.,
2020, 2021). Another mechanism to form high spin merging binary BHs is the CHE
scenario (Mandel & de Mink, 2016; de Mink & Mandel, 2016; Marchant et al.,
2016). In the next section we consider the CHE scenario, by means of a simple
semi-analytical model that lets us understand what limits the masses and spins
that can be produced by this channel.

Figure 2: Dimensionless spins corresponding to the angular momentum and mass
of a critically rotating star undergoing CHE at different evolutionary stages.
Values are estimated from the non-rotating models described in Section 2.

##  3 Semi-analytical model for chemically homogeneous evolution

To explore the possible parameter space in chirp mass and effective spin
covered by CHE, we construct a simple semi-analytical model. This model relies
on a few physical quantities obtained from stellar evolution models described
in Section 2 and provides clear insight into the physics that determine the
boundaries in BBH properties produced by CHE. For simplicity, we consider only
equal-mass binaries and assume solid-body rotation. Through multiple
evolutionary stages we need to account for the evolution of the angular
momentum of each component, as well as their masses and orbital periods. We
describe our assumptions at each stage next.

  1. 1.

TAMS: We ignore evolution during core hydrogen burning and assume as a
starting point that both stars have evolved chemically homogeneous and have
reached the TAMS while remaining tidally synchronized. The parameters that
describe this stage are the orbital period $P_{\mathrm{TAMS}}$ and the mass of
each star $M_{\mathrm{TAMS}}$. Given the assumption of synchronicity the
angular momentum of each star at this stage is equal to

$\displaystyle
J_{\mathrm{TAMS}}=\frac{2\pi}{P_{\mathrm{TAMS}}}I_{\mathrm{TAMS}}.$ (8)

The orbital period at TAMS cannot be arbitrarily small. We consider as a limit
the period at which the stars would be in contact222As a less restrictive
constraint one could consider the limit at which a contact binary at TAMS
would be just below the limit of L2 overflow.. For a star of mass $M_{1}$ with
a companion of mass $M_{2}$ this period can be computed from Kepler’s third
law and the separation at contact,

$\displaystyle
P_{\mathrm{contact}}=\sqrt{\frac{4\pi^{2}}{G(M_{1}+M_{2})}a_{\mathrm{contact}}^{3}}.$
(9)

The separation at contact is determined using the fit for the Roche lobe size
of Eggleton (1983) and the stellar radius of the CHE star at TAMS,

$\displaystyle
a_{\mathrm{contact}}=\frac{0.6q^{2/3}+\ln(1+q^{1/3})}{0.49q^{2/3}}R_{\mathrm{1,TAMS}},\quad
q\equiv\frac{M_{1}}{M_{2}}.$ (10)

  2. 2.

He-ZAMS: After the main sequence, a CHE star will contract until it reaches
thermal equilibrium at the beginning of core-helium burning, the He-ZAMS. We
assume tidal synchronization stops during this short contraction phase and
that the angular frequency of the star evolves by conserving the angular
momentum at TAMS. If contraction is homologous, meaning the radial position of
each mass element contracts by the same fraction $\alpha$, then the angular
frequency of the star increases as $\Omega\propto\alpha^{-2}$. In contrast,
ignoring changes in surface opacity or luminosity, the critical rotation rate
increases as $\Omega_{\mathrm{crit}}\propto\alpha^{-3/2}$. The contracting
star is then expected to evolve closer to critical rotation
($\Omega/\Omega_{\mathrm{crit}}\propto\alpha^{-1/2}$). In practice contraction
is not exactly homologous, and the critical rotation rate will be modified as
the luminosity and surface opacity change, but this illustrates that during
contraction after TAMS stars can evolve towards critical rotation.

To determine the angular momentum at the He-ZAMS we then need to account for
the maximum angular momentum $J_{\mathrm{He-ZAMS,max}}$ that the star can
possibly have while rotating at its critical speed. We simply limit the
angular momentum at this stage to this value,

$\displaystyle J_{\mathrm{He-ZAMS}}=\min(J_{\mathrm{TAMS}},J_{\text{He-ZAMS,
max}}).$ (11)

The loss of excess angular momentum needs to be accompanied by mass loss.
However, if ejection happens through the formation of a decretion disk
(similar to Be-stars, e.g. Rivinius et al. 2013) then significant amounts of
angular momentum can be removed with a small amount of mass ejected. We ignore
the uncertain mass loss and orbital period evolution during the contraction
from TAMS to He-ZAMS,

$\displaystyle M_{\mathrm{He-ZAMS}}=M_{\mathrm{TAMS}},\quad P_{\mathrm{He-
ZAMS}}=P_{\mathrm{TAMS}}.$ (12)

  3. 3.

Core helium burning: During the post-MS evolution the components will undergo
mass loss, possibly as Wolf-Rayet stars. We parameterize the total fraction of
mass lost by a free parameter $f_{\mathrm{WR}}$ defined as

$\displaystyle f_{\mathrm{WR}}\equiv\frac{M_{\mathrm{He-TAMS}}}{M_{\mathrm{He-
ZAMS}}}.$ (13)

This is such that $f_{\mathrm{WR}}=1$ represents no mass loss, while
$f_{\mathrm{WR}}=0$ would indicate that the entire stellar mass is removed.

As mass is lost the orbital period of the binary will also change. Under the
approximation where the stellar winds remove orbital angular momentum equal to
the specific orbital angular momentum of the corresponding component, the
orbital period satisfies

$\displaystyle\frac{P_{\mathrm{He-TAMS}}}{P_{\mathrm{He-
ZAMS}}}=\left(\frac{M_{\mathrm{He-ZAMS}}}{M_{\mathrm{He-
TAMS}}}\right)^{2}=f_{\mathrm{WR}}^{-2}.$ (14)

Since $0<f_{\mathrm{WR}}<1$, mass loss will lead to orbital expansion.

To determine the evolution of each star’s angular momentum in connection to
mass loss, we assume mass loss removes spin angular momentum equal to that of
spherical shells of the stellar radius rotating at the surface angular
frequency,

$\displaystyle\dot{J}=i_{\mathrm{surf}}\Omega\dot{M},\quad
i_{\mathrm{surf}}=\frac{2}{3}R^{2}.$ (15)

Writing the moment of inertia of the star in terms of $k$ (which is the
squared ratio between the gyration radius and the stellar radius) we obtain
that

$\displaystyle\frac{\dot{J}}{J}=\frac{2}{3k}\left(\frac{\dot{M}}{M}\right).$
(16)

Our simulations indicate that during core helium burning $k$ covers the range
$\sim 0.08-0.1$. Ignoring the time dependence of $k$, the previous equation
can be integrated to provide the angular momentum of the star at the He-TAMS,

$\displaystyle\frac{J_{\mathrm{He-TAMS}}}{J_{\mathrm{He-
ZAMS}}}=\left(\frac{M_{\mathrm{He-TAMS}}}{M_{\mathrm{He-
ZAMS}}}\right)^{2/(3k)}=f_{\mathrm{WR}}^{2/(3k)}.$ (17)

This can be related to the dimensionless spin of the star, i.e., the spin of a
BH which a star would produce through direct collapse, $a=Jc/M^{2}G$. In terms
of $f_{\mathrm{WR}}$ we have that

$\displaystyle\frac{a_{\mathrm{He-TAMS}}}{a_{\mathrm{He-
ZAMS}}}=f_{\mathrm{WR}}^{-2+2/(3k)}.$ (18)

This equation indicates that the final compact object spin is reduced under
our assumption of solid-body rotation as long as $k<1/3$. Owing to the central
concentration of mass in stars, this is generally the case.

  4. 4.

Carbon depletion and BBH formation: Similar to the evolution from TAMS to He-
ZAMS, we limit the angular momentum of a star to its maximum value at critical
rotation,

$\displaystyle J_{\text{C. dep}}=\min(J_{\mathrm{He-TAMS}},J_{\text{C. dep,
max}}),$ (19)

and we keep the mass and orbital period from the He-TAMS,

$\displaystyle M_{\text{C. dep}}=M_{\mathrm{He-TAMS}},\quad P_{\text{C.
dep}}=P_{\mathrm{He-TAMS}}.$ (20)

We assume direct collapse to a BH, ignoring any loss of mass or angular
momentum at stages beyond carbon depletion. For the mass range we are
concerned with this is supported both by theory (e.g. Fryer 1999) as well as
observations of BH binaries with massive companions such as Cyg X-1 (Mirabel &
Rodrigues, 2003) and VFTS-243 (Shenar et al., 2022). The resulting mass and
spin of the BHs is finally computed as

$\displaystyle M_{\mathrm{BH}}=M_{\text{C. dep}},\quad
a_{\mathrm{BH}}=\min\left(1,\frac{J_{\text{C. dep}}c}{M_{\text{C.
dep}}^{2}G}\right),$ (21)

while the orbital period of the BBH is taken to be that at carbon depletion.
The delay time until BBH merger is then computed from the results of Peters
(1964) for the case of a circular orbit. As a simple model for the impact of
(P)PISNe we ignore any BH that would form within the predicted PISN gap,
taking the lower limit of the gap to be $M_{\mathrm{BH}}=43.94M_{\odot}$ from
Marchant et al. (2019).

Our model then relies on three input parameters: the component masses at TAMS,
the orbital period at TAMS, and the fraction of mass lost after the main
sequence ($M_{\mathrm{TAMS}}$, $P_{\mathrm{TAMS}}$ and (1-$f_{\mathrm{WR}}$)
respectively). The model is semi-analytical as it does require results from
detailed calculations of stellar evolution. Specifically, we require, as a
function of mass: the radius and moment of inertia at TAMS
($R_{\mathrm{TAMS}}$ and $I_{\mathrm{TAMS}}$), the maximum angular momenta at
He-ZAMS and carbon depletion ($J_{\text{He-ZAMS, max}}$ and $J_{\text{C. dep,
max}}$), and the squared ratio between the gyration and surface radius $k$.
All of these quantities, which are determined by MESA calculations in our
work, are metallicity dependent.

Figure 3: Resulting BH masses and spins of CHE stars using a semi-analytical
model (see Section 3). The model assumes a mass ratio of unity, so two black
holes with identical properties are produced. We consider the outcome for
three masses at TAMS ($25M_{\odot}$, $40M_{\odot}$ and $60M_{\odot}$) and
three orbital periods at TAMS scaled from the period at which the binary would
be in contact. Each line represents the results for a fixed
$M_{\mathrm{TAMS}}$ and $P_{\mathrm{TAMS}}$ with a variable $f_{\mathrm{WR}}$
(where $M_{\mathrm{BH}}=f_{\mathrm{WR}}M_{\mathrm{TAMS}}$). Cases where the
resulting BBH would be too wide to merge within the age of the Universe are
marked with a dotted line. Note that for $f_{\mathrm{WR}}$ close to unity
lines with different values of $P_{\text{TAMS}}$ converge to the same spins at
any given mass. The effect of (P)PISNe is ignored in this Figure. Precise
boundaries depend on the metallicity used to determine stellar properties in
the semi-analytical model used, which in this case where taken from
calculations at $Z_{\odot}/50$.

We focus on predictions at the metallicity $Z=Z_{\odot}/50$ to describe the
qualitative results. Results from the other computed metallicities
($Z_{\odot}/10$ and $Z_{\odot}/250$) are qualitatively similar. As we find
that $k$ covers values between $\sim 0.08-0.1$ during core helium burning, we
adopt $k=0.09$. By considering different values for $M_{\mathrm{TAMS}}$ and
$P_{\mathrm{TAMS}}$ we can explore what sets the highest and lowest spins
produced for different BH masses. This is done in Figure 3, where for each
choice of $M_{\text{TAMS}}$ and $P_{\text{TAMS}}$ we vary continuously
$f_{\mathrm{WR}}$ between unity and zero. The right-most point of each curve
represents the case of no mass loss ($f_{\mathrm{WR}}=1$). The most apparent
feature in Figure 3 is an exclusion region for high spins and high BH masses.
This limit is set by critical rotation at carbon depletion (i.e. the
$\Gamma-\Omega$ limit, see Section 2), and determines the final BH spin as a
function of BH mass, independent of the initial orbital period, when WR mass
loss is negligible. The $\Gamma-\Omega$ limit is not restricted to CHE
evolution, but should be a common feature of evolutionary channels where BHs
are formed from helium stars.

More subtle are the conditions that define the minimum spins that can be
achieved at any given mass from CHE. There are two ways in which angular
momentum could be reduced: either reach TAMS at a higher orbital period
(leading to slower rotation due to synchronization), or lose mass to remove
angular momentum after TAMS. Since we still need the resulting BBH to merge,
we can question which of these two options results in a smaller final spin. Is
it better to start at contact at TAMS and then lose as much mass as possible
before the orbit widens too much? Or can lower spins be achieved by starting
with an orbit at TAMS as wide as possible to yield a merging BBH without any
additional mass loss? Figure 3 shows that the former is the case. For each
$M_{\mathrm{TAMS}}$ considered, the lowest spin achievable while still
resulting in a merging BBH comes from binaries that are at contact in TAMS. An
important caveat is that our choice of $f_{\mathrm{WR}}$ is arbitrary, while
in practice the fraction of mass lost will be determined by the metallicity
and mass at the He-ZAMS.

The upper and lower bounds on the spins discussed so far are illustrated in
Figure 4. Two additional bounds are included, one accounting for the (P)PISNe
gap following the results of Marchant et al. (2019), and another one coming
from the limits of CHE evolution. CHE is expected to be a process favored at
higher masses (Maeder, 1987), such that we cannot have an arbitrarily low
$M_{\text{TAMS}}$. In figure 4 we illustrate this by drawing an exclusion
region that corresponds to $M_{\text{TAMS}}\leq 30M_{\odot}$. These boundaries
can be compared to the results of detailed simulations to see if this simple
model can reproduce the predicted features, which we do in the following
section. We point out however that various physical effects are not included
in this model. In particular, mass loss is expected to increase at higher
masses and metallicities. This also implies that at the high metallicities
where we would expect significant post-MS mass loss, mass loss will have
potentially widened the orbit of the stars during the MS itself. This
restricts the minimum period at TAMS and thus the lower limit on the final
spins. Despite this, even in detailed models the impact of mass loss on the
orbital period is approximate, and orbital widening can be reduced or even
inverted by tidal coupling of the wind (Brookshaw & Tavani, 1993; MacLeod &
Loeb, 2020; Schrøder et al., 2021).

Figure 4: Physical processes setting the boundaries for BH masses and spins
for the CHE formation scenario. The four colored areas indicate excluded
regions due to four different physical processes. In this figure an assumed
limit for the $M_{\text{TAMS}}$ of CHE stars is set at $30M_{\odot}$, which
sets the boundary at the left of the Figure. Precise boundaries depend on the
metallicity used to determine stellar properties in the semi-analytical model
used, which in this case where taken from calculations at $Z_{\odot}/50$.

##  4 Comparison to population synthesis results

The semianalytical model presented in the previous Section provides an overall
idea of the limits of CHE in terms of masses and effective spins. However, as
we vary multiple parameters without accounting for their possible correlation
(such as the fraction of mass lost after the MS and the period at TAMS), it is
instructive to compare our results with detailed population synthesis
calculations. In this Section we re-analyze the population synthesis results
of  by including information on stellar spins and compare them against our
results.

###  4.1 The population synthesis simulations of dB20

The work of dB20 already provided predictions for the rates and the mass
distribution of merging BBHs formed through CHE. These results were based on a
large grid of MESA models, exploring the parameter space of initial masses,
orbital periods and metallicity. For computational reasons the simulations of
dB20 were restricted to a mass ratio of $q=1$, which was taken to be
representative of the evolution of binaries in the range $q=0.8-1$. In order
to provide rate predictions this grid of detailed binary simulations was
combined with the theoretical predictions of Taylor & Kobayashi (2015) for the
cosmic star formation history, accounting in particular for the distribution
of formed stars at different metallicities. The results of dB20 are provided
in the form of a comoving box simulation, indicating the formation redshift of
all BBHs formed within it through CHE.

Consistent with the results of Marchant et al. (2016), for any given
metallicity the MESA simulations of dB20 show a tight correlation between the
final orbital period and mass of a binary. This is a consequence of wind mass
loss, which at high metallicities lowers the total mass and, under the
assumption of Jeans’ mode mass loss (stellar winds removing orbital angular
momentum equal to that of each component), result in orbital widening. At a
fixed metallicity and chirp mass, the small spread in orbital periods is
connected to the narrow range in initial orbital periods that are expected to
lead to CHE (Mandel & de Mink, 2016; Marchant et al., 2016; Song et al., 2016;
Riley et al., 2021). For BH formation through direct collapse, this leads to a
similar relationship between the period and the chirp mass of the binary BH
(see top panel of Figure 5).

Figure 5: (Top) Final orbital period and chirp mass in the MESA simulations of
dB20, assuming direct collapse into a BH (with its mass at carbon depletion)
and ignoring mass loss due to (P)PISNe. Each dot indicates the final result of
one CHE simulation. Systems above the dashed gray line have merger times
longer than the age of the Universe and are marked with a lighter color.
Results are shown for their grids with metallicities of $\log_{10}Z=-2.375$,
$2.5$, $-2.75$, $-3$, $-3.5$, $-4$ and $-5$, with the labels indicating
approximately their values in terms of the solar metallicity (taken here to be
$Z_{\odot}=0.017$ as in Grevesse et al. (1996) for consistency with the
definition used in dB20). (Bottom) same as the top panel but after including
the impact of (P)PISNe as described in Marchant et al. (2019).

Including the effect of (P)PISNe produces a more complex relationship between
the chirp mass and the orbital period of BBHs formed through CHE. Mass loss
due to PPISNe was included in the results of dB20 by interpolating the pre-
PPISN to BH mass relationship of Marchant et al. (2019). The change in orbital
properties was determined by assuming mass loss was instantaneous, followed by
circularisation of the orbit at constant orbital angular momentum prior to BH
formation. The bottom panel of Figure 5 shows the resulting final orbital
periods and chirp masses after the inclusion of mass loss due to PPISNe. As
the relationship between pre-PPISNe mass and BH mass is non-monotonic, the
simple relationship between $M_{\text{chirp}}$ and the final orbital period
does not hold anymore. Owing to the ocurrence of PISNe, no BHs are produced
between $44M_{\odot}$ and $123M_{\odot}$, which constitutes the PISN mass gap
according to Marchant et al. (2019).

The resulting spins of BBHs formed through CHE were not discussed by dB20. In
all simulations that are unaffected by (P)PISNe we have assumed all angular
momentum and mass collapse into the BH, resulting in a black hole spin of

$\displaystyle a_{\text{BH}}=\frac{Jc}{M^{2}G},$ (22)

where $J$ and $M$ are the final pre-collapse spin angular momentum and mass.
Measurements of individual spins in merging binary BHs have large errors, but
the effective spin is much better constrained. As dB20 only consider equal
mass binaries with aligned spins, $\chi_{\mathrm{eff}}$ is equal to
$a_{\text{BH}}$ in their simulations. For cases undergoing PPISNe we use the
angular momentum profile of the star at the end point of the calculation,
which in dB20 corresponds to either core-carbon depletion or the onset of
collapse due to pair-instability. Taking the final expected BH mass
$M_{\mathrm{BH}}$ from Marchant et al. (2019) we consider only the angular
momentum $J(M_{\text{BH}})$ contained within the corresponding mass
coordinate. The final spin is then computed as

$\displaystyle a_{\text{BH,
PPISN}}=\frac{J(M_{\text{BH}})c}{M_{\text{BH}}^{2}G}.$ (23)

This assumes the mass ejection from PPISN is instantaneous. One caveat of this
approach is that the resulting BH spin can exceed unity, indicating the
potential formation of a GRB through the collapsar model (Woosley, 1993;
MacFadyen & Woosley, 1999). For simplicity we simply cap the dimensionless
spin at unity. More detailed models to relate the angular momentum
distribution of a stellar model to the final mass and spin of the BH are
available (e.g., Batta & Ramirez-Ruiz 2019), which consistently produce sub-
critical spins.

The resulting chirp masses and effective spins from the MESA models of dB20
are shown in Figure 6. Without the inclusion of (P)PISNe, these two quantities
follow clear relationships at each metallicity, with higher chirp masses
resulting in lower spins. This is a consequence of wind mass loss after TAMS,
lowering the spin of the star before BH formation. At each metallicity the
lowest spins produced for merging BBHs correspond to the point where the
resulting delay time is equal to the age of the Universe, an effect we
explained in Section 3. At very low metallicities ($\leq Z_{\odot}/200$) the
MESA simulations produce low spins that break from the general trend. We have
identified that in these cases, owing to a smaller relative difference between
the size of the stars at TAMS and the He-ZAMS, stars remain synchronous after
the end of the main sequence. This leads to angular momentum being removed
from each component by tides on the He main sequence, as stars spin up during
the contraction phase after core-hydrogen depletion. This effect is
potentially spurious, as the model used for tidal evolution is based on the
dynamical tide process of Zahn (1977) using pre-computed properties of main-
sequence stars (see Paxton et al. 2015 for details). Recently Ma & Fuller
(2023) have shown that tidal synchronization likely only operates in stripped
stars for very short orbital periods, which are forbidden in CHE evolution as
the orbit needs to be wide enough to fit the stars during their main sequence
evolution.

Figure 6: Same as Figure 5 but showing the resulting chirp masses and
effective spins of BBHs formed through CHE in the MESA simulations of dB20.

For this work we have also recomputed the detection probabilities from dB20
using updated noise curves. We have used noise curves corresponding to the O3
measured sensitivity, and the O4 and O5/A+ expected sensitivities as described
in Abbott et al. (2020a). In particular, for O3 we have used the measured
noise curve based on three months of data from the LIGO-Hanford detector,
while for O4 we have used the pessimistic noise curve estimate with a range
for binary neutron star detection of 160 Mpc333These noise curves are
available at https://dcc.ligo.org/LIGO-T2000012/public. We have also included
the effect of spins on detectability by using the IMPRPhenomB waveform model
with non-precessing spins from Ajith et al. (2011).

###  4.2 Rates of high spin mergers from chemically homogeneous evolution

Using the insight developed from the semi-analytical model of Section 3 we re-
analyze the rate predictions of dB20 for CHE evolution while considering the
predicted spins for these models. Before discussing rates for specific
detector sensitivities we first compare the properties of our semi-analytical
model with the results of dB20 for a detector with infinite sensitivity (which
we refer to as “cosmological” rates). Figure 7 illustrates how the
cosmological BBH mergers are distributed in effective spins and chirp masses.
As PPISNe produce a pile-up of systems near the PISN gap (Belczynski et al.,
2016), we see that the distribution in Figure 7 peaks at
$M_{\mathrm{chirp}}\sim 38M_{\odot}$. An additional pile-up is seen at the top
of the distribution, coming from the spins of BHs being limited by the
critical rotation of their progenitor stars. For comparison, we also consider
the predicted region for CHE evolution shown in Figure 4. Although the semi-
analytical model qualitatively reproduces various features of the detailed
simulations, it is shifted significantly towards lower effective spins. The
mergers with $\chi_{\text{eff}}\lesssim 0.2$ and $30\lesssim
M_{\mathrm{chirp}}\lesssim 40$ correspond to the very metal poor systems that
undergo tidal synchronization after TAMS. As discussed in Section 4.1 the low
spins in these very metal-poor simulations could be spurious.

Figure 7: Distribution of the rate of merging BBHs per unit chirp mass and
unit effective spin that is observable from Earth with a detector with
infinite sensitivity, derived from the simulations of dB20. We refer to these
rates as “cosmological”. Numbers are normalized to the maximum value of the
differential rate
$\text{d}R/\text{d}M_{\mathrm{chirp}}/\mathrm{d}\chi_{\mathrm{eff}}$. The two
panels highlight predicted systems below and above the PISN gap. The white
dashed contour corresponds to the predicted region below the PISN gap that
arises from CHE, computed with our semi-analytical model in Figure 4

One of the reasons why the simulations of dB20 and our semi-analytical model
make a different prediction for the minimum spins produced by CHE was
discussed in Section 3. The lowest spins predicted arise from binaries that
are at contact at TAMS and lose a significant amount of mass during the post-
MS. But systems that would undergo significant post-MS mass loss would also be
expected to be widened by mass loss before TAMS, so the two free parameters in
the semi-analytical model are not really independent. However, another source
of discrepancy comes from angular momentum coupling and how close to solid-
body rotation the MESA simulations of dB20 are. Differential rotation also
modifies how the $\Gamma-\Omega$ limit translates into an upper limit on BH
spins, as it allows for a larger angular momentum budget in the stellar core
without exceeding critical rotation at the surface.

To test the impact of angular momentum transport in our results, we have
recomputed a simulation from dB20 where we replaced the Tayler-Spruit dynamo
with the modified version from Fuller et al. (2019), as implemented by Fuller
& Lu (2022). We also compute a model with an artificially high angular
momentum diffusivity in order to enforce solid-body rotation throughout. The
simulation chosen corresponds to the merging BBH with the lowest chirp mass
produced in the MESA grids of dB20, with $M_{\mathrm{chirp}}=16M_{\odot}$ and
$\chi_{\mathrm{eff}}=0.44$. The resulting BH spin from these simulations is
shown in Figure 8. Throughout most of the evolution of the binary, all three
simulations match closely, indicating the original calculation was near solid-
body rotation. However, the simulations diverege during core-helium burning,
with the original calculation from dB20 resulting in a much higher black hole
spin. The simulation using the modified Tayler-Spruit dynamo of Fuller et al.
(2019) remains almost identical to the one where solid-body rotation is
enforced, producing a BBH with $\chi_{\mathrm{eff}}=0.25$.

Figure 8: Mass of each component versus instantaneous dimensionless spin
during the evolution of an equal mass binary undergoing CHE. Shown are the
original simulation from dB20 using the Tayler-Spruit (normal TS), one
recomputed using the modified version of the Tayler-Spruit dynamo from Fuller
et al. (2019) (TSF) and one where solid-body rotation is enforced. These
simulations where computed with the same setup and MESA version (11701) used
by dB20.

Keeping in mind that the results of dB20 might overestimate the spins of the
resulting merging BBHs, we analyze the range of spins covered by their
simulations of CHE. To dissect how different systems contribute, we
distinguish the rate between regions of low mass
($M_{\mathrm{chirp}}<30M_{\odot}$), high mass
($30M_{\odot}<M_{\mathrm{chirp}}<60M_{\odot}$), and systems above the PISN gap
($M_{\mathrm{chirp}}>60M_{\odot}$). Similarly, we distinguish between systems
with low ($\chi_{\text{eff}}<0.5$) and high ($\chi_{\text{eff}}>0.5$)
effective spins. The resulting merger rates are given in Table 1. Focusing on
systems below the PISN gap, we find that the dominant contribution ($41\%$)
comes from the high-spin and low $M_{\mathrm{chirp}}$ regime, with high-spin
systems constituting the majority of the rate ($61\%$) below the PISN gap.
This matches the general view that CHE produces mostly high-spin BHs.

Table 1: Predicted detection rates (per year of detector time) for BBH mergers produced from CHE, derived from the simulations of dB20. Row marked as ”Cosmological” represents the full rate observable from Earth with a detector of infinite sensitivity. Column marked ”Full” shows the full rate for each detector, while the other columns represent subsets of the full population filtered by their effective spins and masses. |  | $\chi_{\mathrm{eff}}<0.5,$ | $\chi_{\mathrm{eff}}>0.5,$ | $\chi_{\mathrm{eff}}<0.5,$ | $\chi_{\mathrm{eff}}>0.5,$ |
---|---|---|---|---|---|---
Detector | Full | $\displaystyle\frac{M_{\mathrm{chirp}}}{M_{\odot}}<30$ | $\displaystyle\frac{M_{\mathrm{chirp}}}{M_{\odot}}<30$ | $\displaystyle 30<\frac{M_{\mathrm{chirp}}}{M_{\odot}}<60$ | $\displaystyle 30<\frac{M_{\mathrm{chirp}}}{M_{\odot}}<60$ | $\displaystyle\frac{M_{\mathrm{chirp}}}{M_{\odot}}>60$
Cosmological | 8600 | 1200 | 2800 | 1500 | 1300 | 1900
O3 | 80 | 31 | 13 | 20 | 7.5 | 7.3
O4 | 290 | 100 | 61 | 68 | 33 | 30
O5/A+ | 1700 | 420 | 550 | 320 | 260 | 170

Most high-spin systems come from low metallicity ($Z\lesssim Z_{\odot}/10$)
environments that experience little mass loss and have short merger times (see
Figures 5 and 6). One way to visualize this is by considering the evolution of
the comoving merger rate as a function of redshift, which is done in Figure 9.
We can see that at redshift zero the predicted volumetric rate is dominated by
low-spin systems below the PISN gap, accounting for 84% of the local rate. It
is only above redshift $1.5$ that high-spin systems constitute the majority of
the comoving rate. Observations done with current detectors are not sensitive
beyond a redshift of $z\sim 1$, so we expect significant differences between
the rate of high-spin mergers in the observable Universe, against those
actually observed with current detectors. Another important feature of Figure
9 is that systems predicted above the PISN gap become dominant at redshift
$\simeq 5$. Predictions above the PISN gap are one of the most uncertain
features of the model, as they require the formation of binary stars in excess
of $100M_{\odot}$ for each component (Heger & Woosley, 2002). However, this
makes mergers above the gap an interesting target for 3rd generation detectors
which will be sensitive to such massive mergers at redshifts beyond 10 (Hall &
Evans, 2019).

Figure 9: Evolution of the BBH merger rate and of different populations of
merger events predicted by dB20 for CHE. (Top) Fraction of the total rate for
specific constraints on the chirp mass and effective spins of BBH mergers.
(Bottom) Comoving merger rate corresponding to the same constraints in chirp
mass and effective spins. The dotted line indicates the cosmological rate
resulting from the addition of all others.

As current ground based GW detectors are not sensitive to high redshift
mergers, we expect the observed population of BBHs formed through CHE to be
dominated by low-spin systems. Figure 10 shows the detectable rate
distribution as a function of $M_{\mathrm{chirp}}$ and $\chi_{\text{eff}}$
with the measured sensitivity for O3 as well as projected sensitivities for O4
and O5 (which is planned to include the A+ upgrade). Owing to the lack of
sensitivity at redshifts above unity, both for O3 and O4 we find that the
majority of observed BBHs would populate a narrow band going from
$M_{\mathrm{chirp}}\simeq 18M_{\odot}$ and $\chi_{\mathrm{eff}}\simeq 0.45$
down to $M_{\mathrm{chirp}}\simeq 38$ and $\chi_{\mathrm{eff}}\simeq 0.25$.
This corresponds to the limit for which, at any given metallicity, BBHs that
are too wide to merge in the age of the Universe are produced. Compared to the
currently observed population we find that this band lies at higher
$\chi_{\mathrm{eff}}$ than the majority of observed sources. As shown in Table
1, for the sensitivity of O3 we would expect 30% of observed sources produced
by CHE below the PISN gap to have $\chi_{\mathrm{eff}}>0.5$, while only a
single source has been detected with both median component masses below
$60M_{\odot}$ and a median effective spin above $0.5$ (GW200308_173609). The
predicted $\chi_{\mathrm{eff}}$ band of observable CHE mergers would be
lowered by about a half if we assumed the enhanced Tayler-Spruit dynamo of
Fuller et al. (2019).

As the sensitivity of detectors increases in the future, the contribution from
highly-spinning mergers produced through CHE should become more important in
observations. As shown in Figure 10 and Table 1, high-spin
($\chi_{\text{eff}}>0.5$) and low-spin ($\chi_{\text{eff}}>0.5$) CHE BBHs are
predicted make similar contributions to observations during the fifth
observing run, with 52% of mergers detected below the PISN gap corresponding
to the high-spin population. This indicates that direct constraints on the
evolution of the merger rate as a function of redshift (including
subpopulations of merging BBHs) will be possible before the onset of third-
generation detections. For instance, current observations show an indication
that the distribution of $\chi_{\mathrm{eff}}$ broadens with redshift
(Biscoveanu et al., 2022).

Our predicted observable rates for merging BHs formed by CHE would require
that almost all GW sources observed to date are a consequence of CHE (see also
Stevenson & Clarke 2022). Predicted rates are also in direct contradiction for
mergers with both components above the PISN gap, which have not been observed
so far. The third observing run of the LVK collaboration, for which detections
are described in GWTC-2 and GWTC-3 (Abbott et al., 2020b, 2021), accumulated a
total of 275 days of observations with at least two detectors. According to
our results for O3, this would translate into $5.5$ detections of mergers
above the gap. The lack of such detections could point to weaknesses in the
CHE model, but also in our understanding of mass loss rates at very low
metallicities and of the distribution of orbital properties of very massive
binaries. However, we do find that BHs formed above the gap should have spins
$a<0.25$, owing to the limit set on their progenitors by critical rotation.
This should be the case for any BH formed above the gap from a stripped-
envelope star (not just through CHE) and should serve to clearly distinguish
them from objects formed through hierarchical mergers, which are expected to
have spins $a\gtrsim 0.7$ (Pretorius, 2005; Berti & Volonteri, 2008; Gerosa &
Berti, 2017).

Figure 10: Observable rate of BBH mergers considering the detector sensitivity
of O3, O4 and O5 for the LIGO detectors. The distribution of the rate is shown
as a function of $M_{\mathrm{chirp}}$ and $\chi_{\text{eff}}$. Data points
indicate all observations made so far as reported in the current gravitational
wave transient catalogues (GWTC-1, GWTC-2 and GWTC-3). For each observation
the median values of $M_{\mathrm{chirp}}$ and $\chi_{\mathrm{eff}}$ are shown
together with their 90% credible intervals.

##  5 Conclusions

In this work we describe an upper limit on the spins of merging BBHs formed
through binary evolution. This limit arises from the critical rotation rate of
stripped star progenitors before their collapse, preventing the formation of
BHs with spins near unity for BH masses above $\gtrsim 25M_{\odot}$. To study
how this limit arises through binary evolution we have considered BBH
formation through CHE, a process that allows for the formation of rapidly
rotating BH progenitors. We have analyzed the BH spins produced by CHE with a
semianalytical model, as well as by reanalyzing the population synthesis
results of dB20 by including information on BH spins formed through CHE. We
compare our results to current observations and provide predictions for
upcoming observing runs of ground-based GW detectors.

Based on our semi-analytical model, we identified three different processes
(in addition to the PISNe gap) that restrict the range of masses and effective
spins that can be achieved through CHE (see Figure 4). The lowest achievable
spins come from CHE binaries with delay times close to the age of the
Universe, making this population favored in current observations as they merge
at low redshifts. Another restriction is set by the minimum mass for CHE to
operate, resulting in a boundary with increasing $\chi_{\text{eff}}$ as a
function of $M_{\text{chirp}}$. Finally, we find that CHE binaries that
experience little mass loss after the main sequence have spins determined by
the critical rotation rate of their progenitors, which is determined by a
combination of gravity and radiation pressure for massive, luminous stars.
Identifying any of these features in GW observations would provide important
insight into stellar evolution processes that remain poorly constrained.
Although we have studied only the evolution of BBH progenitors through CHE,
the spin limit due to critically rotating BH progenitors should be shared by
any formation channel for BBHs that involves the collapse of stars without a
hydrogen envelope.

Comparing our analytical models with the simulations of dB20 we find that,
although qualitatively the expected boundaries in $M_{\text{chirp}}$ and
$\chi_{\text{eff}}$ are present, the final spins are higher in the detailed
binary simulations. In part we attribute this to angular momentum transport,
as simulations with more efficient coupling between their cores and outer
layers can lower the final spins by a factor of $\sim 2$. Coupling is expected
to impact both the lower and upper limits of spins that can be produced by
CHE, with stronger coupling resulting in lower final spins owing to increased
loss of angular momentum through stellar winds. The simulations of dB20
predict that during the ongoing O4 run (and also previous observing runs) the
majority of detected BBHs formed through CHE would have effective spins
smaller than $0.5$. If the angular momentum coupling is strong, the majority
of sources should have spins below $\sim 0.25$. Only in future science runs
will the high-spin population (which merges promptly at high redshifts) become
dominant.

The upper limit on BH spins we have identified also has consequences for the
formation of gamma-ray bursts from collapsing stars through the collapsar
model (Woosley, 1993; MacFadyen & Woosley, 1999). If BHs with spins near unity
are required to produce a highly relativistic jet through the Blanford-Znajek
process (Blandford & Znajek, 1977; McKinney, 2005), then the progenitors of
lGRBs formed by stripped stars would be limited to pre-collapse masses
$\lesssim 20M_{\odot}$. Only broad-lined Ic supernovae have been associated
with lGRB events from massive star collapses (Galama et al., 1998), possibly
indicating that only stripped stars are their progenitors. The precise limits
on progenitor masses are sensitive to angular momentum coupling, and future GW
observations can help to identify it. Constraining which stellar progenitors
can produce lGRBs is also critical to link observations of high-spin merging
BBHs with the observed rates of lGRBs (Bavera et al., 2022).

Finally, a detection with both BHs unambiguously above the PISN gap is still
missing and in contradiction with the results of dB20. This could be related
to our lack of knowledge of the distribution of binary properties at both very
high masses and low metallicities. Perhaps the stellar binaries we predict to
form merging BBHs above the gap do not exist in nature. Another big
uncertainty concerns the location of the PISN mass gap which depends strongly
on the value of the ${}^{12}\text{C}(\alpha,\gamma)^{16}\text{O}$ reaction
rate (Takahashi, 2018; Farmer et al., 2019). Recent calculations using updated
values for this rate (Mehta et al., 2022; Farag et al., 2022) have pushed the
mass gap to the range $\sim 60M_{\odot}-140M_{\odot}$ in BH masses, shifting
both boundaries upwards by $\sim 20M_{\odot}$. This can further lower
predicted rates for BBHs above the PISN gap with a corresponding increase in
the amount of detections below it. In any case, we expect the spins of BHs
formed above the gap to be significantly restricted by the limit of critical
rotation on their progenitors, with spins below $0.25$, providing a clear
distinction from BHs formed through hierarchical mergers.

###### Acknowledgements.

PM acknowledges support from the FWO junior fellowship number 12ZY520N and the
FWO senior fellowship number 12ZY523N. IM is a recipient of the Australian
Research Council Future Fellowship FT190100574. Analysis of results relied on
the Julia programming language (Bezanson et al., 2017) and the Makie plotting
package (Danisch & Krumbiegel, 2021).

## References

  * Abbott et al. (2016) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Physical Review Letters, 116, 061102
  * Abbott et al. (2020a) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2020a, Living Reviews in Relativity, 23, 3
  * Abbott et al. (2019) Abbott, B. P., Abbott, R., Abbott, T. D., Abraham, S., & et al. 2019, Physical Review X, 9, 031040
  * Abbott et al. (2020b) Abbott, R., Abbott, T. D., Abraham, S., & et al. 2020b, arXiv e-prints, arXiv:2010.14527
  * Abbott et al. (2021) Abbott, R., Abbott, T. D., Acernese, F., & et al. 2021, arXiv e-prints, arXiv:2111.03606
  * Abbott et al. (2023) Abbott, R., Abbott, T. D., Acernese, F., & et al. 2023, Physical Review X, 13, 011048
  * Adamcewicz et al. (2023) Adamcewicz, C., Lasky, P. D., & Thrane, E. 2023, arXiv e-prints, arXiv:2307.15278
  * Aerts et al. (2019) Aerts, C., Mathis, S., & Rogers, T. M. 2019, ARA&A, 57, 35
  * Ajith et al. (2011) Ajith, P., Hannam, M., Husa, S., et al. 2011, Phys. Rev. Lett., 106, 241101
  * Angulo et al. (1999) Angulo, C., Arnould, M., Rayet, M., et al. 1999, Nuclear Physics A, 656, 3
  * Antonelli et al. (2023) Antonelli, A., Kritos, K., Ng, K. K. Y., Cotesta, R., & Berti, E. 2023, arXiv e-prints, arXiv:2306.11088
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481
  * Batta & Ramirez-Ruiz (2019) Batta, A. & Ramirez-Ruiz, E. 2019, arXiv e-prints, arXiv:1904.04835
  * Bavera et al. (2020) Bavera, S. S., Fragos, T., Qin, Y., et al. 2020, A&A, 635, A97
  * Bavera et al. (2022) Bavera, S. S., Fragos, T., Zapartas, E., et al. 2022, A&A, 657, L8
  * Bavera et al. (2021) Bavera, S. S., Fragos, T., Zevin, M., et al. 2021, A&A, 647, A153
  * Belczynski et al. (2014) Belczynski, K., Buonanno, A., Cantiello, M., et al. 2014, ApJ, 789, 120
  * Belczynski et al. (2016) Belczynski, K., Heger, A., Gladysz, W., et al. 2016, A&A, 594, A97
  * Belczynski et al. (2020) Belczynski, K., Klencki, J., Fields, C. E., et al. 2020, A&A, 636, A104
  * Belczynski et al. (2022) Belczynski, K., Romagnolo, A., Olejak, A., et al. 2022, ApJ, 925, 69
  * Berti & Volonteri (2008) Berti, E. & Volonteri, M. 2008, ApJ, 684, 822
  * Bezanson et al. (2017) Bezanson, J., Edelman, A., Karpinski, S., & Shah, V. B. 2017, SIAM Review, 59, 65
  * Biscoveanu et al. (2022) Biscoveanu, S., Callister, T. A., Haster, C.-J., et al. 2022, ApJ, 932, L19
  * Blandford & Znajek (1977) Blandford, R. D. & Znajek, R. L. 1977, MNRAS, 179, 433
  * Bouffanais et al. (2019) Bouffanais, Y., Mapelli, M., Gerosa, D., et al. 2019, ApJ, 886, 25
  * Brookshaw & Tavani (1993) Brookshaw, L. & Tavani, M. 1993, ApJ, 410, 719
  * Callister et al. (2021) Callister, T. A., Haster, C.-J., Ng, K. K. Y., Vitale, S., & Farr, W. M. 2021, ApJ, 922, L5
  * Chugunov et al. (2007) Chugunov, A. I., Dewitt, H. E., & Yakovlev, D. G. 2007, Phys. Rev. D, 76, 025028
  * Cyburt et al. (2010) Cyburt, R. H., Amthor, A. M., Ferguson, R., et al. 2010, ApJS, 189, 240
  * Danisch & Krumbiegel (2021) Danisch, S. & Krumbiegel, J. 2021, Journal of Open Source Software, 6, 3349
  * de Mink et al. (2009) de Mink, S. E., Cantiello, M., Langer, N., et al. 2009, A&A, 497, 243
  * de Mink & Mandel (2016) de Mink, S. E. & Mandel, I. 2016, MNRAS, 460, 3545
  * du Buisson et al. (2020) du Buisson, L., Marchant, P., Podsiadlowski, P., et al. 2020, MNRAS, 499, 5941
  * Eggleton (1983) Eggleton, P. P. 1983, ApJ, 268, 368
  * Eker et al. (2018) Eker, Z., Bakış, V., Bilir, S., et al. 2018, MNRAS, 479, 5491
  * Fabry et al. (2022) Fabry, M., Marchant, P., & Sana, H. 2022, A&A, 661, A123
  * Farag et al. (2022) Farag, E., Renzo, M., Farmer, R., Chidester, M. T., & Timmes, F. X. 2022, ApJ, 937, 112
  * Farmer et al. (2019) Farmer, R., Renzo, M., de Mink, S. E., Marchant, P., & Justham, S. 2019, ApJ, 887, 53
  * Ferguson et al. (2005) Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585
  * Fishbach & Holz (2017) Fishbach, M. & Holz, D. E. 2017, ApJ, 851, L25
  * Franciolini & Pani (2022) Franciolini, G. & Pani, P. 2022, Phys. Rev. D, 105, 123024
  * Fryer (1999) Fryer, C. L. 1999, ApJ, 522, 413
  * Fuller & Lu (2022) Fuller, J. & Lu, W. 2022, MNRAS, 511, 3951
  * Fuller & Ma (2019) Fuller, J. & Ma, L. 2019, ApJ, 881, L1
  * Fuller et al. (2019) Fuller, J., Piro, A. L., & Jermyn, A. S. 2019, MNRAS, 485, 3661
  * Galama et al. (1998) Galama, T. J., Vreeswijk, P. M., van Paradijs, J., et al. 1998, Nature, 395, 670
  * Galaudage et al. (2021) Galaudage, S., Talbot, C., Nagar, T., et al. 2021, ApJ, 921, L15
  * Gerosa & Berti (2017) Gerosa, D. & Berti, E. 2017, Phys. Rev. D, 95, 124046
  * Gerosa & Fishbach (2021) Gerosa, D. & Fishbach, M. 2021, Nature Astronomy, 5, 749
  * Grevesse et al. (1996) Grevesse, N., Noels, A., & Sauval, A. J. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 99, Cosmic Abundances, ed. S. S. Holt & G. Sonneborn, 117
  * Hall & Evans (2019) Hall, E. D. & Evans, M. 2019, Classical and Quantum Gravity, 36, 225002
  * Hannam et al. (2013) Hannam, M., Brown, D. A., Fairhurst, S., Fryer, C. L., & Harry, I. W. 2013, ApJ, 766, L14
  * Heger & Woosley (2002) Heger, A. & Woosley, S. E. 2002, ApJ, 567, 532
  * Hoy et al. (2022) Hoy, C., Mills, C., & Fairhurst, S. 2022, Phys. Rev. D, 106, 023019
  * Iglesias & Rogers (1996) Iglesias, C. A. & Rogers, F. J. 1996, ApJ, 464, 943
  * Irwin (2012) Irwin, A. W. 2012, FreeEOS: Equation of State for stellar interiors calculations
  * Itoh et al. (1996) Itoh, N., Hayashi, H., Nishikawa, A., & Kohyama, Y. 1996, ApJS, 102, 411
  * Izzard et al. (2004) Izzard, R. G., Ramirez-Ruiz, E., & Tout, C. A. 2004, MNRAS, 348, 1215
  * Jermyn et al. (2023) Jermyn, A. S., Bauer, E. B., Schwab, J., et al. 2023, ApJS, 265, 15
  * Jermyn et al. (2021) Jermyn, A. S., Schwab, J., Bauer, E., Timmes, F. X., & Potekhin, A. Y. 2021, ApJ, 913, 72
  * Kippenhahn et al. (2013) Kippenhahn, R., Weigert, A., & Weiss, A. 2013, Stellar Structure and Evolution
  * Köhler et al. (2015) Köhler, K., Langer, N., de Koter, A., et al. 2015, A&A, 573, A71
  * Langer (1997) Langer, N. 1997, in Astronomical Society of the Pacific Conference Series, Vol. 120, Luminous Blue Variables: Massive Stars in Transition, ed. A. Nota & H. Lamers, 83
  * Ma & Fuller (2023) Ma, L. & Fuller, J. 2023, arXiv e-prints, arXiv:2305.08356
  * MacFadyen & Woosley (1999) MacFadyen, A. I. & Woosley, S. E. 1999, ApJ, 524, 262
  * MacLeod & Loeb (2020) MacLeod, M. & Loeb, A. 2020, ApJ, 902, 85
  * Maeder (1987) Maeder, A. 1987, A&A, 178, 159
  * Maeder & Meynet (2000) Maeder, A. & Meynet, G. 2000, A&A, 361, 159
  * Mandel & de Mink (2016) Mandel, I. & de Mink, S. E. 2016, MNRAS, 458, 2634
  * Mandel & O’Shaughnessy (2010) Mandel, I. & O’Shaughnessy, R. 2010, Classical and Quantum Gravity, 27, 114007
  * Mandel & Smith (2021) Mandel, I. & Smith, R. J. E. 2021, ApJ, 922, L14
  * Marchant et al. (2016) Marchant, P., Langer, N., Podsiadlowski, P., Tauris, T. M., & Moriya, T. J. 2016, A&A, 588, A50
  * Marchant et al. (2019) Marchant, P., Renzo, M., Farmer, R., et al. 2019, ApJ, 882, 36
  * McKinney (2005) McKinney, J. C. 2005, ApJ, 630, L5
  * Mehta et al. (2022) Mehta, A. K., Buonanno, A., Gair, J., et al. 2022, ApJ, 924, 39
  * Miller-Jones et al. (2021) Miller-Jones, J. C. A., Bahramian, A., Orosz, J. A., et al. 2021, Science, 371, 1046
  * Mirabel & Rodrigues (2003) Mirabel, I. F. & Rodrigues, I. 2003, Science, 300, 1119
  * Murguia-Berthier et al. (2020) Murguia-Berthier, A., Batta, A., Janiuk, A., et al. 2020, ApJ, 901, L24
  * Nitz et al. (2023) Nitz, A. H., Kumar, S., Wang, Y.-F., et al. 2023, ApJ, 946, 59
  * Paxton et al. (2011) Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3
  * Paxton et al. (2013) Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4
  * Paxton et al. (2015) Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15
  * Paxton et al. (2018) Paxton, B., Schwab, J., Bauer, E. B., et al. 2018, ApJS, 234, 34
  * Paxton et al. (2019) Paxton, B., Smolec, R., Schwab, J., et al. 2019, ApJS, 243, 10
  * Périgois et al. (2023) Périgois, C., Mapelli, M., Santoliquido, F., Bouffanais, Y., & Rufolo, R. 2023, arXiv e-prints, arXiv:2301.01312
  * Peters (1964) Peters, P. C. 1964, Physical Review, 136, 1224
  * Podsiadlowski et al. (2004) Podsiadlowski, P., Mazzali, P. A., Nomoto, K., Lazzati, D., & Cappellaro, E. 2004, ApJ, 607, L17
  * Poniatowski et al. (2021) Poniatowski, L. G., Sundqvist, J. O., Kee, N. D., et al. 2021, A&A, 647, A151
  * Potekhin & Chabrier (2010) Potekhin, A. Y. & Chabrier, G. 2010, Contributions to Plasma Physics, 50, 82
  * Poutanen (2017) Poutanen, J. 2017, ApJ, 835, 119
  * Pretorius (2005) Pretorius, F. 2005, Phys. Rev. Lett., 95, 121101
  * Qin et al. (2018) Qin, Y., Fragos, T., Meynet, G., et al. 2018, A&A, 616, A28
  * Qin et al. (2022) Qin, Y., Wang, Y.-Z., Bavera, S. S., et al. 2022, ApJ, 941, 179
  * Qiu Cheng et al. (2023) Qiu Cheng, A., Zevin, M., & Vitale, S. 2023, arXiv e-prints, arXiv:2307.03129
  * Riley et al. (2021) Riley, J., Mandel, I., Marchant, P., et al. 2021, MNRAS, 505, 663
  * Rivinius et al. (2013) Rivinius, T., Carciofi, A. C., & Martayan, C. 2013, A&A Rev., 21, 69
  * Rodriguez et al. (2019) Rodriguez, C. L., Zevin, M., Amaro-Seoane, P., et al. 2019, Phys. Rev. D, 100, 043027
  * Rodriguez et al. (2016) Rodriguez, C. L., Zevin, M., Pankow, C., Kalogera, V., & Rasio, F. A. 2016, ApJ, 832, L2
  * Rogers & Nayfonov (2002) Rogers, F. J. & Nayfonov, A. 2002, ApJ, 576, 1064
  * Roulet et al. (2021) Roulet, J., Chia, H. S., Olsen, S., et al. 2021, Phys. Rev. D, 104, 083010
  * Safarzadeh et al. (2020) Safarzadeh, M., Farr, W. M., & Ramirez-Ruiz, E. 2020, ApJ, 894, 129
  * Sanyal et al. (2015) Sanyal, D., Grassitelli, L., Langer, N., & Bestenlehner, J. M. 2015, A&A, 580, A20
  * Saumon et al. (1995) Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713
  * Schrøder et al. (2021) Schrøder, S. L., MacLeod, M., Ramirez-Ruiz, E., et al. 2021, arXiv e-prints, arXiv:2107.09675
  * Shenar et al. (2022) Shenar, T., Sana, H., Mahy, L., et al. 2022, Nature Astronomy, 6, 1085
  * Song et al. (2016) Song, H. F., Meynet, G., Maeder, A., Ekström, S., & Eggenberger, P. 2016, A&A, 585, A120
  * Spruit (2002) Spruit, H. C. 2002, A&A, 381, 923
  * Stevenson & Clarke (2022) Stevenson, S. & Clarke, T. A. 2022, MNRAS, 517, 4034
  * Suijs et al. (2008) Suijs, M. P. L., Langer, N., Poelarends, A. J., et al. 2008, A&A, 481, L87
  * Tagawa et al. (2021) Tagawa, H., Haiman, Z., Bartos, I., Kocsis, B., & Omukai, K. 2021, MNRAS, 507, 3362
  * Takahashi (2018) Takahashi, K. 2018, ApJ, 863, 153
  * Talbot & Thrane (2018) Talbot, C. & Thrane, E. 2018, ApJ, 856, 173
  * Taylor & Kobayashi (2015) Taylor, P. & Kobayashi, C. 2015, MNRAS, 452, L59
  * Timmes et al. (2000) Timmes, F. X., Hoffman, R. D., & Woosley, S. E. 2000, ApJS, 129, 377
  * Torres et al. (2010) Torres, G., Andersen, J., & Giménez, A. 2010, A&A Rev., 18, 67
  * van Son et al. (2022) van Son, L. A. C., de Mink, S. E., Renzo, M., et al. 2022, ApJ, 940, 184
  * Wong et al. (2021) Wong, K. W. K., Breivik, K., Kremer, K., & Callister, T. 2021, Phys. Rev. D, 103, 083021
  * Woosley (1993) Woosley, S. E. 1993, ApJ, 405, 273
  * Woosley (2017) Woosley, S. E. 2017, ApJ, 836, 244
  * Zahn (1977) Zahn, J.-P. 1977, A&A, 57, 383
  * Zaldarriaga et al. (2018) Zaldarriaga, M., Kushnir, D., & Kollmeier, J. A. 2018, MNRAS, 473, 4174
  * Zevin et al. (2021) Zevin, M., Bavera, S. S., Berry, C. P. L., et al. 2021, ApJ, 910, 152
  * Zevin et al. (2017) Zevin, M., Pankow, C., Rodriguez, C. L., et al. 2017, ApJ, 846, 82

##  Appendix A Details on microphysical inputs of MESA simulations

Our MESA calculations rely on various inputs regarding microphysical
processes, which we list here. Opacities for various temperature and density
regimes are taken from OPAL (Iglesias & Rogers 1996), low-temperature
opacities from Ferguson et al. (2005) and Compton-scattering opacities from
Poutanen (2017). The equation of state is determined from a combination of
Rogers & Nayfonov (2002), Irwin (2012), Saumon et al. (1995), Timmes et al.
(2000), Potekhin & Chabrier (2010) and Jermyn et al. (2021), with different
equations of state used at different conditions for temperature, density and
composition (see Jermyn et al. 2023 for details on how these are patched
together). Nuclear reaction rates are taken from Angulo et al. (1999) and
Cyburt et al. (2010), with screening accounted for as in Chugunov et al.
(2007). Thermal neutrino emission is determined as in Itoh et al. (1996).

##  Appendix B Fit to maximum BH spin

The values of the maximum dimensionless spin derived from our simulation are
given in Table 2. These correspond to the value of $a_{\mathrm{max}}$ at
carbon depletion in Figure 2. For convenience, we provide the following fit to
our results:

$\displaystyle a_{\mathrm{max}}(x)=10^{c_{0}+c_{1}x+c_{2}x^{2}},$ (24)

where $x=\log_{10}(M_{\mathrm{BH}}/M_{\odot})$. The dependence on metallicity
is included through $c_{0}$, $c_{1}$ and $c_{2}$,

$\displaystyle c_{0}$ $\displaystyle=$ $\displaystyle 0.13\sqrt{z}$ (25)
$\displaystyle c_{1}$ $\displaystyle=$ $\displaystyle 0.74-0.19\sqrt{z}$ (26)
$\displaystyle c_{2}$ $\displaystyle=$ $\displaystyle-0.58+0.12\sqrt{z},$ (27)

with $z=\log_{10}(0.1Z_{\odot}/Z)$. This fit has not been tested beyond the
range of masses and metallicities explored in Section 2. It is also important
to point out that these predictions rely on an understanding of the radii of
very massive helium stars near the Eddington limit. Our results have been
computed using hydrostatic stellar evolution models, but an accurate
prediction of the stellar radius in this regime requires the inclusion of
hydrodynamics (e.g. Poniatowski et al. 2021).

Table 2: Maximum black hole dimensionless spin $a_{\mathrm{max}}$ as a function of black hole mass and progenitor metallicity. These results correspond to the values at carbon depletion of Figure 2. | $a_{\mathrm{max}}$
---|---
$\log_{10}\left(M_{\mathrm{BH}}/M_{\odot}\right)$ | $Z_{\odot}/10$ | $Z_{\odot}/50$ | $Z_{\odot}/250$
0.8 | 1.74 | 1.97 | 2.02
1.0 | 1.35 | 1.55 | 1.60
1.2 | 1.03 | 1.20 | 1.25
1.4 | 0.753 | 0.905 | 0.945
1.6 | 0.517 | 0.649 | 0.684
1.8 | 0.301 | 0.423 | 0.453
2.0 | 0.139 | 0.240 | 0.269
2.2 | 0.0583 | 0.121 | 0.148
2.4 | 0.0271 | 0.0514 | 0.0619