# Black Hole Mergers from Star Clusters with Top-Heavy Initial Mass Functions

Newlin C. Weatherford,11affiliationmark:  22affiliationmark:  Giacomo
Fragione,11affiliationmark:  22affiliationmark:  Kyle
Kremer,33affiliationmark:  44affiliationmark:  Sourav
Chatterjee,55affiliationmark:  Claire S. Ye,11affiliationmark:
22affiliationmark:  Carl L. Rodriguez,66affiliationmark:  Frederic A.
Rasio11affiliationmark:  22affiliationmark:  1Center for Interdisciplinary
Exploration & Research in Astrophysics (CIERA), Evanston, IL 60202, USA
2Department of Physics & Astronomy, Northwestern University, Evanston, IL
60202, USA  3TAPIR, California Institute of Technology, Pasadena, CA 91125,
USA  4The Observatories of the Carnegie Institution for Science, Pasadena, CA
91101, USA  5Tata Institute of Fundamental Research, Homi Bhabha Road, Mumbai
400005, India  6McWilliams Center for Cosmology, Department of Physics,
Carnegie Mellon University, Pittsburgh, PA 15213, USA

###### Abstract

Recent observations of globular clusters (GCs) provide evidence that the
stellar initial mass function (IMF) may not be universal, suggesting
specifically that the IMF grows increasingly top-heavy with decreasing
metallicity and increasing gas density. Non-canonical IMFs can greatly affect
the evolution of GCs, mainly because the high end determines how many black
holes (BHs) form. Here we compute a new set of GC models, varying the IMF
within observational uncertainties. We find that GCs with top-heavy IMFs lose
most of their mass within a few Gyr through stellar winds and tidal stripping.
Heating of the cluster through BH mass segregation greatly enhances this
process. We show that, as they approach complete dissolution, GCs with top-
heavy IMFs can evolve into ‘dark clusters’ consisting of mostly BHs by mass.
In addition to producing more BHs, GCs with top-heavy IMFs also produce many
more binary BH (BBH) mergers. Even though these clusters are short-lived,
mergers of ejected BBHs continue at a rate comparable to, or greater than,
what is found for long-lived GCs with canonical IMFs. Therefore these
clusters, although they are no longer visible today, could still contribute
significantly to the local BBH merger rate detectable by LIGO/Virgo,
especially for sources with higher component masses well into the BH mass gap.
We also report that one of our GC models with a top-heavy IMF produces dozens
of intermediate-mass black holes (IMBHs) with masses $M>100\,{\rm M_{\odot}}$,
including one with $M>500\,{\rm M_{\odot}}$. Ultimately, additional
gravitational wave observations will provide strong constraints on the stellar
IMF in old GCs and the formation of IMBHs at high redshift.

##  1\. Introduction

Among the densest environments in the Universe, globular clusters (GCs) are
ideal laboratories to investigate the importance of stellar dynamics in
forming and evolving compact objects and compact binaries (see, e.g., Heggie &
Hut, 2003). Frequent dynamical encounters between cluster members are
fundamental in creating and explaining the existence of a number of exotic
populations, such as X-ray binaries (e.g., Clark, 1975; Verbunt et al., 1984;
Giesler et al., 2018; Kremer et al., 2018), radio pulsars (e.g., Lyne et al.,
1987; Sigurdsson & Phinney, 1995; Ivanova et al., 2008; Ye et al., 2019), and
gravitational wave (GW) sources (e.g., Rodriguez et al., 2015; Askar et al.,
2017; Banerjee, 2017; Fragione & Kocsis, 2018; Kremer et al., 2019).

Many results of GC modeling rely on the assumption that the stellar initial
mass function (IMF) has the form of a canonical Kroupa (2001) IMF.
Observations, however, suggest the IMF may not be universal. For example,
Milky Way GCs with low central densities appear deficient in low-mass stars
(De Marchi et al., 2007), while GCs in the Andromeda galaxy exhibit a trend
between metallicity and mass-to-light ratio that only a non-canonical, top-
heavy IMF could explain (Haghi et al., 2017). Ultra-compact dwarf galaxies
have large dynamical mass-to-light ratios and appear to contain an
overabundance of LMXB sources (Dabringhausen et al., 2009). Our Galactic
center provides an extreme example of a non-canonical IMF; stars formed there
in the past few Myr have masses consistent with a top-heavy mass function
(Bartko et al., 2010). All these observations can be explained if the stellar
IMF becomes increasingly top-heavy with decreasing metallicity and increasing
gas density of the pre-GC cloud (Marks et al., 2012). This is theoretically
expected given the Jeans mass instability in molecular clouds (Larson, 1998)
and self-regulation of accretion onto forming stars (Adams & Fatuzzo, 1996;
Adams & Laughlin, 1996). Metallicity plays a decisive role by regulating line-
emission cooling in the collapsing gas cloud and radiation pressure against
stellar accretion (De Marchi et al., 2017).

The IMF has long been known to strongly impact the dynamical evolution and
survival of GCs. Already Chernoff & Weinberg (1990) found enhanced mass loss
rates in cluster models with flatter IMFs, due to increased winds from high-
mass stars and faster halo expansion and evaporation. More recent results show
that GC models with top-heavy IMFs dissolve much faster than models with
canonical IMFs (Chatterjee et al., 2017; Giersz et al., 2019). In particular,
clusters with top-heavy IMFs produce more numerous and more massive black
holes (BHs). Crucially, this promotes the BH-burning process, in which strong
dynamical encounters with BHs provide energy to stellar populations, inflating
the cluster halo (Mackey et al., 2007, 2008; Breen & Heggie, 2013; Kremer et
al., 2020a). Taken to an extreme by a top-heavy IMF, this mechanism will force
rapid and unstable evaporation through the tidal boundary and early cluster
dissolution (Giersz et al., 2019).

The impact on cluster evolution of varying BH abundance has recently been
further analyzed via a combination of analytical calculations and $N$-body
simulations (Breen & Heggie, 2013; Wang, 2020). However, these studies only
consider idealized star clusters with two components (stars and BHs) and no
stellar evolution, which plays a crucial role in cluster disruption (e.g.,
Chernoff & Weinberg, 1990). Apart from the few models with non-canonical IMFs
examined by Chatterjee et al. (2017) and Giersz et al. (2019), there is no
systematic study of the IMF’s influence on BH populations generated via
stellar evolution. In this Letter, we extend the grid of cluster models in the
CMC Cluster Catalog (Kremer et al., 2020b), focusing on the role of a varying
IMF. In particular, we study how non-canonical IMFs shape the evolution of
GCs, their BH populations, and the number of dynamically produced BBH mergers.

The Letter is organized as follows. In Section 2, we describe the parameters
of the numerical models we evolve. In Section 3, we analyze the structural
evolution of clusters and describe the properties of their black hole
populations, including BBH mergers. Finally, in Section 4, we discuss the
implications of our findings and lay out our conclusions.

##  2\. Cluster models

To evolve our cluster models, we use CMC (Cluster Monte Carlo), a Hénon-type
Monte Carlo code that includes the newest prescriptions for wind-driven mass
loss, compact object formation, and pulsational-pair instabilities (see Hénon,
1971a, b; Joshi et al., 2000, 2001; Fregeau et al., 2003; Chatterjee et al.,
2010, 2013; Pattabiraman et al., 2013; Rodriguez et al., 2015; Kremer et al.,
2020b, and references therein).

We generate $15$ independent cluster simulations varying the initial number of
particles (single stars plus binaries) $N_{\rm i}$ and initial total mass
$M_{\rm i}$, with uniform initial cluster virial radius ($r_{v}=1\,\rm{pc}$),
metallicity ($Z=0.1\,{\rm Z}_{\odot}$), and Galactocentric distance111Assuming
a Milky Way-like potential (e.g., Dehnen & Binney, 1998).
($R_{\rm{gc}}=8\,\rm{kpc}$). These choices of $r_{v}$, $R_{\rm gc}$, and $Z$
are known to result in models that closely match typical Milky Way GCs when
using a canonical IMF (Kremer et al., 2020b).

We assume that all models are initially described by a King profile with
concentration $W_{0}=5$ (King, 1962). Primary stellar masses are sampled from
the Kroupa (2001) multi-component power-law IMF,

$\xi(m)\propto\begin{cases}{m}^{-1.3}&\text{$0.08\leq m/\mathrm{M}_{\odot}\leq
0.5$}\\\ {m}^{-2.3}&\text{$0.5\leq m/\mathrm{M}_{\odot}\leq 1.0$}\\\
{m}^{-\alpha_{3}}&\text{$1.0\leq m/\mathrm{M}_{\odot}\leq
150.0$}\,.\end{cases}$ (1)

We choose three different values for $\alpha_{3}=(1.6,2.3,3.0)$, corresponding
to the 95% confidence interval around the canonical value, $\alpha_{3}=2.3$
(Kroupa, 2001). We fix the primordial stellar binary fraction to $f_{b}=5\%$
and draw secondary masses from a uniform distribution in mass ratio (e.g.,
Duquennoy & Mayor, 1991). Binary orbital periods are sampled from a log-
uniform distribution, with orbital separations ranging from near contact to
the hard/soft boundary, while binary eccentricities are thermal (e.g., Heggie,
1975). For details on stellar evolution in CMC, see Kremer et al. (2020b). We
compute GW recoil kicks for BH merger products following the methods described
in Rodriguez et al. (2019) and references therein. We assume all BHs have zero
natal spin while BH merger products are assigned new spins of $\sim 0.7$
(e.g., Berti et al., 2007; Tichy & Marronetti, 2008; Lousto et al., 2010),
which are then taken into account if they merge again.

Table 1Model initial conditions and data on cumulative BH formation and BBH mergers from $0\leq t\leq 13\,{\,\rm Gyr}$; upper-IMF slope $\alpha_{3}$, initial number of particles $N_{i}$, initial mass $M_{i}$, total BHs formed ${N_{\rm BH}}$, total BBH mergers $N_{\rm BBH}^{\rm(merg)}$, ejected BBH mergers $N_{\rm BBH,ejected}^{\rm(merg)}$, BBH mergers with a ‘mass-gap’ BH component ($M_{\rm BH}>40.5\,{\rm M_{\odot}}$) formed via stellar collisions $N_{\rm BBH,gap,coll}^{\rm(merg)}$, multiple-generation BBH mergers with a mass-gap BH component $N_{\rm BBH,gap,2G+}^{\rm(merg)}$, and total BBH mergers at redshifts $z<1$, $N_{\rm BBH}^{\rm(merg)}(z<1)$, assuming that all GCs were born $13\,{\,\rm Gyr}$ ago. Note that the BBH merger totals in the last 5 columns for models 1-6 are all lower limits since we do not consider mergers resulting from dynamics in ‘dark clusters’ (see Section 3.3). |  |  |  |  | BBH mergers
---|---|---|---|---|---
| $\alpha_{3}$ |  $N_{\rm i}$ ($10^{5}$) |  $M_{\rm i}$ ($10^{5}\,\,\rm M_{\odot}$) | $N_{\rm BH}$ | $N_{\rm BBH}^{\rm(merg)}$ | $N_{\rm BBH,ejected}^{\rm(merg)}$ | $N_{\rm BBH,gap,coll}^{\rm(merg)}$ | $N_{\rm BBH,gap,2G+}^{\rm(merg)}$ | $N_{\rm BBH}^{\rm(merg)}(z<1)$
1 | 1.6 | 4 | 8.5 | 11263 | 95 | 58 | 1 | 7 | 7
2 | 1.6 | 8 | 17.0 | 22539 | 222 | 132 | 7 | 26 | 19
3 | 1.6 | 16 | 34.2 | 45358 | 434 | 195 | 32 | 83 | 14
4 | 1.6 | 1.1 | 2.4 | 3224 | 20 | 15 | 0 | 2 | 0
5 | 1.6 | 2.3 | 4.8 | 6319 | 57 | 42 | 1 | 3 | 5
6 | 1.6 | 4.5 | 9.6 | 12764 | 111 | 73 | 5 | 10 | 11
7 | 2.3 | 4 | 2.4 | 1114 | 38 | 17 | 1 | 3 | 5
8 | 2.3 | 8 | 4.8 | 2252 | 91 | 36 | 0 | 9 | 13
9 | 2.3 | 16 | 9.7 | 4518 | 233 | 106 | 1 | 23 | 42
10 | 3.0 | 4 | 1.6 | 113 | 6 | 0 | 0 | 1 | 0
11 | 3.0 | 8 | 3.3 | 247 | 9 | 3 | 0 | 1 | 2
12 | 3.0 | 16 | 6.6 | 499 | 32 | 6 | 1 | 2 | 1
13 | 3.0 | 5.9 | 2.4 | 188 | 9 | 1 | 0 | 1 | 0
14 | 3.0 | 11.8 | 4.8 | 341 | 28 | 3 | 1 | 1 | 2
15 | 3.0 | 23.6 | 9.7 | 744 | 71 | 9 | 2 | 5 | 7

Table 1 summarizes the model parameters considered in our work. Each
simulation evolves to a final time $T_{\rm H}\approx 14$ Gyr, unless the
cluster disrupts. In these cases, which occur only for our models with
$\alpha_{3}=1.6$, simulation terminates once there are fewer than 300
particles per CPU, the default particle count used to compute local densities
throughout the cluster. This typically corresponds to a few$\times 10^{4}$
stars in the GC, beyond the point where the assumptions of the Monte Carlo
method – namely spherical symmetry and a relaxation timescale longer than the
dynamical timescale – start to break down (Chatterjee et al., 2017). However,
the exact cut-off point in the simulation does not affect the general
evolutionary trends of the dissolution process. Further accurate evolution of
the remaining ‘dark star cluster’ would require a switch to direct $N$-body
methods (e.g., Banerjee & Kroupa, 2011).

##  3\. Results

In this Section, we discuss how different IMFs affect the evolution of the
cluster mass, size, and BH population. We do not find significant qualitative
differences between the cases where we keep the initial number of particles
constant (models 1-3, 7-9, and 10-12 in Table 1) or the initial total mass
constant (models 4-9 and 13-15).

Figure 1.— Evolution of the total mass (top-left), half-mass radius $r_{h}$
(bottom-left), number of BHs $N_{\rm BH}$ (top-right), and fractional number
of BHs $N_{\rm BH}/N$ (bottom-right) for cluster models with different values
of $\alpha_{3}$ (blue $1.6$, black $2.3$, red $3.0$). Different line styles
represent different initial number of particles (dotted $4\times 10^{5}$,
dashed $8\times 10^{5}$, solid $1.6\times 10^{6}$). The simulations terminate
at $T_{H}=14\,{\,\rm Gyr}$, except for $\alpha_{3}=1.6$, where we cut off the
plot just before $r_{h}$ starts to drop sharply (around this point,
assumptions in our Monte Carlo methods stop being valid; see text). In the
top-left panel, the four mass loss phases discussed in the text are
highlighted via shaded intervals and the dominant mass loss mechanisms
indicated. The bifurcation of the third phase is due to decreased duration of
this phase for $\alpha_{3}=1.6$.

###  3.1. Cluster evolution

Variations in the high end of the IMF can dramatically affect cluster
evolution and survival, especially by setting the number of BHs that are
formed and retained, determining the degree to which BH-burning regulates the
cluster’s energy reservoir (Chatterjee et al., 2017). In this process, BHs
quickly mass-segregate into the cluster core on a sub-${\,\rm Gyr}$ timescale,
forming a central BH population that undergoes frequent phases of core
collapse and re-expansion (Morscher et al., 2015). During these events, the
BHs mix with the rest of the cluster and provide energy to passing stars in
scattering interactions (Breen & Heggie, 2013). Cumulatively, these
interactions inflate the cluster halo and force faster evaporation through the
tidal boundary (Chatterjee et al., 2017; Giersz et al., 2019; Wang, 2020).

In Figure 1, we show the evolution of the total mass (top-left panel) for star
clusters with different $N_{i}$ and $\alpha_{3}$. Clusters evolve very
differently depending on the choice of $\alpha_{3}$. The models with
$\alpha_{3}=2.3$ and $3$ survive until $T_{\rm H}$, while the models with
$\alpha_{3}=1.6$ appear to disrupt at about $t\sim 1\,{\,\rm Gyr}$ with faster
disruption occurring for smaller $N_{i}$. 222While these models may not
survive with significant mass beyond a few ${\,\rm Gyr}$, models with larger
$N_{i}$, larger Galactocentric distance, and/or smaller virial radius could
feasibly survive much longer. The different fates of the cluster models are
understandable given the mechanisms that power cluster mass loss at different
epochs:

  1. 1.

Initially, mass is lost primarily via stellar winds from massive stars (e.g.,
Chernoff & Weinberg, 1990). This relatively quiet first phase is visible early
in the top-left panel of Figure 1, terminating around $t\approx 3\,{\,\rm
Myr}$.

  2. 2.

The second phase of mass loss, extending from $3\,{\,\rm Myr}\lesssim
t\lesssim 60\,{\,\rm Myr}$, is driven primarily by Type II supernovae (SNe)
accompanying compact object formation. This phase most clearly depends on the
choice of $\alpha_{3}$. For $\alpha_{3}=3$, SNe-driven mass loss from few
massive stars is negligible compared to ongoing stellar wind mass loss among
the far more numerous low-mass stars. For $\alpha_{3}=1.6$, however, the SNe-
driven mass loss phase is nearly catastrophic, causing the cluster to lose
close to half its original mass by the time compact object formation slows
($t\approx 60\,{\,\rm Myr}$). Ejections of lighter compact objects due to
natal kicks also contribute to mass loss in this phase starting at $t\approx
6\,{\,\rm Myr}$, causing a plateau in ${N_{\rm BH}}$ (top-right panel).

  3. 3.

A ‘third phase’ of mass loss begins at the end of rapid compact object
formation. Really just a return of the first phase, this period of mass loss
is once again dominated by stellar winds, this time from lower-mass stars. For
the models where $\alpha_{3}=2.3$ or $3$, stellar winds continue to account
for most of the mass loss until about $6\,{\,\rm Gyr}$. The wind-dominated
mass loss phase ends much earlier for the $\alpha_{3}$ = 1.6 models – around
$t\approx 200\,{\,\rm Myr}$.

  4. 4.

The fourth and final phase of mass loss is driven by evaporation through the
tidal boundary. This phase depends strongly on the IMF due to its influence on
the number of BHs formed, and thereby the degree to which BH-burning inflates
the halo. In the case of extreme BH-burning, the cluster quickly expands to
fill its Roche lobe within the assumed Galactic potential and starts to
tidally disrupt, as seen in the models with $\alpha_{3}=1.6$. Note that
similar accelerating mass loss is also seen at late times in the models with
$\alpha_{3}=3$. In these models, however, the runaway process of halo
evaporation and tidal contraction is not stimulated by extreme BH-burning.
Instead, these top-light, low-mass models simply start out with much smaller
tidal radii. Hence they see significant evaporation rates relatively quickly
compared to canonical models ($\alpha_{3}=2.3$), despite forming far fewer
BHs.

The evolution of the half-mass radius $r_{h}$ (bottom-left panel, Figure 1)
also reflects this general picture. In all models, $r_{h}$ expands as a
consequence of stellar mass loss. Enhanced BH-burning in models with
$\alpha_{3}=1.6$ drives faster cluster expansion such that these models
overflow their Roche lobes after a few hundred ${\,\rm Myr}$. Due to more
energetic BH-burning, these clusters lose about $80\%$ of their mass by $t\sim
1$ Gyr. Models with $\alpha_{3}=3$ do not exhibit rapid expansion, evolving
more gradually as a result of reduced massive star formation and the lack of a
significant central BH population.

Figure 2.— Evolution of the number of binary black holes $N_{\rm BBH}$ bound
to the star clusters with different values of $\alpha_{3}$ (blue $1.6$, black
$2.3$, red $3.0$) and different initial number of particles ($N_{i}$).
Different line styles represent different $N_{i}$, as in Figure 1. Since BBHs
continuously form and disrupt in central scattering interactions, we apply a
rolling average over every $10^{4}$ time-steps to smooth the curves at low
$N_{\rm BBH}$.

###  3.2. Black hole population

In the top-right panel of Figure 1, we show the evolution of the number of BHs
${N_{\rm BH}}$. Clusters with top-heavy IMFs produce more high-mass stars and
thereby more BHs. As described above, natal kicks start to eject newly-formed
BHs after $t\approx 6\,{\,\rm Myr}$, causing ${N_{\rm BH}}$ to plateau for
each model. At this point, we find that ${N_{\rm BH}}$ in models with
$\alpha_{3}=1.6$ is $\sim 10$ times higher than in models with
$\alpha_{3}=2.3$ and $\sim 100$ times higher than in models with
$\alpha_{3}=3$ (see also the cumulative number of BHs formed in column 4 of
Table 1). During this phase of roughly constant ${N_{\rm BH}}$, dynamical
friction causes BHs to segregate to the cluster center. This phase’s duration
depends on $\alpha_{3}$ since clusters with top-heavy IMFs have higher average
stellar mass and correspondingly longer segregation timescales. Once the BHs
have segregated to the core, they are gradually ejected via strong dynamical
encounters.

In the bottom-right panel of Figure 1, we also show the evolution of the
fractional number of BHs ${N_{\rm BH}}/N$. In models with $\alpha_{3}=2.3$ or
$3$, ${N_{\rm BH}}/N$ decreases in time as strong encounters preferentially
eject massive objects (e.g., BHs) from the cluster core. In models with
$\alpha_{3}=1.6$, however, ${N_{\rm BH}}/N$ increases in time due to earlier
tidal evaporation. By the moment these clusters retain only $\sim 20\%$ of
their original mass, ${N_{\rm BH}}/N\approx 3$–$4\%$, corresponding to a
significant fraction of cluster mass (up to $75\%$ by the time $r_{h}$ peaks).
Note that the number and mass fraction of BHs at late times depends on the
cluster dissolution process beyond this point, which our Monte Carlo methods
are not designed to address. It is nevertheless clear that clusters with top-
heavy IMFs can evolve through a ‘dark cluster’ stage during which their total
mass is dominated by BHs (e.g., Banerjee & Kroupa, 2011).

Regardless of age, individual BHs are only directly detectable if they reside
in binaries, either via detection of GWs from BBH mergers or through the
presence of non-BH (nBH) companions. While dense star clusters are expected to
be efficient factories of BBH and BH–nBH binaries, their numbers typically
remain small as their dynamical assembly competes with their disruption and
ejection (e.g., Downing et al., 2010, 2011; Morscher et al., 2015; Kremer et
al., 2018). However, the choice of IMF greatly impacts BH binary formation.

In Figure 2, we plot the evolution of the number of BBHs $N_{\rm BBH}$.
Clusters with top-heavy IMFs produce more BBHs via enhanced formation of high-
mass stars and BHs. $N_{\rm BBH}$ plateaus a little later than ${N_{\rm BH}}$
(top-right panel, Figure 1), a delay reflecting the dynamical assembly of
BBHs. An equilibrium between $N_{\rm BBH}$ and ${N_{\rm BH}}$ is apparent
during this plateau, with of order one BBH for every 100 bound BHs in the
cluster. Subsequently, $N_{\rm BBH}$ decreases more rapidly for clusters with
higher ${N_{\rm BH}}$, as evident by the slope past the plateau. This is
unsurprising since BH-burning, enhanced for clusters with higher ${N_{\rm
BH}}$, is characterized by periodic collapse and re-expansion of the central
BH population, a process that disrupts and ejects many BBHs (Morscher et al.,
2015; Chatterjee et al., 2017).

Figure 3.— Cumulative binary black hole mergers with respect to redshift (and
time) for star clusters with different values of $\alpha_{3}$ (blue $1.6$,
black $2.3$, red $3.0$) and different initial number of particles ($N_{i}$).
Different line styles represent different $N_{i}$, as in Figure 1. All model
GCs were assumed to be born $13\,{\,\rm Gyr}$ ago and appropriate redshifts
were then computed using Astropy’s (Astropy Collaboration et al., 2013)
cosmology calculator under the flat $\Lambda$CDM model with $H_{0}=69.6\,{\rm
km/s/Mpc}$ and $\Omega_{\rm matter}=0.286$ (see, e.g., Bennett et al., 2014).
Values for $\alpha_{3}=1.6$ are lower limits (see Section 3.3).

###  3.3. BBH merger rates

We report in Table 1, column 5, the cumulative numbers of BBH mergers $N_{\rm
BBH}^{\rm(merg)}$ in each model through $13\,{\,\rm Gyr}$. It is readily
apparent that each reduction in $\alpha_{3}$ generally increases $N_{\rm
BBH}^{\rm(merg)}$, with some nuance introduced depending on which values of
$\alpha_{3}$ are compared and whether $N_{i}$ or $M_{i}$ is held constant in
the comparison. Most notably, $N_{\rm BBH}^{\rm(merg)}$ in short-lived
clusters with $\alpha_{3}=1.6$ is comparable to or greater than the number in
their longer-lived cousins with $\alpha_{3}=3$. This is largely due to the
sheer number of merging BBHs they eject (Table 1, column 6).

Even though clusters with top-heavy IMFs rapidly dissolve, the large number of
merging BBHs they eject suggests they could contribute significantly to the
total BBH merger rate in the local Universe. Indeed, we show in Figure 3 that
the distribution of BBH mergers across redshift $z$ is roughly comparable
between the models with $\alpha_{3}=1.6$ and $2.3$, even for $z<1$. The number
of BBH mergers in each model at $z<1$ is listed in column 9 of Table 1: a
total of 40 across models 1-3 with $\alpha_{3}=1.6$ vs 60 across models 7-9
with $\alpha_{3}=2.3$. This comparison is quite conservative since the merger
counts for the models with $\alpha_{3}=1.6$ are lower limits. First, we only
consider BBHs ejected prior to the simulation’s end, ignoring the potentially
substantial contribution to the $z<1$ merger rate from the large number of
bound BHs left in the ‘dark phase’ of these clusters’ evolution. Many of these
BHs could merge after ejection as binaries. Second, to compute redshift we
assume the clusters were uniformly born $13\,{\,\rm Gyr}$ ago, ignoring
entirely the possibility that clusters with top-heavy IMFs may be born more
recently (El-Badry et al., 2019). It is therefore plausible that clusters born
with top-heavy IMFs could still contribute significantly to the local BBH
merger rate even if they are dissolved by the present day (Fragione & Kocsis,
2018).

Since the above comparisons of $z<1$ merger counts are model-dependent, it is
useful to more generally estimate the BBH merger rate expected in clusters
with top-heavy IMFs. We do so using the total number of BBH mergers from each
model (which are still lower limits in the case of $\alpha_{3}=1.6$). The rate
due to clusters with $\alpha_{3}=x$ at a given redshift $z$ can be
approximated as $\Gamma_{x}(z)\approx\langle N_{\rm
BBH}^{\rm(merg)}/M_{i}\rangle\rho_{\rm SF}(z)f_{\rm SF}f_{\rm x}$. Here,
$\langle N_{\rm BBH}^{\rm(merg)}/M_{i}\rangle$ is the mean number of BBH
mergers per initial cluster mass, $\rho_{\rm SF}$ is the cosmological density
of the star formation rate, $f_{\rm SF}$ is the fraction of the star formation
rate assumed to occur in star clusters, and $f_{x}$ is the fraction of
clusters born with $\alpha_{3}\approx x$. To mitigate uncertainties in the
latter three terms, we simply compute $\langle N_{\rm
BBH}^{\rm(merg)}/M_{i}\rangle$ for each value of $\alpha_{3}$ and express the
estimated merger rates from clusters with non-canonical IMFs as ratios with
respect to the better-studied rates from clusters with canonical IMFs.

To compute $\langle N_{\rm BBH}^{\rm(merg)}/M_{i}\rangle$ for each
$\alpha_{3}$, we extract the functional dependence of $N_{\rm
BBH}^{\rm(merg)}$ on $M_{i}$ from columns 3 and 5 of Table 1. Though the
number of BBH mergers per cluster scales roughly linearly with the present-day
cluster mass, at least in clusters with canonical IMFs (e.g., Rodriguez et
al., 2015; Kremer et al., 2020b), we find a more general power law of the form
$N_{\rm BBH}^{\rm(merg)}=a(M_{i})^{b}$ fits the data slightly better given
$\chi^{2}$ goodness-of-fit tests. Specifically,

$N_{\rm BBH}^{\rm(merg)}\begin{cases}\geq(9.5\pm
1.0)\Bigr{(}\frac{M_{i}}{10^{5}{\rm M_{\odot}}}\Bigr{)}^{1.09\pm
0.04}&\text{$\alpha_{3}=1.6$}\\\ =(11.8\pm 0.5)\Bigr{(}\frac{M_{i}}{10^{5}{\rm
M_{\odot}}}\Bigr{)}^{1.31\pm 0.02}&\text{$\alpha_{3}=2.3$}\\\ =(1.9\pm
0.7)\Bigr{(}\frac{M_{i}}{10^{5}{\rm M_{\odot}}}\Bigr{)}^{1.6\pm
0.2}&\text{$\alpha_{3}=3.0$}\,,\end{cases}$ (2)

where the $1\sigma$ uncertainties on the fit parameters are computed assuming
Poisson uncertainties on $N_{\rm BBH}^{\rm(merg)}$ and the $\geq$ symbol
indicates a lower limit for top-heavy IMFs.

Assuming the distribution of GC birth masses takes the form $dN_{\rm
cluster}/dM_{i}\propto M_{i}^{-2}$ (Lada & Lada, 2003), then

$\displaystyle\biggr{\langle}\frac{N_{\rm
BBH}^{\rm(merg)}}{M_{i}}\biggr{\rangle}$ $\displaystyle=\ \frac{a}{(10^{5}{\rm
M_{\odot}})^{b}}\frac{\int_{M_{L}}^{M_{H}}M_{i}^{b-3}dM_{i}}{\int_{M_{L}}^{M_{H}}M_{i}^{-2}dM_{i}}$
(3) $\displaystyle=\begin{cases}\frac{aM_{H}\ln(M_{H}/M_{L})}{(10^{5}{\rm
M_{\odot}})^{2}(M_{H}/M_{L}-1)}&\text{$b=2$}\\\
\frac{aM_{L}^{b-1}}{(2-b)(10^{5}{\rm
M_{\odot}})^{b}}\Bigr{[}\frac{1-(M_{L}/M_{H})^{2-b}}{1-M_{L}/M_{H}}\Bigr{]}&\text{$b\neq
2$}\,,\end{cases}$

where $M_{L}$ and $M_{H}$ are the lower and upper limits of the mass function
for clusters capable of producing BBH mergers. Both bounds are somewhat
arbitrary, but the lower bound has a greater impact on the merger rate
calculation since the integrands in Equation 3 scale inversely with $M_{i}$.
To avoid extrapolating our fit to cluster masses that rarely produce BBHs, we
set the lower bound for each $\alpha_{3}$ as the $M_{i}$ that produces an
average of two stars with $M>25\,{\rm M_{\odot}}$, given the assumed IMF. This
is roughly the minimum cluster mass needed to produce a single BBH merger
since most progenitors with $M>25\,{\rm M_{\odot}}$ will collapse to a BH.
Under this definition, $M_{L}/(100{\rm M_{\odot}})\approx(2,9,80)$ for
$\alpha_{3}=(1.6,2.3,3)$, respectively. Although the maximum cluster mass
could vary with $\alpha_{3}$ in principle, observations do not well-constrain
this value, so we naïvely set the upper bound for all $\alpha_{3}$ to be
$M_{H}=10^{7}{\rm M_{\odot}}$ (Harris et al., 2014). With these assumptions,
we find that $\langle N_{\rm
BBH}^{\rm(merg)}/M_{i}\rangle\approx(59,40,10)/(10^{6}\,{\rm M_{\odot}})$,
respectively. The rates for GCs with $\alpha_{3}=1.6$ and $3$, relative to the
rates in canonical GCs with $\alpha_{3}=2.3$, are then roughly given by

$\frac{\Gamma_{1.6}}{\Gamma_{2.3}}\gtrsim 1.5\,\frac{f_{{\rm GC},1.6}}{f_{{\rm
GC},2.3}}\,,\ \ \ \ \ \frac{\Gamma_{3.0}}{\Gamma_{2.3}}\approx
0.25\,\frac{f_{{\rm GC},3.0}}{f_{{\rm GC},2.3}}\,.$ (4)

Note again that the above ratio for clusters with $\alpha_{3}=1.6$ is a lower
limit. Both ratios are also very approximate since we neither consider
interdependence of $M_{i}$ and $\alpha_{3}$ nor any cumulative dependency of
these ratios on the choices of initial parameters, such as virial radius,
Galactocentric distance, metallicity, and binary fraction. Furthermore,
$f_{x}$ is notoriously uncertain and could depend on redshift, metallicity,
and time of birth. Nevertheless, these estimates qualitatively suggest that if
$f_{1.6}>f_{2.3}\gtrsim f_{3.0}$, then the contribution to the cosmological
BBH merger rate from clusters could be a bit higher than currently estimated
(e.g., by Kremer et al., 2020b). In turn, if $f_{3.0}>f_{2.3}\gtrsim f_{1.6}$,
then those rate estimates could be too high already. Overall, uncertainties in
the IMF may contribute about an order of magnitude to the uncertainty on the
dynamical BBH merger rate (see also Chatterjee et al., 2017).

Finally, it is worth commenting briefly on black hole-neutron star (BH-NS) and
binary neutron star (BNS) mergers. In clusters that retain a BH population,
NSs are pushed out of the dense, BH-filled core due to mass segregation. Hence
the merger likelihood is significantly smaller for BH-NS and BNS mergers than
for BBH mergers in most GCs (Ye et al., 2020). It is therefore unsurprising
that all our models produce insignificant numbers of BH-NS and BNS mergers.
While our models with $\alpha_{3}=3$ – and therefore very few formed BHs – do
exhibit about four times more frequent binary NS mergers than the models with
$\alpha_{3}=2.3$, the numbers are still far too low to account for the
LIGO/Virgo estimated rates (Ye et al., 2020; Fragione & Banerjee, 2020).

Figure 4.— Cumulative density functions (CDFs) of the component masses for
binary black hole mergers in star clusters with different values of
$\alpha_{3}$ (blue $1.6$, black $2.3$, red $3.0$) and different initial number
of particles ($N_{i}$). Different line styles represent different $N_{i}$, as
in Figure 1.

###  3.4. BH and BBH merger masses

The average stellar mass is higher in clusters with top-heavy IMFs and lower
in clusters with top-light IMFs. So, it is unsurprising that the average
masses of BHs, BBH merger components, and BBH merger products in our models
all increase with decreasing $\alpha_{3}$. Averaging across all our models for
each $\alpha_{3}$ value, we find a mean BH mass at formation of
$(16,12,9)\,\rm{M_{\odot}}$ for $\alpha_{3}=(1.6,2.3,3)$, respectively. The
corresponding mean masses for BBH merger components and merger products are
$(38,26,21)\,\rm{M_{\odot}}$ and $(74,51,40)\,\rm{M_{\odot}}$. This trend is
further exhibited in Figure 4 showing the mass distribution of BBH merger
components. The vertical jump in the distributions at $M=40.5\,{\rm
M_{\odot}}$ is due to the pile-up at the start of the ‘mass gap’ of BHs formed
from massive stars via pulsational pair-instability supernovae (PPISNe). CMC
assumes this gap begins at $M>40.5\,{\rm M_{\odot}}$ (for details, see
Belczynski et al., 2016), extending to around $120\,{\rm M_{\odot}}$. Hence,
the first vertical jump indicates that clusters with top-heavy IMFs form
significantly more BHs via PPISNe than clusters with canonical IMFs. A second,
smaller jump in the distributions for top-heavy IMFs at $M\approx 77\,{\rm
M_{\odot}}$ arises from 2nd-generation mergers with a component produced in an
earlier merger between two such PPISNe-generated BHs.

The anti-correlation between the average mass of BBH merger components and
$\alpha_{3}$ arises for three reasons. First, clusters with top-heavy IMFs
produce higher-mass stars, increasing not only the number of BHs formed, but
also their average mass. Second, for initial $N$ held constant, clusters with
top-heavy IMFs also have greater total mass and higher collision rates in
their densely populated cores. This collision rate enhancement is obvious in
Table 1, column 7, listing the number of mass gap BBH merger components formed
from the stellar product of collisions (e.g., Di Carlo et al., 2020; Kremer et
al., 2020c). While the column only lists collisionally formed merger
components in the mass gap, the overall collision rate scales similarly with
decreasing $\alpha_{3}$. Finally, the increased merger frequency in clusters
with top-heavy IMFs itself increases the average mass of the components by
increasing the chances that they will have already experienced a merger (e.g.,
Fragione et al., 2020a, b; Rodriguez et al., 2020).

Indeed, the cumulative number of 2nd-generation or higher (2G+) BBH mergers is
also enhanced in clusters with top-heavy IMFs, as seen in Table 1, column 8,
listing the number of 2G+ mergers with a component in the mass gap. For
clusters with top-heavy IMFs, this is identical to the total number of 2G+ BBH
mergers, though 2G+ mergers with neither component in the mass gap are common
for higher $\alpha_{3}$. Overall, 2G+ merger totals roughly double with each
decrease in $\alpha_{3}$ from $3$ to $2.3$ to $1.6$, while the numbers of 2G+
mergers with a mass gap component roughly quadruple.

Though we find no trend between initial $N$ and merger component or product
mass for $\alpha_{3}=2.3$ or $3$, there is a clear correlation for
$\alpha_{3}=1.6$. In this case, for $N_{i}=(4,8,16)\times 10^{5}$, the average
BBH merger product mass is $(54,62,79)\,\rm{M_{\odot}}$, respectively. This
rise is in part due to 2G+ mergers, which account for $(7\%,12\%,18\%)$ of the
BBH mergers in these respective models. In general, we find that clusters with
top-heavy IMFs and high initial $N$ are especially good at efficiently
producing many 2G+ mergers and many BHs in the mass gap.

The increased collision and 2G+ merger rates in clusters with top-heavy IMFs
can also result in the formation of ‘intermediate-mass’ BHs (IMBHs), which we
define as BHs with masses exceeding $100\,{\rm M_{\odot}}$. In particular, our
highest-$N$ model with $\alpha_{3}=1.6$ (model 3 in Table 1) produce 46 IMBHs,
including one with mass $M=537\,{\rm M_{\odot}}$. This particular IMBH formed
in the merger of a $122\,{\rm M_{\odot}}$ IMBH with a $426\,{\rm M_{\odot}}$
IMBH, which itself formed in the merger of two IMBHs with $M\approx 200\,{\rm
M_{\odot}}$. While we intend to explore IMBH-formation in GCs with top-heavy
IMFs more thoroughly in future work, we for now direct the reader’s attention
to our collaboration’s recent study on IMBH formation in GCs (González et al.,
2020).

##  4\. Discussion and Conclusions

Observations provide evidence that the stellar IMF may not be universal. A
non-canonical IMF can greatly affect a star cluster’s dynamical evolution,
especially since its high end determines how many BHs form within, regulating
the cluster’s energy budget and dynamical clock (Breen & Heggie, 2013;
Chatterjee et al., 2017; Kremer et al., 2020a; Wang, 2020). In this Letter, we
have extended the CMC Cluster Catalog (Kremer et al., 2020b) to examine how a
varying IMF affects the evolution of star clusters and their BH populations.

We have shown that massive star clusters with top-heavy IMFs (low
$\alpha_{3}$) are likely to lose most of their mass within a few Gyr, assuming
they have low-to-average mass and Galactocentric distance for typical Milky
Way GCs. The rapid mass loss during dissolution occurs in stages, first driven
by stellar winds and dynamical ejections, then by evaporation of the halo
through the tidal boundary. Extensive BH-burning enhances the latter stage in
clusters with top-heavy IMFs, which produce many BHs. Such clusters evolve
through a point where they consist mostly of BHs by mass (and up to at least
$3\%$ by number).333Just like the fictional Maw cluster (Anderson, 1994), this
‘dark cluster’ stage of a GC born with a top-heavy IMF may be short-lived due
to the rapid pace of tidal evaporation. Further study with direct $N$-body
methods (e.g., Banerjee & Kroupa, 2011) is required to fully understand the
evolution of these clusters. Initializing direct $N$-body simulations with the
pre-dissolution states of CMC models could be especially useful in such an
evaluation across the full cluster mass distribution.

We note that the processes described above are also affected by the choices of
initial cluster metallicity and natal kick distribution (Chatterjee et al.,
2017). BHs in metal-rich clusters have lower mass and do not inject as much
energy into the BH-burning process as BHs in metal-poor clusters. Thus, metal-
rich clusters typically have higher densities and dispersion velocities,
therefore processing the BH population on shorter timescales and disrupting
more binaries. Meanwhile, high natal kicks will eject most BHs from the
cluster during formation. In such a case, BHs and BBHs are expected to be
small in number regardless of the IMF. Top-heavy clusters with higher natal
kicks also live longer; while they experience less tidal mass loss driven by
BH-burning, these clusters lose even more mass due to kick-driven dynamical
ejections of BHs (Chatterjee et al., 2017). For a recent study of top-heavy
cluster evolution featuring different natal kick assumptions, see Haghi et al.
(2020, published during review of this paper), who examine lower-mass clusters
with direct $N-$body simulations incorporating gas expulsion physics. Notably,
they find that top-heavy clusters, albeit drastically reduced in mass, may
well survive to the present day if born with masses above $M_{i}\gtrsim
7\times 10^{5}\,{\rm M_{\odot}}$. For times not too close to disruption, where
CMC’s assumptions start to be challenged (see Section 2), our results are
encouragingly compatible with those of Haghi et al. (2020).

Regardless of their long-term evolution and stability, we have also shown that
clusters with top-heavy IMFs – and correspondingly high BH production – may
contribute significantly to the present-day binary BH merger rate. Even though
these clusters rapidly lose most of their mass within a few Gyr, mergers from
ejected BBHs continue to contribute at later times (Fragione & Kocsis, 2018),
at rates comparable to or greater than those for clusters with canonical
Kroupa IMFs. The rate of 2nd-generation mergers with component masses in the
mass gap may be especially enhanced in top-heavy GCs, motivating the existence
of more GW190521-like mergers (Abbott et al., 2020a). In addition, the
enhancement of collision rates and multiple-generation mergers in top-heavy
GCs may also lead to the formation of IMBHs and even IMBH-IMBH mergers, as
demonstrated in one of our models.

In general, we have shown that the high-mass slope of the cluster birth IMF
may significantly impact the exact contribution to the cosmological BBH merger
rate due to cluster dynamics. Specifically, if a large fraction of clusters
were born with top-heavy IMFs, the cluster-dynamics merger rate may be
somewhat enhanced relative to recent estimates (e.g., Kremer et al., 2020b).
In turn, if a large fraction of clusters were born with top-light IMFs, the
cluster-dynamics merger rate may be significantly reduced. With future
observations of gravitational waves providing unique information on the BBH
merger rate (Abbott et al., 2020b), it may be possible to leverage this
understanding to better constrain the IMFs of old star clusters.

## Acknowledgements

This work was supported by NSF grant AST-1716762 and through the computational
resources and staff contributions provided for the Quest high-performance
computing facility at Northwestern University. NW acknowledges support from
the CIERA Riedel Family Graduate Fellowship as well as the NSF GK-12
Fellowship Program under Grant DGE-0948017. GF acknowledges support from a
CIERA Fellowship. KK is supported by an NSF Astronomy and Astrophysics
Postdoctoral Fellowship under award AST-2001751. SC acknowledges support from
the Department of Atomic Energy, Government of India, under Project No.
12-R&D-TFR-5.02-0200.

## References

  * Abbott et al. (2020a) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020a, Phys. Rev. Lett., 125, 101102
  * Abbott et al. (2020b) Abbott, R., et al. 2020b, arXiv e-prints, arXiv:2010.14527
  * Adams & Fatuzzo (1996) Adams, F. C., & Fatuzzo, M. 1996, ApJ, 464, 256
  * Adams & Laughlin (1996) Adams, F. C., & Laughlin, G. 1996, ApJ, 468, 586
  * Anderson (1994) Anderson, K. J. 1994, Star Wars: The Jedi Academy Trilogy, Vol. 1, Jedi Search (New York: Bantam Spectra)
  * Askar et al. (2017) Askar, A., Szkudlarek, M., Gondek-Rosińska, D., Giersz, M., & Bulik, T. 2017, MNRAS, 464, L36
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Banerjee (2017) Banerjee, S. 2017, MNRAS, 467, 524
  * Banerjee & Kroupa (2011) Banerjee, S., & Kroupa, P. 2011, ApJL, 741, L12
  * Bartko et al. (2010) Bartko, H., Martins, F., Trippe, S., et al. 2010, ApJ, 708, 834
  * Belczynski et al. (2016) Belczynski, K., Heger, A., Gladysz, W., et al. 2016, A&A, 594, A97
  * Bennett et al. (2014) Bennett, C. L., Larson, D., Weiland, J. L., & Hinshaw, G. 2014, ApJ, 794, 135
  * Berti et al. (2007) Berti, E., Cardoso, V., Gonzalez, J. A., et al. 2007, Phys. Rev. D, 76, 064034
  * Breen & Heggie (2013) Breen, P. G., & Heggie, D. C. 2013, MNRAS, 432, 2779
  * Chatterjee et al. (2010) Chatterjee, S., Fregeau, J. M., Umbreit, S., & Rasio, F. A. 2010, ApJ, 719, 915
  * Chatterjee et al. (2017) Chatterjee, S., Rodriguez, C. L., & Rasio, F. A. 2017, ApJ, 834, 68
  * Chatterjee et al. (2013) Chatterjee, S., Umbreit, S., Fregeau, J. M., & Rasio, F. A. 2013, MNRAS, 429, 2881
  * Chernoff & Weinberg (1990) Chernoff, D. F., & Weinberg, M. D. 1990, ApJ, 351, 121
  * Clark (1975) Clark, G. 1975, ApJ, 199, L143
  * Dabringhausen et al. (2009) Dabringhausen, J., Kroupa, P., & Baumgardt, H. 2009, MNRAS, 394, 1529
  * De Marchi et al. (2017) De Marchi, G., Panagia, N., & Beccari, G. 2017, ApJ, 846, 110
  * De Marchi et al. (2007) De Marchi, G., Paresce, F., & Pulone, L. 2007, ApJL, 656, L65
  * Dehnen & Binney (1998) Dehnen, W., & Binney, J. 1998, MNRAS, 294, 429
  * Di Carlo et al. (2020) Di Carlo, U. N., Mapelli, M., Bouffanais, Y., et al. 2020, MNRAS, 497, 1043
  * Downing et al. (2010) Downing, J. M. B., Benacquista, M. J., Giersz, M., & Spurzem, R. 2010, MNRAS, 407, 1946
  * Downing et al. (2011) —. 2011, MNRAS, 416, 133
  * Duquennoy & Mayor (1991) Duquennoy, A., & Mayor, M. 1991, A&A, 500, 337
  * El-Badry et al. (2019) El-Badry, K., Quataert, E., Weisz, D. R., Choksi, N., & Boylan-Kolchin, M. 2019, MNRAS, 482, 4528
  * Fragione & Banerjee (2020) Fragione, G., & Banerjee, S. 2020, ApJL, 901, L16
  * Fragione & Kocsis (2018) Fragione, G., & Kocsis, B. 2018, Phys. Rev. Lett., 121, 161103
  * Fragione et al. (2020a) Fragione, G., Loeb, A., & Rasio, F. A. 2020a, ApJL, 895, L15
  * Fragione et al. (2020b) —. 2020b, ApJL, 902, L26
  * Fregeau et al. (2003) Fregeau, J. M., Gürkan, M. A., Joshi, K. J., & Rasio, F. A. 2003, ApJ, 593, 772
  * Giersz et al. (2019) Giersz, M., Askar, A., Wang, L., et al. 2019, MNRAS, 487, 2412
  * Giesler et al. (2018) Giesler, M., Clausen, D., & Ott, C. D. 2018, MNRAS, 477, 1853
  * González et al. (2020) González, E., Kremer, K., Chatterjee, S., et al. 2020, arXiv e-prints, arXiv:2012.10497
  * Haghi et al. (2017) Haghi, H., Khalaj, P., Hasani Zonoozi, A., & Kroupa, P. 2017, ApJ, 839, 60
  * Haghi et al. (2020) Haghi, H., Safaei, G., Zonoozi, A. H., & Kroupa, P. 2020, ApJ, 904, 43
  * Harris et al. (2014) Harris, W. E., Morningstar, W., Gnedin, O. Y., et al. 2014, ApJ, 797, 128
  * Heggie & Hut (2003) Heggie, D., & Hut, P. 2003, The Gravitational Million-Body Problem: A Multidisciplinary Approach to Star Cluster Dynamics (Cambridge: Cambridge Univ. Press)
  * Heggie (1975) Heggie, D. C. 1975, MNRAS, 173, 729
  * Hénon (1971a) Hénon, M. 1971a, Astrophysics and Space Science, 13, 284
  * Hénon (1971b) —. 1971b, Astrophysics and Space Science, 14, 151
  * Ivanova et al. (2008) Ivanova, N., Heinke, C. O., Rasio, F. A., Belczynski, K., & Fregeau, J. M. 2008, MNRAS, 386, 553
  * Joshi et al. (2001) Joshi, K. J., Nave, C. P., & Rasio, F. A. 2001, ApJ, 550, 691
  * Joshi et al. (2000) Joshi, K. J., Rasio, F. A., Zwart, S. P., & Portegies Zwart, S. 2000, ApJ, 540, 969
  * King (1962) King, I. 1962, AJ, 67, 471
  * Kremer et al. (2018) Kremer, K., Chatterjee, S., Rodriguez, C. L., & Rasio, F. A. 2018, ApJ, 852, 29
  * Kremer et al. (2020a) Kremer, K., Ye, C. S., Chatterjee, S., Rodriguez, C. L., & Rasio, F. A. 2020a, in Star Clusters: From the Milky Way to the Early Universe, ed. A. Bragaglia, M. Davies, A. Sills, & E. Vesperini, Vol. 351, 357
  * Kremer et al. (2019) Kremer, K., Rodriguez, C. L., Amaro-Seoane, P., et al. 2019, Phys. Rev. D, 99, 063003
  * Kremer et al. (2020b) Kremer, K., Ye, C. S., Rui, N. Z., et al. 2020b, ApJS, 247, 48
  * Kremer et al. (2020c) Kremer, K., Spera, M., Becker, D., et al. 2020c, ApJ, 903, 45
  * Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231
  * Lada & Lada (2003) Lada, C. J., & Lada, E. A. 2003, ARA& A, 41, 57
  * Larson (1998) Larson, R. B. 1998, MNRAS, 301, 569
  * Lousto et al. (2010) Lousto, C. O., Campanelli, M., Zlochower, Y., & Nakano, H. 2010, Classical and Quantum Gravity, 27, 114006
  * Lyne et al. (1987) Lyne, A., Brinklow, A., Middleditch, J., et al. 1987, Nature, 328, 399
  * Mackey et al. (2007) Mackey, A. D., Wilkinson, M. I., Davies, M. B., & Gilmore, G. F. 2007, MNRAS, 379, L40
  * Mackey et al. (2008) —. 2008, MNRAS, 386, 65
  * Marks et al. (2012) Marks, M., Kroupa, P., Dabringhausen, J., & Pawlowski, M. S. 2012, MNRAS, 422, 2246
  * Morscher et al. (2015) Morscher, M., Pattabiraman, B., Rodriguez, C., Rasio, F. A., & Umbreit, S. 2015, ApJ, 800, 9
  * Pattabiraman et al. (2013) Pattabiraman, B., Umbreit, S., Liao, W.-k., et al. 2013, ApJS, 204, 15
  * Rodriguez et al. (2015) Rodriguez, C. L., Morscher, M., Pattabiraman, B., et al. 2015, Phys. Rev. Lett., 115, 051101
  * Rodriguez et al. (2019) Rodriguez, C. L., Zevin, M., Amaro-Seoane, P., et al. 2019, Phys. Rev. D, 100, 043027
  * Rodriguez et al. (2020) Rodriguez, C. L., Kremer, K., Grudić, M. Y., et al. 2020, ApJL, 896, L10
  * Sigurdsson & Phinney (1995) Sigurdsson, S., & Phinney, E. S. 1995, ApJS, 99, 609
  * Tichy & Marronetti (2008) Tichy, W., & Marronetti, P. 2008, Phys. Rev. D, 78, 081501
  * Verbunt et al. (1984) Verbunt, F., van Paradijs, J., & Elson, R. 1984, MNRAS, 210, 899
  * Wang (2020) Wang, L. 2020, MNRAS, 491, 2413
  * Ye et al. (2020) Ye, C. S., Fong, W.-f., Kremer, K., et al. 2020, ApJL, 888, L10
  * Ye et al. (2019) Ye, C. S., Kremer, K., Chatterjee, S., Rodriguez, C. L., & Rasio, F. A. 2019, ApJ, 877, 122


# Indications for a Nonzero Lepton Asymmetry in the Early Universe

Anne-Katherine Burns<EMAIL_ADDRESS>Tim M.P. Tait<EMAIL_ADDRESS>Department
of Physics and Astronomy, University of California, Irvine, CA 92697 USA
Mauro Valli<EMAIL_ADDRESS>C.N. Yang Institute for Theoretical
Physics, Stony Brook University, Stony Brook, NY 11794, USA

###### Abstract

The recent measurement of helium-4 from the near-infrared spectroscopy of
extremely metal-poor galaxies by the Subaru Survey may point to a new puzzle
in the Early Universe. We exploit this new helium measurement together with
the percent-level determination of primordial deuterium, to assess indications
for a non-vanishing lepton asymmetry during the Big Bang Nucleosynthesis (BBN)
era, paying particular attention to the role of uncertainties in the nuclear
reaction network. A cutting-edge Bayesian analysis of BBN data jointly with
information from the Cosmic Microwave Background suggests the existence of a
nonzero lepton asymmetry at around the $2\sigma$ level, providing a hint for
cosmology beyond $\Lambda$CDM. We discuss conditions for a large total lepton
asymmetry to be consistently realized in the Early Universe.

††preprint: UCI-TR-2022-10††preprint: YITP-SB-2022-22

##  I Introduction

Cosmological observations from the Early Universe provide an invaluable probe
of Physics Beyond the Standard Model (BSM). Observations of the Cosmic
Microwave Background (CMB), epitomized by the Planck mission Aghanim  _et al._
(2020) and further developed e.g. by the ACT Aiola  _et al._ (2020) and SPT
Dutcher  _et al._ (2021) collaborations, paint a picture of a Universe that is
dominated by non-baryonic dark energy and dark matter, well-described by the
$\Lambda$CDM model Baumann  _et al._ (2016); Komatsu (2022); Chang  _et al._
(2022). Equipped with the CMB inference of the small cosmological baryonic
abundance, $\Omega_{B}\sim 4\%$, the theory of Big Bang Nucleosynthesis (BBN)
within the Standard Model (SM) of Particle Physics is highly predictive, and
confronted with accurate measurements of primeval elements such as the mass
density fraction of helium-4, ${\rm Y_{P}}$, and the relative abundance of
deuterium, ${\rm D/H}$, offers important constraints on New Physics (NP)
Sarkar (1996); Steigman (2007); Iocco  _et al._ (2009); Pospelov and Pradler
(2010) active during the first few minutes of the lifetime of the Universe
Serpico and Raffelt (2004); Boehm  _et al._ (2013); Sabti  _et al._ (2020);
Giovanetti  _et al._ (2021).

At present, measurements of deuterium in quasar absorption spectra provide the
best proxy for the determination of a primordial abundance. The most recent
measurements from damped Lyman-$\alpha$ systems achieve better than $1\%$
precision Cooke  _et al._ (2016); Riemer-Sørensen  _et al._ (2017); Cooke  _et
al._ (2018), yielding a weighted average of ${\rm D/H}\times 10^{5}=2.547\pm
0.025$ Zyla  _et al._ (2020). This remarkable precision appears to be in
tension with the SM at about the 2$\sigma$ level Pitrou  _et al._ (2021a),
although this remains under debate Pisanti  _et al._ (2021); Yeh  _et al._
(2021) in light of the uncertainties plaguing our understanding of the key
nuclear reactions involved. This highlights the primary importance to assess
the impact of uncertainties in the nuclear network rates on the predictions
from BBN Pitrou  _et al._ (2021b). A notable recent advance in this direction
is the improved determination of the $\textrm{D}(p,\gamma)^{3}\textrm{He}$
rate by the LUNA collaboration Mossa  _et al._ (2020), which has an important
impact on BBN constraints from primordial deuterium on various NP scenarios
Sabti  _et al._ (2021).

Figure 1: 68% and 90% joint probability distribution of the primordial
chemical potential of neutrinos, $\mu_{\nu}$, normalized to the neutrino
temperature $T_{\nu}$ and the number of extra relativistic degrees of freedom
in the Early Universe, $\Delta N_{\rm eff}$, from a state-of-the-art analysis
of BBN and CMB data. The red and blue contours indicate the results for two
different sets of nuclear uncertainties adopted in the analysis as explained
in the text. Magenta dashed lines highlight the $\Lambda$CDM prediction.

The recent near-infrared observation of 10 extremely metal-poor galaxies
(EMPGs) by the Subaru Survey Matsumoto  _et al._ (2022) points to even more
puzzling mysteries. Spectroscopic observations of EMPGs provide a crucial
input to the inference of ${\rm Y_{P}}$, because they host the gas of nebulae
resembling extraordinarily pristine environments which allow for a more
accurate extrapolation of the helium density to zero metallicity. Combined
with the pre-existing data from 3 EMPGs and 51 metal-poor galaxies Hsyu  _et
al._ (2020) and measurements of the He $\lambda$10830 infrared emission line
(relevant for parameter-degeneracy breaking Aver  _et al._ (2015)), the 10
Subaru EMPGs yield a determination of primordial helium-4 of ${\rm
Y_{P}}=0.2379^{+0.0031}_{-0.0030}$, in sharp contrast with the PDG value ${\rm
Y_{P}}=0.245\pm 0.003$ Zyla  _et al._ (2020), and well below the SM prediction
Pitrou  _et al._ (2021a); Pisanti  _et al._ (2021); Yeh  _et al._ (2021),
naively a 3$\sigma$-level discrepancy.

Ref. Matsumoto  _et al._ (2022) took the first steps toward an interpretation
of this helium anomaly in terms of a BSM fit where the standard theory has
been extended by extra-relativistic degrees of freedom, $\Delta N_{\rm eff}$,
as well as a nonzero electron neutrino asymmetry, $\xi_{\nu_{e}}$, while
simply anchoring $\Omega_{B}h^{2}$ to the most precise determination derived
by Planck Aghanim  _et al._ (2020). In view of the situation depicted for both
helium and deuterium, in this Letter we revise the inference of a lepton
asymmetry $\xi_{\nu}$ in the Early Universe, as well as on $\Delta N_{\rm
eff}$, paying attention to the details of a joint likelihood analysis of BBN
and CMB data as recently carefully formulated in Sabti  _et al._ (2020);
Giovanetti  _et al._ (2021). To this aim, we perform a Bayesian analysis that
takes into account the theory uncertainties pivotal for unbiased conclusions
based on the use of the new public code for state-of-the-art investigations of
(BSM) physics in the Early Universe – PRyMordial – presented in detail in a
companion paper A.-K. Burns, T. M.P. Tait, and M. Valli .

Figure 1 shows our key result: a nonzero asymmetry in the lepton sector is
favored by current data at $90\%$ probability, irrespective of nuclear
modeling. At the same time, with somewhat more dependence on the treatment of
the BBN nuclear network rates, the data may provide some hint for new sources
of relativistic degrees of freedom. In the following, we detail the outcome of
the present analysis and comment on the possible NP implications.

##  II Primordial Lepton Asymmetries

Electric charge neutrality of the Early Universe does not allow for a large
primordial asymmetry in the charged lepton sector, which is constrained to be
(at most) of the order of the baryon-to-photon ratio $\eta_{B}\equiv
n_{B}/n_{\gamma}\sim\mathcal{O}(10^{-10})$ Kolb and Turner (1990); Rubakov and
Gorbunov (2017). Nevertheless, a large cosmic asymmetry can be hidden in the
neutrino sector Simha and Steigman (2008):

$\eta_{L}\equiv\frac{1}{n_{\gamma}}\sum_{i=e,\mu,\tau}(n_{\nu_{i}}-n_{\bar{\nu}_{i}})\simeq\frac{\pi^{2}}{33\zeta(3)}(\xi_{\nu_{e}}+\xi_{\nu_{\mu}}+\xi_{\nu_{\tau}})\
,$ (1)

where $n_{\gamma}$ is the photon number density, $n_{\nu_{i}}$ the flavor $i$
neutrino density, and $\xi_{\nu_{i}}\equiv\mu_{\nu_{i}}/T_{\nu_{i}}$ are the
degeneracy parameters defined as the chemical potential for each neutrino
normalized to its temperature, which encode the relevant lepton asymmetries
today. Eq. (1) assumes $T_{\nu_{i}}/T_{\gamma}=(4/11)^{1/3}$, which is a good
approximation given the modest impact of non-instantaneous neutrino decoupling
and tiny departures from the Fermi-Dirac distributions in relativistic freeze-
out Akita and Yamaguchi (2020); Froustey  _et al._ (2020); Bennett  _et al._
(2021). It is further relevant that non-zero neutrino chemical potentials play
a marginal role in SM neutrino decoupling Grohs  _et al._ (2016); Escudero
Abenza (2020).

Eq. (1) further implements the condition $|\xi_{\nu_{i}}|<1$, as
$\mathcal{O}(1)$ degeneracy parameters were probed by early-stage CMB
observations about two decades ago Kinney and Riotto (1999); Lesgourgues and
Pastor (1999), and now are robustly Oldengott and Schwarz (2017) ruled out
(irrespective of the lepton flavor Abazajian  _et al._ (2002); Castorina  _et
al._ (2012)). In fact, a nonzero chemical potential for the $i$-flavored
neutrino would yield a contribution to the total radiation density (relative
to photons) of

$\frac{\Delta\rho_{\rm
rad}}{\rho_{\gamma}}\simeq\frac{15}{4\pi^{2}}\left(\frac{4}{11}\right)^{4/3}\xi_{\nu_{i}}^{2}\
,$ (2)

and would increase the expansion rate of the Universe, resulting in a positive
shift of $N_{\rm eff}$ that would delay the time of matter-radiation equality
which is tightly constrained by the CMB acoustic peaks.

From the Planck constraint on $N_{\rm eff}$ adopting the likelihood analysis
including TTTEEE and low-$\ell$ measurements, as well as baryonic acoustic
oscillations (BAO) and lensing data, and assuming a flat prior on $\rm Y_{\rm
P}$, one may derive a simple upper-bound on the degeneracy parameters purely
driven by the CMB. In particular, for $N_{\rm eff}=2.97\pm 0.29$ (68%
probability interval) Aghanim  _et al._ (2020, 2019a), considering the SM
prediction $N_{\rm eff}=3.044$ (known better than the per-mille level) Akita
and Yamaguchi (2020); Froustey  _et al._ (2020); Bennett  _et al._ (2021), the
1$\sigma$ upper-bound is:

$\xi_{\nu_{e}}^{2}+\xi_{\nu_{\mu}}^{2}+\xi_{\nu_{\tau}}^{2}\lesssim 0.5\ ,$
(3)

implying the conservative constraint $|\xi_{\nu_{i}}|\lesssim 0.71$, valid for
each flavor individually (see also Ref. Barenboim  _et al._ (2017)). Since the
onset of neutrino oscillations is expected to occur around $T_{\nu}\sim 10$
MeV, flavor equilibration in the muon-tau sector is predicted to be complete
by the time of neutrino decoupling ($T_{\nu}\sim 2$ MeV) Dolgov  _et al._
(2002); Pastor  _et al._ (2009), and the conservative CMB bound of Eq. (3)
becomes slightly tighter for the 2nd and 3rd generation $\nu$ asymmetries:

$|\xi_{\nu_{\mu,\tau}}|\lesssim 0.5\ .$ (4)

BBN can place stronger constraints on the electron-neutrino asymmetry (by
about an order of magnitude Simha and Steigman (2008); Fields  _et al._
(2020)) largely because an electron-neutrino asymmetry at the time of BBN
affects the $\beta$ equilibrium of weak interactions controlling the neutron-
proton conversion. A positive (negative) value of $\xi_{\nu_{e}}$ acts through
the equilibrium reactions $n\,\nu_{e}\leftrightarrow p\,e^{-}$,
$p\,\bar{\nu}_{e}\leftrightarrow n\,e^{+}$ and neutron decay to reduce
(enhance) the neutron-to-proton ratio:

$(n_{n}/n_{p})\big{|}_{\rm
eq.}\simeq\exp\left(-\mathcal{Q}/T_{\gamma}-\xi_{\nu_{e}}\right)\ ,$ (5)

where $\mathcal{Q}\equiv m_{n}-m_{p}=1.293$ MeV is the neutron-proton mass
difference. While light primordial abundances like deuterium are particularly
sensitive to $\Omega_{B}$, from Eq. (5) it follows that helium-4 – which
depends crucially on the amount of neutrons at the time where deuterium is no
longer photo-dissociated Kolb and Turner (1990); Rubakov and Gorbunov (2017) –
can be regarded as a sensitive primordial leptometer.

In Ref. Pitrou  _et al._ (2018) a combined analysis of the helium and
deuterium PDG values together with a Gaussian prior on $\Omega_{B}h^{2}$ from
the CMB yields the precise determination: $\xi_{\nu}=0.001\pm 0.016$,
consistent with zero. By assuming full equilibration of lepton flavor
asymmetries due to neutrino oscillations, this inference is more stringent
than the bound outlined in Eq. (4). Nevertheless, a recent state-of-the-art
investigation in Ref. Froustey and Pitrou (2022) indicates that the degree to
which full flavor equilibration is realized during the BBN era sensitively
depends on the PMNS mixing angle $\theta_{13}$ and on the initially generated
values of the degeneracy parameters. In the following, we revisit the
determination of $\xi_{\nu}$ in light of the newly measured helium-4 mass
fraction from EMPGs as reported in Ref. Matsumoto  _et al._ (2022). While
$\xi_{\nu}=\xi_{\nu_{e,\mu,\tau}}$ may be achieved in the Early Universe, one
should bear in mind that the conservative interpretation of our main finding
in Fig. 1 applies only to $\xi_{\nu}=\xi_{\nu_{e}}$, i.e. the primordial
electron-neutrino asymmetry probed by BBN via $\beta$ equilibrium, Eq. (5).

##  III Methodology

Figure 2: Probability density function (p.d.f.) for the primordial light
elements analyzed in this study. In the left panel, the p.d.f. for helium-4,
$Y_{\rm P}$, as precisely predicted in the SM according to two different set
of nuclear uncertainties and adopting the determination of the cosmological
baryon-to-photon ratio from the fit of CMB data within $\Lambda$CDM (color
code similar to Figure 1). In the same panel, the outcome from the joint fit
to BBN and CMB likelihoods in the BSM scenario where $\xi_{\nu}$ and $\Delta
N_{\rm eff}$ are consistently allowed to differ from their $\Lambda$CDM limit.
In the right panel, the same set of p.d.f.s is shown for the deuterium. In
both panels, vertical dark green bands correspond to the 1$\sigma$ interval
for the BBN measurements employed in the analysis. In the left one, the PDG
2021 recommended value for helium-4 is also reported, in agreement with the SM
prediction.

Our computation of BBN abundances via the PRyMordial A.-K. Burns, T. M.P.
Tait, and M. Valli  code, proceeds through 3 main steps:

  1. 1)

Solving for the thermal background comprising the electron-photon plasma and
neutrino species, to pin down the scale factor of the universe as a function
of $T_{\gamma}$ and time, as well as the relation $T_{\nu}(T_{\gamma})$;

  2. 2)

Computing the amplitudes for neutron-proton conversion via the weak effective
theory, determining the initial conditions for the BBN nuclear network;

  3. 3)

Evolving the primordial abundances down to the end of BBN according to a set
of reactions that are key for precise helium-4 and deuterium predictions.

For 1) we base our computation on the approach proposed in Escudero (2019) and
further developed in Escudero Abenza (2020) (see also Chu  _et al._ (2022)).
It consists in solving the Boltzmann equations for the electron-photon plasma
and neutrinos assuming a thermal distribution for the species, including NLO
QED corrections for the plasma Bennett  _et al._ (2020) as well as non-
instantaneous decoupling effects for the neutrino sector Escudero Abenza
(2020). For our purposes, it suffices to describe the neutrino sector by a
common temperature $T_{\nu}$, yielding the SM prediction $N_{\rm eff}=3.045$,
differing from the most refined prediction in Akita and Yamaguchi (2020);
Froustey  _et al._ (2020); Bennett  _et al._ (2021) only at the per-mille
level, well within current and future observational sensitivity Ade  _et al._
(2019); Abazajian  _et al._ (2022); Aiola  _et al._ (2022). A non-zero
chemical potential for neutrinos would influence our analysis of the thermal
background via Eq. (2). If full neutrino equilibration is achieved, we find
a-posteriori a contribution to the radiation density that would be totally
negligible. Nevertheless, in our BSM analysis we also account for the
possibility of a nonzero lepton asymmetry
$|\xi_{\nu_{\mu,\tau}}|\gg|\xi_{\nu_{e}}|$ by varying $\Delta N_{\rm eff}$.
Note that from our bound in Eq. (4), a muon-tau neutrino asymmetry can induce
a maximal shift $\Delta N_{\rm eff}\sim 0.1$.

Moving to 2), we compute $n\leftrightarrow p$ matrix elements beyond the Born
approximation Sirlin (1967), namely including isospin-breaking contributions
like finite-mass Brown and Sawyer (2001) and QED Czarnecki  _et al._ (2004);
Ivanov  _et al._ (2017) corrections, as well as finite-temperature effects
Lopez  _et al._ (1997), following the implementation carried out in Pitrou
_et al._ (2018). Most importantly, we evaluate weak-interaction rates
integrating over nucleon thermal distributions with chemical potential
$\mu_{\mathcal{Q}}\equiv\mu_{n}-\mu_{p}=-\mu_{\nu_{e}}\neq 0$.

Finally, regarding 3), we proceed evolving the abundances according to the
network of thermonuclear reactions comprising the main processes listed in
Table 1 of Ref. Iliadis and Coc (2020) (plus
${}^{3}\textrm{He}(p,\gamma)^{4}\textrm{He}$, taken from Serpico  _et al._
(2004)), yielding state-of-the-art predictions for $Y_{\rm P}$ and ${\rm
D/H}$. In particular, for the radiative neutron capture rate we adopt the MCMC
result of Ref. Ando  _et al._ (2006), while in the treatment of the other 10
key reactions we distinguish two approaches:

  * $\circ$

PRIMAT driven: Nuclear rates are implemented according to the statistical
determination of Refs Descouvemont  _et al._ (2004); Iliadis  _et al._ (2016);
Iñesta Gómez  _et al._ (2017); de Souza  _et al._ (2019a, b); Moscoso  _et
al._ (2021), i.e. following theoretical ab-initio energy modeling tuned to
datasets for which an estimate of systematic errors is available Pitrou  _et
al._ (2018, 2021a, 2021b).

  * $\circ$

NACRE II driven: Nuclear rates are interpolated from the updated NACRE
compilation Xu  _et al._ (2013), comprising charged-particle-induced
reactions; for $\textrm{D}(p,\gamma)^{3}\textrm{He}$ we use the LUNA result
reported in Mossa  _et al._ (2020); for
${}^{7}\textrm{Be}(n,p)^{7}\textrm{Li}$ we adopt the baseline of Ref. Fields
_et al._ (2020).

We perform a Bayesian analysis of Early Universe data constructing the
cosmological test statistic:

$\textrm{TS}_{\rm cosmo}\equiv-2(\log\mathcal{L}_{\rm
CMB}+\log\mathcal{L}_{\rm BBN})\ ;$ (6)

the CMB likelihood explicitly reads:

$\log\mathcal{L}_{\rm
CMB}=-\frac{1}{2}\Delta\vec{v}^{\,T}\,\mathcal{C}^{-1}_{\rm
CMB}\,\Delta\vec{v}\ ,$ (7)

with $\Delta\vec{v}\equiv\vec{v}^{\,\rm th}-\vec{v}$, $\vec{v}=({\rm
Y_{P}},\Omega_{B}h^{2},N_{\rm eff})^{T}\ $, using mean and standard-deviation
values from the TTTEEE + low-$\ell$ \+ BAO + lensing Planck run varying also
${\rm Y_{P}}$, $N_{\rm eff}$ Aghanim  _et al._ (2020, 2019a) and also
retrieving correlations in $\mathcal{C}_{\rm CMB}$ from Aghanim  _et al._
(2019b). The BBN likelihood of our study corresponds to:

$\log\mathcal{L}_{\rm BBN}=-\frac{1}{2}\sum_{X}\left(\frac{X^{\rm
th}-X}{\sigma_{X}}\right)^{2}\ ,$ (8)

where $X=\\{{\rm Y_{P}},{\rm D/H}$}, and we use the measurements: ${\rm
Y_{P}}=0.2379(31)$ Matsumoto  _et al._ (2022), ${\rm D/H}=0.00002547(25)$ Zyla
_et al._ (2020).

The parameters we infer are varied according to uniform priors: $-2\leq\Delta
N_{\rm eff}\leq 2$, $\ -0.2\leq\xi_{\nu}\leq 0.2$, $1\leq(\eta_{B}\times
10^{10})\leq 10$ (using $\eta_{B}\times 10^{10}\simeq
273.748\,\,\Omega_{B}h^{2}$). We marginalize over the neutron lifetime and the
adopted nuclear uncertainties. From the PDG analysis Zyla  _et al._ (2020) we
assign the Gaussian prior: $\tau_{n}=(879.4\pm 0.6$) s to the neutron
lifetime. For the uncertainties in the nuclear rates, we assign log-normal
distributions following the method detailed in Ref. Coc  _et al._ (2014),
varying a total of 12 additional nuisance parameters.

We perform an MCMC analysis via the emcee Foreman-Mackey  _et al._ (2013)
package, using 60 walkers with 2100 steps each, discarding the first 700 steps
of each walker as burn-in. From the best-fit values minimizing
$\textrm{TS}_{\rm cosmo}$ we also compute for each scenario the Information
Criterion Akaike (1974); Gelman  _et al._ (2013)
$IC\equiv-2\log\hat{\mathcal{L}}_{\rm BBN}+2k-1$, $k$ being the number of BSM
parameters and accounting for the CMB information as an extra constraint in
the fit. Then, we evaluate the $IC$ difference with respect to the SM
prediction of the primordial light abundances within a given approach: $\Delta
IC\sim\mathcal{O}(1)$ ($\sim\mathcal{O}(10)$) provides positive (strong)
support in favor of NP beyond $\Lambda$CDM according to the canonical scales
of evidence Kass and Raftery (1995).

##  IV Results

Scenario | Approach | $\ \rm Y_{P}\times 10\ $ | $\ {\rm D/H}\times 10^{5}\ $ | $\ \Delta N_{\rm eff}\ $ | $\ \xi_{\nu}\ $ | $\ \eta_{B}\times 10^{10}\ $ | $\ \Delta IC\ $
---|---|---|---|---|---|---|---
SM prediction | PRIMAT driven | 2.4715(14) | 2.439(36) | – | – | 6.137(38) | –
NACRE II driven | 2.4706(16) | 2.51(10) | – | – | 6.137(38) | –
$\Delta N_{\rm eff}$ BSM fit | PRIMAT driven | 2.472(11) | 2.471(44) | 0.02(18) | – | 6.096(64) | 2
NACRE II driven | 2.453(14) | 2.46(11) | -0.28(23) | – | 6.088(67) | 0
($\Delta N_{\rm eff},\xi_{\nu}$) BSM fit | PRIMAT driven | 2.393(38) | 2.475(44) | 0.27(24) | 0.039(18) | 6.120(67) | 8
NACRE II driven | 2.383(41) | 2.47(11) | -0.01(27) | 0.036(19) | 6.114(65) | 5
Table 1: 68% probability interval for the posterior distribution of the main
observables and parameters in the scenarios considered in this work. For the
BSM fits, improvement with respect to the SM is given by $\Delta IC>0$, see
text for more details.

In Figure 1 we report the main result of our study: the 68% and 90%
probability region for the primordial lepton asymmetry $\xi_{\nu}$ and the
extra-relativistic degrees of freedom $\Delta N_{\rm eff}$ as determined by
TScosmo, Eq. (6), corresponding to the two approaches to thermonuclear rates
described in the previous section. From the $\Lambda$CDM limit highlighted in
the same figure, we can conclude that a BSM fit to a dataset that includes the
newly measured EMPGS by Subaru Matsumoto  _et al._ (2022) favors at present a
non-vanishing asymmetry in the neutrino sector.

In Figure 1 we also observe that, dependent on the approach to nuclear
uncertainties, a shift of $N_{\rm eff}$ of $\mathcal{O}(1)$ can be
simultaneously favored by current cosmological data. Note that the size of the
shift in the number of relativistic degrees of freedom can be interpreted
within the 68% probability region as the result of a large neutrino asymmetry
in the muon-tau sector in case flavor equilibration has not been fully
realized.

To further investigate the different outcome from each approach, we show in
Figure 2 the posterior probability density function (p.d.f.) for the BBN
observables ${\rm Y_{P}}$ and ${\rm D/H}$. We report both the result from the
BSM fit varying both $\xi_{\nu}$ and $N_{\rm eff}$ as well as the one from the
SM prediction, obtained fixing the BSM parameters to 0 and replacing the CMB
likelihood with the Gaussian prior: $\Omega_{B}h^{2}=0.02242\pm 0.00014$, from
the $\Lambda$CDM Planck analysis (TTTEEE + low-$\ell$ \+ BAO + lensing)
Aghanim  _et al._ (2020, 2019a). In the same figure, we also highlight with
vertical dark green bands the measurements adopted in our BBN analysis via Eq.
(8), and report the PDG 2021 value ${\rm Y_{P}}=0.245(3)$ Zyla  _et al._
(2020), in optimal agreement with the analysis of Ref. Xu  _et al._ (2013)
that comprises the set studied also in Matsumoto  _et al._ (2022) without the
new EMPGs from Subaru.

Figure 2 neatly highlights two tensions in the limit where BSM physics is not
accounted for:

  * •

A discrepancy at the 3$\sigma$ level between the SM prediction of ${\rm
Y_{P}}$ and the newly inferred helium-4 mass-fraction value, regardless of the
approach taken for the thermonuclear reactions; the tension is fully driven by
the new measurement delivered by Ref. Matsumoto  _et al._ (2022), while the
overall significance also depends on the precision obtained for the inference
of the cosmological baryon abundance within $\Lambda$CDM;

  * •

A tension of about 2$\sigma$ significance between the SM prediction of D/H and
the PDG 2021 recommended measurement Zyla  _et al._ (2020) when the PRIMAT
driven approach is taken for the analysis of the key thermonuclear reactions
involved, in line with recent discussions in the literature Pitrou  _et al._
(2021b).

From Figure 2, it is clear that a shift of $\Delta N_{\rm eff}$ is required
together with $\xi_{\nu}\neq 0$ only when the PRIMAT driven approach is
considered, in order to address the discrepancy consequently present in the
fit in relation to the observed primordial deuterium abundance. In the same
figure it is also evident how the PDG 2021 recommended measurement of the
helium-4 mass fraction is in perfect agreement with the SM prediction, and our
inference for a nonzero degeneracy parameter $\xi_{\nu}$ is the consequence of
adopting the new ${\rm Y_{P}}$ measurement Matsumoto  _et al._ (2022).

We report in Table 1 the $68\%$ probability interval for the scenarios
discussed so far as well as the one for the BSM fit where only $\Delta N_{\rm
eff}$ is considered. Looking at the $\Delta IC$ values, we conclude that a
joint analysis of BBN + CMB data provides mild to strong evidence for a
scenario with non-vanishing lepton asymmetry. Moreover, within the NACRE II
approach no notable support from data is found for the presence of extra
relativistic degrees of freedom in the Early Universe, whereas a scenario
where only $\Delta N_{\rm eff}$ is varied may be slightly preferred by data
over the SM in the case of the PRIMAT driven approach, partially ameliorating
a potential deuterium anomaly.

##  V Discussion and outlook

We find that the new measurement of primeval helium including 10 additional
EMPGs by the Subaru Survey leads to a tantalizing hint for a large total
lepton asymmetry of primordial origin. From Eq. (1) today’s value
$\eta_{L}\gg\eta_{B}$, and ranges from $\sim 0.01$ to $\sim 0.26$, depending
on the details of neutrino oscillations in the Early Universe as well as on
the initial value for the neutrino degeneracy parameters. Such a finding would
require a cosmology dramatically different from the one inferred extrapolating
the SM, and thus hints toward BSM physics.

While it is possible to imagine many types of NP which could generate such a
situation, there are common factors that any successful explanation must
share. At temperatures above the scale of electroweak symmetry restoration,
electroweak sphalerons equilibrate $B+L$ such that the final total lepton and
baryon asymmetries differ by a ${\cal O}(1)$ factor Khlebnikov and
Shaposhnikov (1988); Harvey and Turner (1990). Thus, for a difference of
orders of magnitude between $\eta_{L}$ and $\eta_{B}$ to persist, it must
either be generated after the sphalerons become inactive (in the SM, at the
electroweak phase transition around temperatures of order 100 GeV) or the
individual flavor asymmetries must be distributed such that the net $L$ is
much smaller than the individual asymmetries March-Russell  _et al._ (1999).
The latter scenario would point to flavor-dependent NP in the lepton sector,
with possible interesting implications for the smallness of $\eta_{B}$ as
well, see for instance Dreiner and Ross (1993); March-Russell  _et al._
(1999); Shu  _et al._ (2007); Gu (2010).

Because equilibration of neutrino species depends both on imprecisely
determined mixing parameters and the assumed initial asymmetry in each flavor
Froustey and Pitrou (2022), mapping the inferred neutrino asymmetries during
BBN into the space of consistent initial conditions at some earlier time is an
interesting inverse problem; it requires also assumptions on the
interpretation of the inference carried out here for $\Delta N_{\rm eff}$, and
it is beyond the scope of this work. Several examples of theories capable of
generating a sufficiently large and persistent lepton-flavored neutrino
asymmetry via variations of the Affleck-Dine mechanism Affleck and Dine (1985)
exist in the literature Casas  _et al._ (1999); McDonald (2000); March-Russell
_et al._ (1999); Kawasaki and Murai (2022).

Looking forward, our study highlights the combined power of BBN and CMB data
to provide a key glimpse of an early epoch of the Universe. Making the most of
its opportunity requires more work, to improve the treatment of nuclear
physics inputs and measurements of the primordial abundances. If the current
indications of a large total lepton asymmetry $\eta_{L}$ survives further
scrutiny, it provides important clues towards our understanding of the
Universe in its early stages, offering at the same time a quite privileged
view on BSM physics.

###### Acknowledgements.

We thank Federico Bianchini, Kevork Abazajian, Manoj Kaplinghat, Rouven Essig
and Peizhi Du for discussion. The work of M.V. is supported by the Simons
Foundation under the Simons Bridge for Postdoctoral Fellowships at SCGP and
YITP, award number 815892. The work of TMPT is supported by the US National
Science Foundation under grant no. PHY-1915005.

## References

  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Aiola  _et al._ (2020) S. Aiola _et al._ (ACT), JCAP 12, 047 (2020), arXiv:2007.07288 [astro-ph.CO] .
  * Dutcher  _et al._ (2021) D. Dutcher _et al._ (SPT-3G), Phys. Rev. D 104, 022003 (2021), arXiv:2101.01684 [astro-ph.CO] .
  * Baumann  _et al._ (2016) D. Baumann, D. Green, J. Meyers, and B. Wallisch, JCAP 01, 007 (2016), arXiv:1508.06342 [astro-ph.CO] .
  * Komatsu (2022) E. Komatsu, (2022), arXiv:2202.13919 [astro-ph.CO] .
  * Chang  _et al._ (2022) C. L. Chang _et al._ , (2022), arXiv:2203.07638 [astro-ph.CO] .
  * Sarkar (1996) S. Sarkar, Rept. Prog. Phys. 59, 1493 (1996), arXiv:hep-ph/9602260 .
  * Steigman (2007) G. Steigman, Ann. Rev. Nucl. Part. Sci. 57, 463 (2007), arXiv:0712.1100 [astro-ph] .
  * Iocco  _et al._ (2009) F. Iocco, G. Mangano, G. Miele, O. Pisanti, and P. D. Serpico, Phys. Rept. 472, 1 (2009), arXiv:0809.0631 [astro-ph] .
  * Pospelov and Pradler (2010) M. Pospelov and J. Pradler, Ann. Rev. Nucl. Part. Sci. 60, 539 (2010), arXiv:1011.1054 [hep-ph] .
  * Serpico and Raffelt (2004) P. D. Serpico and G. G. Raffelt, Phys. Rev. D 70, 043526 (2004), arXiv:astro-ph/0403417 .
  * Boehm  _et al._ (2013) C. Boehm, M. J. Dolan, and C. McCabe, JCAP 08, 041 (2013), arXiv:1303.6270 [hep-ph] .
  * Sabti  _et al._ (2020) N. Sabti, J. Alvey, M. Escudero, M. Fairbairn, and D. Blas, JCAP 01, 004 (2020), arXiv:1910.01649 [hep-ph] .
  * Giovanetti  _et al._ (2021) C. Giovanetti, M. Lisanti, H. Liu, and J. T. Ruderman, (2021), arXiv:2109.03246 [hep-ph] .
  * Cooke  _et al._ (2016) R. J. Cooke, M. Pettini, K. M. Nollett, and R. Jorgenson, Astrophys. J. 830, 148 (2016), arXiv:1607.03900 [astro-ph.CO] .
  * Riemer-Sørensen  _et al._ (2017) S. Riemer-Sørensen, S. Kotuš, J. K. Webb, K. Ali, V. Dumont, M. T. Murphy, and R. F. Carswell, Mon. Not. Roy. Astron. Soc. 468, 3239 (2017), arXiv:1703.06656 [astro-ph.CO] .
  * Cooke  _et al._ (2018) R. J. Cooke, M. Pettini, and C. C. Steidel, Astrophys. J. 855, 102 (2018), arXiv:1710.11129 [astro-ph.CO] .
  * Zyla  _et al._ (2020) P. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020), and 2021 update.
  * Pitrou  _et al._ (2021a) C. Pitrou, A. Coc, J.-P. Uzan, and E. Vangioni, Mon. Not. Roy. Astron. Soc. 502, 2474 (2021a), arXiv:2011.11320 [astro-ph.CO] .
  * Pisanti  _et al._ (2021) O. Pisanti, G. Mangano, G. Miele, and P. Mazzella, JCAP 04, 020 (2021), arXiv:2011.11537 [astro-ph.CO] .
  * Yeh  _et al._ (2021) T.-H. Yeh, K. A. Olive, and B. D. Fields, JCAP 03, 046 (2021), arXiv:2011.13874 [astro-ph.CO] .
  * Pitrou  _et al._ (2021b) C. Pitrou, A. Coc, J.-P. Uzan, and E. Vangioni, Nature Rev. Phys. 3, 231 (2021b), arXiv:2104.11148 [astro-ph.CO] .
  * Mossa  _et al._ (2020) V. Mossa _et al._ , Nature 587, 210 (2020).
  * Sabti  _et al._ (2021) N. Sabti, J. Alvey, M. Escudero, M. Fairbairn, and D. Blas, JCAP 08, A01 (2021), arXiv:2107.11232 [hep-ph] .
  * Matsumoto  _et al._ (2022) A. Matsumoto _et al._ , (2022), arXiv:2203.09617 [astro-ph.CO] .
  * Hsyu  _et al._ (2020) T. Hsyu, R. J. Cooke, J. X. Prochaska, and M. Bolte, Astrophys. J. 896, 77 (2020), arXiv:2005.12290 [astro-ph.GA] .
  * Aver  _et al._ (2015) E. Aver, K. A. Olive, and E. D. Skillman, JCAP 07, 011 (2015), arXiv:1503.08146 [astro-ph.CO] .
  * (28) A.-K. Burns, T. M.P. Tait, and M. Valli, “PRyMordial: The first minutes beyond the Standard Model in a few seconds,” in preparation.
  * Kolb and Turner (1990) E. W. Kolb and M. S. Turner,  _The Early Universe_, Vol. 69 (1990).
  * Rubakov and Gorbunov (2017) V. A. Rubakov and D. S. Gorbunov,  _Introduction to the Theory of the Early Universe: Hot big bang theory_ (World Scientific, Singapore, 2017).
  * Simha and Steigman (2008) V. Simha and G. Steigman, JCAP 08, 011 (2008), arXiv:0806.0179 [hep-ph] .
  * Akita and Yamaguchi (2020) K. Akita and M. Yamaguchi, JCAP 08, 012 (2020), arXiv:2005.07047 [hep-ph] .
  * Froustey  _et al._ (2020) J. Froustey, C. Pitrou, and M. C. Volpe, JCAP 12, 015 (2020), arXiv:2008.01074 [hep-ph] .
  * Bennett  _et al._ (2021) J. J. Bennett, G. Buldgen, P. F. De Salas, M. Drewes, S. Gariazzo, S. Pastor, and Y. Y. Y. Wong, JCAP 04, 073 (2021), arXiv:2012.02726 [hep-ph] .
  * Grohs  _et al._ (2016) E. Grohs, G. M. Fuller, C. T. Kishimoto, M. W. Paris, and A. Vlasenko, Phys. Rev. D 93, 083522 (2016), arXiv:1512.02205 [astro-ph.CO] .
  * Escudero Abenza (2020) M. Escudero Abenza, JCAP 05, 048 (2020), arXiv:2001.04466 [hep-ph] .
  * Kinney and Riotto (1999) W. H. Kinney and A. Riotto, Phys. Rev. Lett. 83, 3366 (1999), arXiv:hep-ph/9903459 .
  * Lesgourgues and Pastor (1999) J. Lesgourgues and S. Pastor, Phys. Rev. D 60, 103521 (1999), arXiv:hep-ph/9904411 .
  * Oldengott and Schwarz (2017) I. M. Oldengott and D. J. Schwarz, EPL 119, 29001 (2017), arXiv:1706.01705 [astro-ph.CO] .
  * Abazajian  _et al._ (2002) K. N. Abazajian, J. F. Beacom, and N. F. Bell, Phys. Rev. D 66, 013008 (2002), arXiv:astro-ph/0203442 .
  * Castorina  _et al._ (2012) E. Castorina, U. Franca, M. Lattanzi, J. Lesgourgues, G. Mangano, A. Melchiorri, and S. Pastor, Phys. Rev. D 86, 023517 (2012), arXiv:1204.2510 [astro-ph.CO] .
  * Aghanim  _et al._ (2019a) N. Aghanim _et al._ (Planck), https://wiki.cosmos.esa.int/planck-legacy-archive/images/4/43/Baseline_params_table_2018_68pc_v2.pdf (2019a).
  * Barenboim  _et al._ (2017) G. Barenboim, W. H. Kinney, and W.-I. Park, Eur. Phys. J. C 77, 590 (2017), arXiv:1609.03200 [astro-ph.CO] .
  * Dolgov  _et al._ (2002) A. D. Dolgov, S. H. Hansen, S. Pastor, S. T. Petcov, G. G. Raffelt, and D. V. Semikoz, Nucl. Phys. B 632, 363 (2002), arXiv:hep-ph/0201287 .
  * Pastor  _et al._ (2009) S. Pastor, T. Pinto, and G. G. Raffelt, Phys. Rev. Lett. 102, 241302 (2009), arXiv:0808.3137 [astro-ph] .
  * Fields  _et al._ (2020) B. D. Fields, K. A. Olive, T.-H. Yeh, and C. Young, JCAP 03, 010 (2020), [Erratum: JCAP 11, E02 (2020)], arXiv:1912.01132 [astro-ph.CO] .
  * Pitrou  _et al._ (2018) C. Pitrou, A. Coc, J.-P. Uzan, and E. Vangioni, Phys. Rept. 754, 1 (2018), arXiv:1801.08023 [astro-ph.CO] .
  * Froustey and Pitrou (2022) J. Froustey and C. Pitrou, JCAP 03, 065 (2022), arXiv:2110.11889 [hep-ph] .
  * Escudero (2019) M. Escudero, JCAP 02, 007 (2019), arXiv:1812.05605 [hep-ph] .
  * Chu  _et al._ (2022) X. Chu, J.-L. Kuo, and J. Pradler, (2022), arXiv:2205.05714 [hep-ph] .
  * Bennett  _et al._ (2020) J. J. Bennett, G. Buldgen, M. Drewes, and Y. Y. Y. Wong, JCAP 03, 003 (2020), [Addendum: JCAP 03, A01 (2021)], arXiv:1911.04504 [hep-ph] .
  * Ade  _et al._ (2019) P. Ade _et al._ (Simons Observatory), JCAP 02, 056 (2019), arXiv:1808.07445 [astro-ph.CO] .
  * Abazajian  _et al._ (2022) K. Abazajian _et al._ (CMB-S4), in  _2022 Snowmass Summer Study_ (2022) arXiv:2203.08024 [astro-ph.CO] .
  * Aiola  _et al._ (2022) S. Aiola _et al._ (CMB-HD), (2022), arXiv:2203.05728 [astro-ph.CO] .
  * Sirlin (1967) A. Sirlin, Phys. Rev. 164, 1767 (1967).
  * Brown and Sawyer (2001) L. S. Brown and R. F. Sawyer, Phys. Rev. D 63, 083503 (2001), arXiv:astro-ph/0006370 .
  * Czarnecki  _et al._ (2004) A. Czarnecki, W. J. Marciano, and A. Sirlin, Phys. Rev. D 70, 093006 (2004), arXiv:hep-ph/0406324 .
  * Ivanov  _et al._ (2017) A. N. Ivanov, R. Höllwieser, N. I. Troitskaya, M. Wellenzohn, and Y. A. Berdnikov, Phys. Rev. D 95, 033007 (2017), arXiv:1701.04613 [hep-ph] .
  * Lopez  _et al._ (1997) R. E. Lopez, M. S. Turner, and G. Gyuk, Phys. Rev. D 56, 3191 (1997), arXiv:astro-ph/9703065 .
  * Iliadis and Coc (2020) C. Iliadis and A. Coc, Astrophys. J. 901, 127 (2020), arXiv:2008.12200 [astro-ph.CO] .
  * Serpico  _et al._ (2004) P. D. Serpico, S. Esposito, F. Iocco, G. Mangano, G. Miele, and O. Pisanti, JCAP 12, 010 (2004), arXiv:astro-ph/0408076 .
  * Ando  _et al._ (2006) S. Ando, R. H. Cyburt, S. W. Hong, and C. H. Hyun, Phys. Rev. C 74, 025809 (2006), arXiv:nucl-th/0511074 .
  * Descouvemont  _et al._ (2004) P. Descouvemont, A. Adahchour, C. Angulo, A. Coc, and E. Vangioni-Flam, Atom. Data Nucl. Data Tabl. 88, 203 (2004), arXiv:astro-ph/0407101 .
  * Iliadis  _et al._ (2016) C. Iliadis, K. Anderson, A. Coc, F. Timmes, and S. Starrfield, Astrophys. J. 831, 107 (2016), arXiv:1608.05853 [astro-ph.SR] .
  * Iñesta Gómez  _et al._ (2017) A. Iñesta Gómez, C. Iliadis, and A. Coc, Astrophys. J. 849, 134 (2017), arXiv:1710.01647 [astro-ph.IM] .
  * de Souza  _et al._ (2019a) R. S. de Souza, C. Iliadis, and A. Coc, Astrophys. J. 872, 75 (2019a), arXiv:1809.06966 [astro-ph.IM] .
  * de Souza  _et al._ (2019b) R. S. de Souza, S. R. Boston, A. Coc, and C. Iliadis, Phys. Rev. C 99, 014619 (2019b), arXiv:1901.04857 [nucl-th] .
  * Moscoso  _et al._ (2021) J. Moscoso, R. S. de Souza, A. Coc, and C. Iliadis, Astrophys. J. 923, 49 (2021), arXiv:2109.00049 [astro-ph.CO] .
  * Xu  _et al._ (2013) Y. Xu, K. Takahashi, S. Goriely, M. Arnould, M. Ohta, and H. Utsunomiya, Nucl. Phys. A 918, 61 (2013), arXiv:1310.7099 [nucl-th] .
  * Aghanim  _et al._ (2019b) N. Aghanim _et al._ (Planck), https://wiki.cosmos.esa.int/planck-legacy-archive/index.php/Cosmological_Parameters (2019b).
  * Coc  _et al._ (2014) A. Coc, J.-P. Uzan, and E. Vangioni, JCAP 10, 050 (2014), arXiv:1403.6694 [astro-ph.CO] .
  * Foreman-Mackey  _et al._ (2013) D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, PASP 125, 306 (2013), arXiv:1202.3665 [astro-ph.IM] .
  * Akaike (1974) H. Akaike, IEEE Transactions on Automatic Control 19, 716 (1974).
  * Gelman  _et al._ (2013) A. Gelman, J. Hwang, and A. Vehtari, arXiv e-prints , arXiv:1307.5928 (2013), arXiv:1307.5928 [stat.ME] .
  * Kass and Raftery (1995) R. E. Kass and A. E. Raftery, Journal of the American Statistical Association 90, 773 (1995), http://dx.doi.org/10.1080/01621459.1995.10476572.
  * Khlebnikov and Shaposhnikov (1988) S. Y. Khlebnikov and M. E. Shaposhnikov, Nucl. Phys. B 308, 885 (1988).
  * Harvey and Turner (1990) J. A. Harvey and M. S. Turner, Phys. Rev. D 42, 3344 (1990).
  * March-Russell  _et al._ (1999) J. March-Russell, H. Murayama, and A. Riotto, JHEP 11, 015 (1999), arXiv:hep-ph/9908396 .
  * Dreiner and Ross (1993) H. K. Dreiner and G. G. Ross, Nucl. Phys. B 410, 188 (1993), arXiv:hep-ph/9207221 .
  * Shu  _et al._ (2007) J. Shu, T. M. P. Tait, and C. E. M. Wagner, Phys. Rev. D 75, 063510 (2007), arXiv:hep-ph/0610375 .
  * Gu (2010) P.-H. Gu, Phys. Rev. D 82, 093009 (2010), arXiv:1005.1632 [hep-ph] .
  * Affleck and Dine (1985) I. Affleck and M. Dine, Nucl. Phys. B 249, 361 (1985).
  * Casas  _et al._ (1999) A. Casas, W. Y. Cheng, and G. Gelmini, Nucl. Phys. B 538, 297 (1999), arXiv:hep-ph/9709289 .
  * McDonald (2000) J. McDonald, Phys. Rev. Lett. 84, 4798 (2000), arXiv:hep-ph/9908300 .
  * Kawasaki and Murai (2022) M. Kawasaki and K. Murai, (2022), arXiv:2203.09713 [hep-ph] .