# The temperature-dependence of carrier mobility is not a reliable indicator
of the dominant scattering mechanism

Alex M. Ganose<EMAIL_ADDRESS>Energy Technologies Area, Lawrence
Berkeley National Laboratory, Berkeley, California 94720, USA  Department of
Chemistry, Molecular Sciences Research Hub, White City Campus, Imperial
College London, Wood Lane, London, UK  Junsoo Park  Energy Technologies Area,
Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA
Anubhav Jain<EMAIL_ADDRESS>Energy Technologies Area, Lawrence Berkeley
National Laboratory, Berkeley, California 94720, USA

##  S1 DFPT+Wannier dataset

In the main manuscript, we present a summary of the temperature dependence of
mobility calculated from density functional perturbation theory combined with
Wannier interpolation (DFPT+Wannier). Our dataset includes 47 compounds and
covers both $n$ and $p$-type doping, and bulk and monolayer structures. The
data was compiled from studies in which DFPT+Wannier calculations of mobility
were performed at multiple temperatures. In order to indicate the dominant
scattering mechanism, we only consider results that include a breakdown of
different scattering types, such as the energy or k-dependence of the
scattering rates separated into different mechanisms. At the time of writing,
our data set includes all DFPT+Wannier calculations that satisfy the above
selection criteria. Mobility data was extracted using WebPlotDigitizer 1. The
temperature ($T$) dependence of mobility ($\mu$) was calculated by fitting the
mobility to

$\mu(T)=T^{a}+b,$ (S1)

where $a$ and $b$ are fitted parameters. The temperature dependence of
mobility for each materials in the dataset are listed in Table.
LABEL:tab:dfpt-all. Furthermore, in Fig. S1, we plot the full temperature
dependent results and include the line of best fit calculated from Eq. S1.

Table S1: Temperature dependence of mobility for all materials calculated by density functional perturbation theory combined with Wannier interpolation (DFPT+Wannier). We indicate whether the doping is $n$ or $p$-type, the dominant scattering mechanism, and the effective mass of majority carriers ($m^{*}$). The scattering mechanisms identified are deformation potential scattering (DP), polar optical phonon scattering (PO), piezoelectric scattering (PI), and polaronic effects (PL). The unit of effective mass is the electron rest mass ($m_{0}$). Formula | Doping | Topology |  $T$ dependence | Mechanism |  $m^{*}$ ($m_{0}$) | Ref.
---|---|---|---|---|---|---
$\alpha$-Te | $p$ | monolayer | -1.15 | DP | 0.24 | 2
$\alpha$-Te | $n$ | monolayer | -1.47 | DP | 0.10 | 2
1T”-MoSe2 | $p$ | monolayer | -1.98 | DP | 0.54 | 3
1T”-MoSe2 | $n$ | monolayer | -1.52 | DP | 0.94 | 3
3C-SiC | $p$ | bulk | -1.90 | PO, DP | 0.45 | 4
3C-SiC | $n$ | bulk | -2.04 | PO, DP | 0.31 | 4
BC3 | $p$ | monolayer | -1.90 | DP | 0.50 | 5
BC3 | $n$ | monolayer | -1.89 | DP | 0.30 | 5
Ba2BiAu | $n$ | bulk | -1.42 | PO | 0.16 | 6
C3N | $n$ | monolayer | -2.03 | DP | 0.45 | 5
C3N | $p$ | monolayer | -2.50 | DP | 0.24 | 5
CH3NH3PbI3 | $n$ | bulk | -1.49 | PO | 0.21 | 7
CsPbI3 | $n$ | bulk | -0.93 | PO | 0.17 | 7
Ga2O3 | $n$ | bulk | -1.75 | PO, DP | 0.25 | 8
GaAs | $n$ | bulk | -1.70 | PO | 0.06 | 9
GaN | $n$ | bulk | -1.91 | PO, PI | 0.20 | 7
GeO2 | $n$ | bulk | -2.73 | PO | 0.29 | 10
GeO2 | $p$ | bulk | -2.96 | PO | 1.47 | 10
InSe | $p$ | monolayer | -0.41 | DP | 3.40 | 11
InSe | $p$ | bulk | -1.55 | PO | 2.60 | 11
InSe | $n$ | monolayer | -1.86 | PO | 0.16 | 11
InSe | $n$ | bulk | -2.13 | PO | 0.10 | 11
MoS2 | $n$ | monolayer | -1.59 | PO | 0.49 | 12
PbTe | $n$ | bulk | -1.30 | PO | 0.05 | 13
Rb3AuO | $n$ | bulk | -1.16 | PO | 0.27 | 14
Rb3AuO | $p$ | bulk | -1.17 | PO | 0.23 | 14
Si | $n$ | bulk | -2.20 | DP | 0.26 | 15
Si | $p$ | bulk | -2.37 | DP | 0.22 | 15
SnSe | $p$ | bulk | -1.86 | PO | 0.22 | 16
Sr2BiAu | $n$ | bulk | -2.20 | PO, DP | 0.21 | 17
Sr2SbAu | $n$ | bulk | -2.05 | PO, DP | 0.21 | 17
SrTiO3 | $n$ | bulk | -3.09 | PO, PL | 0.75 | 18
Tl2O | $n$ | monolayer | -0.59 | PO | 0.15 | 19
Tl2O | $p$ | monolayer | -0.87 | PO, DP | 0.28 | 19
ZnS | $p$ | bulk | -1.26 | PO, DP | 0.55 | 20
ZnSe | $p$ | bulk | -1.05 | PO, DP | 0.75 | 20
antimonene | $n$ | monolayer | -1.33 | DP | 0.21 | 21
antimonene | $p$ | monolayer | -1.60 | DP | 0.46 | 21
arsene | $p$ | monolayer | -1.40 | DP | 0.53 | 21
arsene | $n$ | monolayer | -1.37 | DP | 0.23 | 21
germanene | $p$ | monolayer | -3.14 | DP | 0.02 | 21
germanene | $n$ | monolayer | -2.85 | DP | 0.02 | 21
h-BeO | $n$ | monolayer | -0.70 | PO | 0.79 | 22
phosphorene | $p$ | monolayer | -3.01 | DP | 0.85 | 21
phosphorene | $n$ | monolayer | -3.32 | DP | 0.44 | 21
silicene | $n$ | monolayer | -3.18 | DP | 0.01 | 21
stannene | $n$ | monolayer | -2.93 | DP | 0.03 | 21

Figure S1: Temperature against mobility (pink line) for all materials in the
density functional perturbation theory combined with Wannier interpolation
(DFPT+Wanner) dataset. The line of best fit, calculated from Eq. S1, is shown
in grey.

##  S2 Temperature dependence in a single parabolic band

Figure S2: Increasing the temperature range or carrier concentration weakens
(makes less negative) the temperature dependence of mobility. Calculations
performed in a single parabolic band using polar optical phonon scattering
(PO), acoustic deformation potential scattering (AD), and ionized impurity
scattering (II). a,b The temperature dependence of mobility against carrier
concentration at a low temperatures
($300\text{\,}\mathrm{K}500\text{\,}\mathrm{K}$) and b high temperatures
($800\text{\,}\mathrm{K}1000\text{\,}\mathrm{K}$). c,d The temperature
dependence of mobility against the temperature range in which mobility is
measured at c low carrier concentrations
(${10}^{13}\text{\,}{\mathrm{cm}}^{-3}$) and d high carrier concentrations
(${10}^{18}\text{\,}{\mathrm{cm}}^{-3}$).

##  S3 Impact of band structure on temperature dependence

In the main text, we investigate the impact of non-isotropic and non-parabolic
bands on the $T$ dependence of mobility. To include non-parabolicity, we
generated a series of band structures according to the Kane quasi-linear
dispersion relation 23

$\frac{\hbar^{2}k}{2m^{*}}=\varepsilon(1+\alpha\varepsilon),$ (S2)

where $\hbar$ is the reduced Planck constant, $\varepsilon$ is the band
energy, $k=\left|\mathbf{k}\right|$ is the distance in reciprocal space of k
from the band edge (set to the $\Gamma$-point $[0,0,0]$), and $m^{*}$ is the
effective mass of charge carriers. A larger value of $\alpha$ (Kane parameter)
indicates increased linearity of the band away from the band minimum. Due to
periodicity, the eigenvalues of Bloch functions exhibit a derivative of zero
at the Brillouin zone boundary. To add this physical behaviour and ensure that
the band structures could be smoothly Fourier interpolated (as required by
amset) we added a bump function to the dispersion generated by Eq. S2
according to the procedure outlined in Ref. 24. The bump function becomes
active at k-points greater than $80\text{\,}\mathrm{\char 37\relax}$ of the
distance from the zone center to the zone boundary. Even at the largest
carrier concentrations used in our tests
(${10}^{18}\text{\,}{\mathrm{cm}}^{-3}$), the Fermi–Dirac occupation in the
bump region is extremely small. Accordingly, the transport results will be
entirely dictated by the quasi-linear dispersion region around the band edge.
Calculations included polar optical phonon, acoustic deformation potential,
and ionized impurity scattering and were performed using the amset package 25.
Due to lack of wave function information, the wave function overlaps in the
scattering matrix elements were set to unity (indicting complete overlap).
This is a reasonable approximation for a single near-parabolic band and is
consistent with the approximations made in the simplified models of transport
historically used to determine the temperature dependence of mobility. We
considered both low ($300\text{\,}\mathrm{K}500\text{\,}\mathrm{K}$) and high
($800\text{\,}\mathrm{K}1000\text{\,}\mathrm{K}$) temperature regimes, and low
(${10}^{13}\text{\,}{\mathrm{cm}}^{-3}$) and high
(${10}^{18}\text{\,}{\mathrm{cm}}^{-3}$) carrier concentrations. The full
results for all temperature and doping combinations are presented in Fig. S3.

Figure S3: Effect of Kane non-parabolicity parameter ($\alpha$) on the
temperature ($T$) dependence of mobility calculated using polar optical phonon
scattering (PO), acoustic deformation potential scattering (AD), and ionized
impurity scattering (II). We have considered all combinations of low
($300\text{\,}\mathrm{K}500\text{\,}\mathrm{K}$) and high
($800\text{\,}\mathrm{K}1000\text{\,}\mathrm{K}$) temperature regimes, and low
(${10}^{13}\text{\,}{\mathrm{cm}}^{-3}$) and high
(${10}^{18}\text{\,}{\mathrm{cm}}^{-3}$) carrier concentrations.

We have also investigated the impact of band anisotropy by constructing a
series of band structures with varying effective masses in each principal
direction. We enforce a simple cubic symmetry in which the effective masses
along $\Gamma$–M and $\Gamma$–R are held the same while the mass along
$\Gamma$–X is increased from $0.7\text{\,}\mathrm{\text{ $m_{\textup{ e }}$
}}$ to $30\text{\,}\mathrm{\text{ $m_{\textup{ e }}$ }}$. These anisotropic
band structures are generated by smoothly inverting upward paraboloids into
downward paraboloids after the methods of 26, whereby the derivative
continuity is enforced at the point of inflection. By doing so, we ensure that
the entire band structure is parabolic while still bounded by the Brillouin
zone, crossing the zone boundary orthogonally. We also set the inflection
point to be the halfway point to the zone boundary in each direction. This
ensures that the band structure retains the same directional effective masses
after inversion, save for the sign, such that the entire band structure is
described by three directional effective masses. Scattering calculations
included polar optical phonon, acoustic deformation potential, and ionized
impurity scattering and were performed using the amset package 25. Since such
arbitrary band structures lack wave function information, the wave function
overlaps in the scattering matrix elements were set to unity (indicating
complete overlap). We considered both low
($300\text{\,}\mathrm{K}500\text{\,}\mathrm{K}$) and high
($800\text{\,}\mathrm{K}1000\text{\,}\mathrm{K}$) temperature regimes, and low
(${10}^{13}\text{\,}{\mathrm{cm}}^{-3}$) and high
(${10}^{18}\text{\,}{\mathrm{cm}}^{-3}$) carrier concentrations. The full
results for all temperature and doping combinations are presented in Fig. S4

Figure S4: Effect of a heavy anisotropic effective mass ($m_{\mathrm{h}}^{*}$)
on the temperature ($T$) dependence of mobility calculated using polar optical
phonon scattering (PO), acoustic deformation potential scattering (AD), and
ionized impurity scattering (II). We have considered all combinations of low
($300\text{\,}\mathrm{K}500\text{\,}\mathrm{K}$) and high
($800\text{\,}\mathrm{K}1000\text{\,}\mathrm{K}$) temperature regimes, and low
(${10}^{13}\text{\,}{\mathrm{cm}}^{-3}$) and high
(${10}^{18}\text{\,}{\mathrm{cm}}^{-3}$) carrier concentrations.

##  S4 Additional impacts of phonon frequency on the temperature dependence of
mobility

The temperature dependence of mobility is also compounded by the fraction of
occupied electronic states above the phonon emission threshold (Fig. 3b in the
main text). We calculate this fraction as
$\int_{\varepsilon_{\mathrm{po}}}^{\infty}[D(\varepsilon)\partialderivative{f}{\varepsilon}]d\varepsilon/\int_{0}^{\infty}[D(\varepsilon)\partialderivative{f}{\varepsilon}]d\varepsilon$,
where $D$ is the density of states, $\partialderivative{f}{\varepsilon}$ is
the derivative of the Fermi–Dirac distribution and
$\varepsilon_{\mathrm{po}}=2\pi\hbar\omega_{\mathrm{po}}$ is the energy of the
polar optical phonon mode. As revealed in Fig. 3e in the main text, at
$300\text{\,}\mathrm{K}$ and at low phonon frequencies
($4\text{\,}\mathrm{THz}$) over $80\text{\,}\mathrm{\char 37\relax}$ of the
states that contribute to the mobility originate from above the phonon
emission threshold. As the temperature increases to $500\text{\,}\mathrm{K}$,
the fraction of post-emission states only increases by
$10\text{\,}\mathrm{\char 37\relax}$ (a relative increase of
$12.5\text{\,}\mathrm{\char 37\relax}$). In contrast, when the phonon
frequency is large ($20\text{\,}\mathrm{THz}$), the post-emission states
account for only $10\text{\,}\mathrm{\char 37\relax}$ of the states that
contribute to the mobility at $300\text{\,}\mathrm{K}$. However, as the
temperature increases, the fraction of post-emission states increases rapidly
to $\sim$$30\text{\,}\mathrm{\char 37\relax}$ at $500\text{\,}\mathrm{K}$, a
relative increase of $200\text{\,}\mathrm{\char 37\relax}$. This is critical
because the lifetimes of the post-emission states near the emission onset are
approximately an order of magnitude shorter than those of the pre-emission
states. Accordingly, this dramatic increase in the contribution of post
emissions states results in a considerable decrease in the electron mobility
at high temperatures, leading to a quicker decay of mobility with temperature.

At higher, “post-Debye” temperatures
($800\text{\,}\mathrm{K}1000\text{\,}\mathrm{K}$), phonon occupation
essentially increases linearly as $n\propto T$, resulting in more gradual
decay of electron lifetimes with temperature (Fig. 3c,d in the main text).
Furthermore, at high temperatures, the post-emission states account for over
$50\text{\,}\mathrm{\char 37\relax}$ of all states that contribute to mobility
across all phonon frequencies (Fig. 3e in the main text). Together these
combine to drastically reduce the temperature dependence of mobility,
particularly for large $\omega_{\mathrm{po}}$ (Fig. 3a in the main text). For
example, at a phonon frequency of $20\text{\,}\mathrm{THz}$ the mobility
dependence is reduced from $T^{-3.34}$ between
$300\text{\,}\mathrm{K}500\text{\,}\mathrm{K}$ to $T^{-1.81}$ between
$800\text{\,}\mathrm{K}1000\text{\,}\mathrm{K}$.

Thermoelectric materials generally exhibit optimal performance at “high
doping”. As demonstrated in Fig. 3a in the main text, an increase in the
electron concentration from ${10}^{13}\text{\,}{\mathrm{cm}}^{-3}$ to
${10}^{18}\text{\,}{\mathrm{cm}}^{-3}$ results in a weakening of the
temperature dependence across all phonon frequencies. This behaviour can again
be rationalized by considering the fraction of post-emission states that
contribute to the mobility. For phonon frequencies below
$\sim$$12\text{\,}\mathrm{THz}$ almost all the occupied electronic states
($>$$99\text{\,}\mathrm{\char 37\relax}$) are post-emission (Fig. 3f in the
main text). The lifetimes of post-emission states exhibit a relatively weak
temperature dependence (between $\tau\propto T^{-0.98}$–$T^{-0.81}$, see Fig.
3d in the main text) and result in a mobility trend that is always less
negative than $\mu\propto T^{-1.3}$. At larger phonon frequencies, the
fraction of post-emission states decreases to as low as
$30\text{\,}\mathrm{\char 37\relax}$ at $300\text{\,}\mathrm{K}$, and the
mobility picks up additional temperature dependence due to the potential for
thermal activation of these states.

##  S5 Generation of AMSET dataset

In the main manuscript, we present the results of 23,000 mobility calculations
performed using the amset package. amset enables the explicit calculation of
charge carrier scattering rates and the solution of the Boltzmann transport
equation to obtain transport properties. amset has been shown provide
scattering rates and mobility within $\sim$ $10\text{\,}\mathrm{\char
37\relax}$ of DFPT+Wannier when benchmarked on 23 semiconductors 25. This is
highlighted in Fig. S5, where we demonstrate that amset can reproduce the
experimental mobility of GaAs across the full temperature range.

Figure S5: amset reproduces the experimental mobility of $n$-GaAs. Temperature
against mobility for $n$-GaAs at low doping
(${10}^{13}\text{\,}{\mathrm{cm}}^{-3}$). Experimental measurements (black
triangles, Ref. 27) compared against amset (blue line), density functional
perturbation theory combined with Wannier interpolation (DFPT+Wannier, red
line, Ref. 28) and the constant relaxation time approximation (grey line) with
the lifetime set to ${10}^{-14}\text{\,}\mathrm{s}$.

The main inputs for amset are an electronic band structure and common
materials parameters, including the static and high-frequency dielectric
constants, elastic constant, polar phonon frequency, and deformation
potentials. There are four key differences between the amset calculations
performed in this work and the calculations in Ref. 25. Firstly, the wave
function overlaps that enter the scattering matrix elements can be computed
using the first principles wave functions directly or approximated using the
orbital projections, as outlined in Rode 29. As the wave functions were not
available for the 23,000 materials used in this dataset, we have used the
latter approach when calculating the matrix elements. Secondly, in Ref. 25,
the materials parameters were all obtained from first principles calculations.
In this work, we have used a combination of first principles calculations and
machine learning to obtain the materials parameters for each material.
Thirdly, in this work we use only two values for the deformation potential
(one for the valence bands and one for the conduction bands) which are
constant across all k-points, whereas in Ref. 25 the deformation potential is
a tensor that depends on the band index and k-point. Finally, amset has the
capacity to exploit any anisotropy in material parameters by using full
property tensors. However, due to the difficulty of machine learning a full
tensor, in this work we used isotropically averaged materials parameters if
the first principles result was not available. While these approximations will
result in a loss of accuracy, we stress that our goal is only to simulate a
realistic distribution of transport properties rather than predict exact
values for specific materials. Furthermore, while some anisotropy will be lost
by the averaging of tensor properties, we still account for any anisotropy in
the electronic band structure by use of the full k-dependent group velocities
and eigenvalues.

amset calculations were performed on uniform electronic band structures taken
from the Materials Project.30 Band structures were calculated using density
functional theory (DFT) with the Perdew–Burke–Ernzerhof31 (PBE) exchange
correlation functional as implemented in the Vienna ab initio simulation
package (VASP).32, 33 Corrections for the Hubbard U were implemented for
cobalt, chromium, iron, manganese, molybdenum, nickel, vanadium, and tungsten
in oxides using the standard Materials Project values.34 Calculations were
performed in a plane-wave basis set with scalar relativistic psueodpoentials
and with the interactions between core and valence electrons described using
the projector augmented-wave method (PAW).35, 36 The plane-wave energy cutoff
was set to 520 eV. Several versions of the pymatgen37 input sets were been
used to calculate the band structures. The primary difference of between the
inputs — in the context of amset transport calculations — is the density of
the k-point mesh given in number of k-points per reciprocal volume. The
majority of calculations in this dataset ($80\text{\,}\mathrm{\char
37\relax}$) used a k-point density of $100\text{\,}\textbf{k}{\mathchar
45\relax}\mathrm{p}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{s}\mathrm{/}\mathrm{\SIUnitSymbolAngstrom}^{-3}$,
with the remainder using densities of $500\text{\,}\textbf{k}{\mathchar
45\relax}\mathrm{p}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{s}\mathrm{/}\mathrm{\SIUnitSymbolAngstrom}^{-3}$
($8\text{\,}\mathrm{\char 37\relax}$), $1000\text{\,}\textbf{k}{\mathchar
45\relax}\mathrm{p}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{s}\mathrm{/}\mathrm{\SIUnitSymbolAngstrom}^{-3}$
($8\text{\,}\mathrm{\char 37\relax}$), and $2000\text{\,}\textbf{k}{\mathchar
45\relax}\mathrm{p}\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{s}\mathrm{/}\mathrm{\SIUnitSymbolAngstrom}^{-3}$
($4\text{\,}\mathrm{\char 37\relax}$). The band structures were accessed using
the Materials Project application programming interface (API) available in
pymatgen.37 The corresponding wave functions were not available, so we have
used the orbital projections at each band and k-point to approximate the wave
function overlap. We have considered all materials in the Materials Project
database for which the band gap is greater than $0.05\text{\,}\mathrm{eV}$ and
for which there is a uniform band structure available.

###  S5.1 Materials parameter datasets

To obtain the materials parameters needed to compute scattering rates for all
23,000 materials, we have employed a combination of first principles databases
and machine learning. In this section, we outline the generation and cleaning
of the datasets. The Materials Project database contains 13,421 elastic
constant tensors and 8,570 dielectric tensors (and Born effective charges and
$\Gamma$-point phonon frequencies as these are generated in the same
calculation).30, 38, 39 We used the Born effective charges and phonon
frequencies to calculate the polar phonon frequency, needed for polar optical
phonon scattering, according to the methodology outlined in Ref. 25. In
addition, we generated a dataset of 995 valence and conduction band acoustic
deformation potentials. The 995 materials were a subset of the semiconductor
materials for which elastic constant tensors were available. The deformation
potentials were calculated similarly to the methodology outlined in Ref. 25,
however, we only consider two deformation potentials — one for the valence
band maximum and one for the conduction band minimum — which are constant
across all bands and k-points. Finally, we used the dataset of 6,030 band gaps
calculated using the HSE06 exchange–correlation functional prepared by Jie et
al. 40. In Table S2, we summarise the source and number of materials in each
dataset.

Property |  $N$ materials (raw) |  $N$ materials (cleaned) | RMSE | $r^{2}$ | Ref.
---|---|---|---|---|---
$E_{\mathrm{g}}^{\mathrm{HSE06}}$ | 6,030 | 5,658a | $0.97\text{\,}\mathrm{eV}$ | 0.78 | 40
$C$ | 13,421 | 11,160a,b,c | $26.1\text{\,}\mathrm{GPa}$ | 0.92 | 38
$D_{\mathrm{vbm}}$ | 995 | 954a,d | $1.49\text{\,}\mathrm{eV}$ | 0.40 | This work
$D_{\mathrm{cbm}}$ | 995 | 954a,d | $2.01\text{\,}\mathrm{eV}$ | 0.31 | This work
$\epsilon_{\mathrm{s}}$ | 8,570 | 6,507a,e,f,g,h | 3.41 | 0.66 | 39
$\epsilon_{\infty}$ | 8,570 | 6,507a,e,f,g,h | 1.24 | 0.80 | 39
$\omega_{\mathrm{po}}$ | 8,569 | 6,507a,e,f,g,h | $1.68\text{\,}\mathrm{THz}$ | 0.88 | This work
|  |  |  |  |

Table S2: Summary of datasets used for machine learning of materials
properties. Properties include the band gap calculated using the HSE06
exchange–correlation functional ($E_{\mathrm{g}}^{\mathrm{HSE06}}$), elastic
constant ($C$), acoustic deformation potential at the valence band maximum
($D_{\mathrm{vbm}}$) and conduction band minimum ($D_{\mathrm{cbm}}$), static
($\epsilon_{\mathrm{s}}$) and high-frequency ($\epsilon_{\infty}$) dielectric
constants, and the polar optical phonon frequency ($\omega_{\mathrm{po}}$).
The root mean squared error (RMSE) and $r^{2}$ correlation coefficient from
5-fold cross validation for the trained machine learning models are also
provided

a Removed materials with a formation energy or energy above the convex hull
more than 500 meV.

b Removed materials with GVoigt, GReuss, GVRH, KVoigt, KReuss, or
K${}_{\mathrm{VRH}}<1$.

c Removed materials failing
$\mathrm{G}_{\mathrm{Voigt}}\leq\mathrm{G}_{\mathrm{VRH}}\leq\mathrm{G}_{\mathrm{Reuss}}$
or
$\mathrm{K}_{\mathrm{Voigt}}\leq\mathrm{K}_{\mathrm{VRH}}\leq\mathrm{K}_{\mathrm{Reuss}}$.

d Removed materials with $D_{\mathrm{vbm}}$ or $D_{\mathrm{cbm}}$ $>$
$30\text{\,}\mathrm{eV}$.

e Removed materials with a Materials Project band gap less than
$0.1\text{\,}\mathrm{eV}$.

f Removed materials with any imaginary phonon frequencies less (more negative)
than $-1\text{\,}\mathrm{THz}$.

g Removed materials failing
$$300\text{\,}\mathrm{THz}$\geq\epsilon_{\mathrm{s}}\geq\epsilon_{\infty}\geq$1\text{\,}\mathrm{THz}$$.

h Removed materials with $\omega_{\mathrm{po}}<0$.

Before performing machine learning, each dataset was cleaned to remove
unphysical materials properties. Firstly, due to the complexity of machine
learning tensors, we take the tensor average of the dielectric and elastic
tensors. Furthermore, for all datasets, we only keep materials which have a
formation energy and energy above the convex hull less than
$500\text{\,}\mathrm{meV}$ (as taken from the Materials Project). Additional
cleaning steps for each dataset are detailed in the footnotes of Table S2
along with the number of materials before and after data cleaning.

###  S5.2 Machine learning of materials parameters

Machine learning was performed using the automatminer package.41 automatminer
is a tool for automatically creating machine learning pipelines for materials
science datasets and includes automatic featurization with matminer,42 feature
reduction, and an automatic machine learning (AutoML) backend. automatminer
has been benchmarked on a diverse set of materials science challenges and
shows comparable or better performance than bespoke hand tuned models.41 All
models in this work included crystal structures as input and employed the
“express“ automatminer preset. In addition to the features generated by
automatminer, we also included the band gap calculated using PBE (with a
Hubbard U correction if required) taken from the Materials Project.30

Figure S6: Cross validation predictions for machine learning models trained on
the materials parameter datasets. a Static dielectric constant
($\epsilon_{s}$), b polar phonon frequency ($\omega_{\mathrm{po}}$), c HSE06
band gap ($E_{g}^{\mathrm{HSE06}}$), d elastic constant ($C$), and e acoustic
deformation potentials at the valence band maximum ($D_{\mathrm{vbm}}$) and f
conduction band minimum ($D_{\mathrm{cbm}}$). The root mean squared error
(RMSE) of predictions and $r^{2}$ coefficient are given for each model.

We evaluated the models using 5-fold cross validation with the results
presented in Fig. S6. The mean root mean squared error (RMSE) of predictions
and $r^{2}$ correlation coefficient for the 5-fold cross validation are also
summarised in Table S2. Overall, the models trained on larger datasets with
more than 5,000 materials (elastic, dielectric, phonon frequency and HSE06
gap) are considerably more accurate, with $r^{2}$ values between 0.66–0.92.
The acoustic deformation potential datasets — containing only 954 materials —
appear significantly more difficult, with $r^{2}$ values of 0.31 and 0.40 for
the conduction band and valence band potentials, respectively. Despite this,
we note that the deformation potential is temperature independent and
therefore will not affect the temperature-dependence of mobility and, thus,
the conclusions in the main text. Furthermore, our goal is only to simulate a
realistic distribution of parameters rather than the exact parameter for each
material, so strict accuracy is not essential.

###  S5.3 High-throughput AMSET workflow

The machine learning models described in Sec. S5.2 were used to predict
materials parameters for all 131,653 semiconductors in the Materials Project
database. amset calculations were performed on all materials that satisfied
the following criteria: (i) A uniform band structure with orbital projections
was available; (ii) The Materials Project band gap was greater than
$0.05\text{\,}\mathrm{eV}$; (iii) The (predicted) HSE06 band gap was greater
than $0.05\text{\,}\mathrm{eV}$; (iv) The (predicted) static dielectric
constant ($\epsilon_{\mathrm{s}}$) was greater than the (predicted) high-
frequency dielectric constant ($\epsilon_{\infty}$); (v) All predicted
materials properties were greater than zero. In total this amounted to 23,572
materials. We preferentially used materials properties calculated using first
principles methods if they were available and satisfied the criteria outlined
in the footnotes of Table S2. In the case a first principles property was not
available, we used the materials properties predicted from the machine
learning models. 651 materials had all properties calculated from first
principles, the remainder used a mixture of machine learned and first
principles properties. Finally, we applied a scissor correction to shift the
band gap of the Materials Project band structure to match that of the
(predicted) HSE06 gap.

amset calculations were performed in a high-throughput mode using the
fireworks workflow management software.43 amset version `0.4.3` was used
throughout. We used the default settings with the addition of free carrier
screening of the polar optical phonon matrix element enabled by the
`free_carrier_screening` option. Calculations were performed at both $n-$ and
$p$-type doping concentrations and at temperatures ranging from
$200\text{\,}\mathrm{K}1000\text{\,}\mathrm{K}$. To ensure convergence of
transport properties, the workflow performed several amset calculations with
increasing `interpolation_factor` — the primary setting that controls the
density of the interpolated k-point mesh — until the change in mobility at all
doping and temperatures was less than $5\text{\,}\mathrm{\char 37\relax}$.

##  S6 High doping AMSET screening

Figure S7: amset screening results at high doping
(${10}^{19}\text{\,}{\mathrm{cm}}^{-3}$). Histograms of the temperature ($T$)
dependence of mobility for different scattering mechanisms at a low
($300\text{\,}\mathrm{K}500\text{\,}\mathrm{K}$) and b high
($800\text{\,}\mathrm{K}1000\text{\,}\mathrm{K}$) temperatures obtained from
23,000 materials. Calculations were performed using the amset package and are
separated into polar optical (PO), acoustic deformation potential (AD), and
ionized impurity (II) scattering. The “Overall” distribution gives the
temperature dependence using the total scattering rate ($\tau^{-1}$) obtained
from Matthiessen’s rule
($\tau^{-1}_{\mathrm{overall}}=\tau^{-1}_{\mathrm{po}}+\tau^{-1}_{\mathrm{ad}}+\tau^{-1}_{\mathrm{ii}}$).

##  S7 Explaining the difference between the single parabolic band model and
the AMSET distribution

In the main text, we note that the temperature dependence resulting from the
single parabolic band model acts as the upper bound to the AMSET distribution.
We believe this originates from the tendency of bands in realistic materials
to be less dispersive than a parabolic band ($E\propto k^{<2}$) rather than
more dispersive ($E\propto k^{>2}$). I.e., the dispersion typically becomes
linear away from the band edges. If $k^{<2}$, then group velocity increases
more slowly with energy than for a parabolic band. This makes higher-energy
carriers that get excited at higher temperatures less mobile, leading to a
more negative temperature dependence. In contrast if $k^{>2}$, then higher-
energy carriers would be relatively more mobile, leading to a less negative
temperature dependence.

##  S8 Explaining the variation in temperature dependence at a particular
polar phonon frequency

As noted in the main text, the temperature dependence of mobility when limited
by polar optical phonon scattering exhibits a strong correlation
($r^{2}=0.82$) to the polar phonon frequency, $\omega_{\mathrm{po}}$ (Fig.
S8a). However, despite this correlation, the distribution of $T$ dependencies
at any particular value of $\omega_{\mathrm{po}}$ is rather broad, indicating
that $\omega_{\mathrm{po}}$ alone can not explain a materials $T$ dependence.
This is illustrated in Fig. S8b, in which the $T$ dependencies of all
materials with
$$9.99\text{\,}\mathrm{THz}$<\omega_{\mathrm{po}}<$10.01\text{\,}\mathrm{THz}$$
are depicted as a histogram. While the distribution is centered at value of
$\sim-2.5$, the tails of the distribution cover a wide range from $-3.36$ to
$-1.95$. Other factors responsible for the range of $T$ dependencies can
include the shape of the bands (group velocities) and the behaviour of the
scattering rates across the band structure. We note that the other materials
parameters that enter the polar optical phonon scattering matrix element,
namely the static and high-frequency dielectric constants, are constant across
all temperatures and therefore play no role in controlling the temperature
dependence of mobility. The behaviour of scattering rates across the band
structure is computationally expensive to obtain and therefore cannot be used
as an efficient predictor of the $T$-dependence. Furthermore, as demonstrated
by the correlation plot in Fig. S8c, for materials with
$$9.99\text{\,}\mathrm{THz}$<\omega_{\mathrm{po}}<$10.01\text{\,}\mathrm{THz}$$
there is a moderate correlation ($r^{2}=0.49$) between the temperature
dependence of mobility limited by polar optical scattering and acoustic
deformation potential scattering. This correlation suggests that band
structure features are at least partially responsible for the distribution of
temperature dependencies at a specific $\omega_{\mathrm{po}}$. In this
section, we explore whether common band structure features can be used to
predict the $T$ dependence of mobility.

Figure S8: a Heatmap indicating the correlation between the polar optical (PO)
phonon frequency ($\omega_{\mathrm{po}}$) and the temperature dependence of
mobility (in the range $300\text{\,}\mathrm{K}500\text{\,}\mathrm{K}$) limited
by optical phonon scattering. The $r^{2}=0.808$ coefficient indicates strong
correlation. b Histogram of temperature dependence of mobility limited by
optical phonon scattering for systems where $9.99<\omega_{\mathrm{po}}<10.01$.
c The correlation between mobility limited by polar optical phonon scattering
and acoustic deformation potential (AD) scattering. The $r^{2}=0.493$
coefficient indicates moderate correlation. In a and c, the pink line is the
line of best fit, from which the $r^{2}$ value is calculated.

In Fig. S9 we depict the band structures of the materials with the largest and
smallest $T$-dependence of POP limited mobility with
$\omega_{\mathrm{po}}\sim$10\text{\,}\mathrm{THz}$$. These comprise $p$-Zr3N4
($\mu\propto T^{-3.36}$), $n$-GaI3Se(CH3)2 ($T^{-3.32}$), $n$-PtC4(NCl3)2
($T^{-3.13}$), $n$-Zn(CuO2)2 ($T^{-2.03}$), $n$-Mn2CrSbO6 ($T^{-1.97}$), and
$n$-Zr3N4 ($T^{-1.95}$). There are no obvious trends among the band
structures. For example, both large and small $T$ dependencies are seen in
band structures with dispersive band edges ($p$\- and $n$-type Zr3N4,
respectively). Similarly, most of the materials possess degenerate or nearly
degenerate band edges with multiple band extrema in the Brillouin zone.

Figure S9: Band structures of the materials with the largest and smallest
temperature dependence of mobility limited by polar optical phonon scattering,
where $9.99<\omega_{\mathrm{po}}<10.01$. a $p$-Zr3N4, b $n$-GaI3Se(CH3)2, c
$n$-PtC4(NCl3)2, d $n$-Zn(CuO2)2, e $n$-Mn2CrSbO6, and f $n$-Zr3N4. The shaded
regions indicate whether the transport is controlled by the conduction bands
($n$-type) or valence bands ($p$-type).

To further investigate the variation in the $T$ dependence of mobility at a
single $\omega_{\mathrm{po}}$, we calculate several band structure features
and look for correlations. These features comprise the conductivity effective
mass ($m^{*}_{c}$), the Seebeck effective mass ($m^{*}_{s}$), the Fermi
surface complexity factor44
($N_{v}K^{*}=\left[m^{*}_{s}/m^{*}_{c}\right]^{3/2}$), the smallest band width
of the bands within $0.3\text{\,}\mathrm{eV}$ of the band edge ($\Delta E$),
and the integrated density of states within $0.3\text{\,}\mathrm{eV}$ of the
band edge ($D_{\mathrm{int}}$). The conductivity effective mass was calculated
according to

$m_{c}^{*}(T,n)=\frac{\sigma(T,n)}{e^{2}\tau n},$ (S3)

where $T$ and $n$ are the temperature and carrier concentration set to
$300\text{\,}\mathrm{K}$ and ${10}^{16}\text{\,}{\mathrm{cm}}^{-3}$,
respectively, $\sigma$ is the conductivity calculated within the constant
relaxation time approximation with a lifetime
$\tau=${10}^{-14}\text{\,}\mathrm{s}$$, and $e$ is the electron charge. The
Seebeck effective mass was calculated using the methodology outlined in Ref.
44.

Figure S10: Correlations between band structure features and the temperature
of mobility limited by polar optical phonon scattering, for materials where
$9.99<\omega_{\mathrm{po}}<10.01$. a conductivity effective mass
($m^{*}_{c}$), b Seebeck effective mass ($m^{*}_{s}$), c Fermi surface
complexity factor ($N_{v}K^{*}=\left[m^{*}_{s}/m^{*}_{c}\right]^{3/2}$), d the
smallest band width of the bands within $0.3\text{\,}\mathrm{eV}$ of the band
edge ($\Delta E$), and e the integrated density of states within
$0.3\text{\,}\mathrm{eV}$ of the band edge ($D_{\mathrm{int}}$). The pink line
is the line of best fit from which the $r^{2}$ value is calculated.

In Fig. S10, we plot the correlation between band structure features and the
$T$ dependence of mobility limited by polar optical phonon scattering. For
each feature, we also calculate the $r^{2}$ score against the line of best fit
as a metric of the correlation. All features exhibit poor correlation with $T$
dependence. The integrated density of states (Fig. S10e) exhibits the greatest
correlation with an $r^{2}=0.269$ (i.e., the larger number of states is weakly
correlated with stronger T-dependence). The minimum band width within
$0.3\text{\,}\mathrm{eV}$ of the band edge is the second most correlated
feature with an $r^{2}=0.174$ (Fig. S10d). The remaining features
($m^{*}_{c}$, $m^{*}_{s}$, and Fermi surface complexity factor) all exhibit
$r^{2}$ values below 0.08, indicating essentially no correlation (Fig.
S10a–c).

Figure S11: Cross validation predictions of the machine learning model to
predict temperature dependence of mobility. Mobility was calculated using only
polar optical phonon scattering, for materials where
$9.99<\omega_{\mathrm{po}}<10.01$. The dashed line shows the ideal trend.

Clearly, none of the identified band structure features can account for the
variation in $T$ dependence of mobility at a particular
$\omega_{\mathrm{po}}$. To investigate this further, we trained a machine
learning model to assess whether a combination of band structure features can
predict the mobility behaviour. Machine learning was performed using a simple
regression model (random forest with 100 estimators). The dataset comprised
the 108 materials where
$$9.99\text{\,}\mathrm{THz}$<\omega_{\mathrm{po}}<$10.01\text{\,}\mathrm{THz}$$
and contained the five features outlined above ($m^{*}_{c}$, $m^{*}_{s}$,
$N_{v}K^{*}$, $D_{\mathrm{int}}$, $\Delta E$). Model performance was evaluated
using 5-fold cross validation. The trained models exhibited essentially no
improvement over the single features, with a cross validation $r^{2}$ score of
0.312 (Fig. S11), only marginally larger than for $D_{\mathrm{int}}$
($r^{2}=0.269$). Accordingly, we posit that the $T$ dependence must be
controlled by a complex interplay between band structure effects, wave
function overlaps, and the behaviour of the scattering rates across the band
structure.

## References

  * Rohatgi 2021 Rohatgi, A. Webplotdigitizer: Version 4.5. 2021; https://automeris.io/WebPlotDigitizer
  * Ma et al. 2020 Ma, J.; Meng, F.; He, J.; Jia, Y.; Li, W. Strain-Induced Ultrahigh Electron Mobility and Thermoelectric Figure of Merit in Monolayer $\alpha$-Te. _ACS Appl. Mater. Interfaces_ 2020, _12_ , 43901–43910, DOI: 10.1021/acsami.0c10236
  * Ge et al. 2020 Ge, Y.; Wan, W.; Ren, Y.; Liu, Y. Large Thermoelectric Power Factor of High-Mobility Transition-Metal Dichalcogenides with 1T′′ Phase. _Phys. Rev. Research_ 2020, _2_ , 013134, DOI: 10.1103/PhysRevResearch.2.013134
  * Meng et al. 2019 Meng, F.; Ma, J.; He, J.; Li, W. Phonon-Limited Carrier Mobility and Temperature-Dependent Scattering Mechanism of 3C-SiC from First Principles. _Phys. Rev. B_ 2019, _99_ , 045201, DOI: 10.1103/PhysRevB.99.045201
  * Su et al. 2020 Su, Y.; Cao, S.; Shi, L.-B.; Qian, P. Phonon-Limited Mobility for Novel Two-Dimensional Semiconductors of BC3 and C3N: First-Principles Calculation. _Appl. Surf. Sci._ 2020, _531_ , 147341, DOI: 10.1016/j.apsusc.2020.147341
  * Park et al. 2019 Park, J.; Xia, Y.; Ozoliņš, V. High Thermoelectric Power Factor and Efficiency from a Highly Dispersive Band in Ba2BiAu. _Phys. Rev. Applied_ 2019, _11_ , 014058, DOI: 10.1103/PhysRevApplied.11.014058
  * Poncé et al. 2019 Poncé, S.; Schlipf, M.; Giustino, F. Origin of Low Carrier Mobilities in Halide Perovskites. _ACS Energy Lett._ 2019, _4_ , 456–463, DOI: 10.1021/acsenergylett.8b02346
  * Ma et al. 2020 Ma, J.; Meng, F.; Xu, D.; Hu, R.; Luo, X. Electron Mobility and Mode Analysis of Scattering for $\beta$-Ga2O3 from First Principles. _J. Phys. Condens. Matter_ 2020, _32_ , 465704, DOI: 10.1088/1361-648X/aba8ca
  * Lee et al. 2020 Lee, N.-E.; Zhou, J.-J.; Chen, H.-Y.; Bernardi, M. Ab Initio Electron-Two-Phonon Scattering in GaAs from next-to-Leading Order Perturbation Theory. _Nat. Commun._ 2020, _11_ , 1607, DOI: 10.1038/s41467-020-15339-0
  * Bushick et al. 2020 Bushick, K.; Mengle, K. A.; Chae, S.; Kioupakis, E. Electron and Hole Mobility of Rutile GeO2 from First Principles: An Ultrawide-Bandgap Semiconductor for Power Electronics. _Appl. Phys. Lett._ 2020, _117_ , 182104, DOI: 10.1063/5.0033284
  * Li et al. 2019 Li, W.; Poncé, S.; Giustino, F. Dimensional Crossover in the Carrier Mobility of Two-Dimensional Semiconductors: The Case of InSe. _Nano Lett._ 2019, _19_ , 1774–1781, DOI: 10.1021/acs.nanolett.8b04799
  * Zhao et al. 2018 Zhao, Y.; Dai, Z.; Zhang, C.; Lian, C.; Zeng, S.; Li, G.; Meng, S.; Ni, J. Intrinsic Electronic Transport and Thermoelectric Power Factor in N-Type Doped Monolayer MoS2. _New J. Phys._ 2018, _20_ , 043009, DOI: 10.1088/1367-2630/aab338
  * Cao et al. 2018 Cao, J.; Querales-Flores, J. D.; Murphy, A. R.; Fahy, S.; Savić, I. Dominant Electron-Phonon Scattering Mechanisms in _n_ -Type PbTe from First Principles. _Phys. Rev. B_ 2018, _98_ , 205202, DOI: 10.1103/PhysRevB.98.205202
  * Zhao et al. 2020 Zhao, Y.; Lian, C.; Zeng, S.; Dai, Z.; Meng, S.; Ni, J. Anomalous Electronic and Thermoelectric Transport Properties in Cubic Rb3AuO Antiperovskite. _Phys. Rev. B_ 2020, _102_ , 094314, DOI: 10.1103/PhysRevB.102.094314
  * Poncé et al. 2018 Poncé, S.; Margine, E. R.; Giustino, F. Towards Predictive Many-Body Calculations of Phonon-Limited Carrier Mobilities in Semiconductors. _Phys. Rev. B_ 2018, _97_ , 121201, DOI: 10.1103/PhysRevB.97.121201
  * Ma et al. 2018 Ma, J.; Chen, Y.; Li, W. Intrinsic Phonon-Limited Charge Carrier Mobilities in Thermoelectric SnSe. _Phys. Rev. B_ 2018, _97_ , 205207, DOI: 10.1103/PhysRevB.97.205207
  * Park et al. 2020 Park, J.; Xia, Y.; Ganose, A. M.; Jain, A.; Ozoliņš, V. High Thermoelectric Performance and Defect Energetics of Multipocketed Full Heusler Compounds. _Phys. Rev. Applied_ 2020, _14_ , 024064, DOI: 10.1103/PhysRevApplied.14.024064
  * Zhou and Bernardi 2019 Zhou, J.-J.; Bernardi, M. Predicting Charge Transport in the Presence of Polarons: The beyond-Quasiparticle Regime in SrTiO3. _Phys. Rev. Research_ 2019, _1_ , 033138, DOI: 10.1103/PhysRevResearch.1.033138
  * Huang et al. 2020 Huang, Q.; Ma, J.; Xu, D.; Hu, R.; Luo, X. High Thermoelectric Figure of Merit in Monolayer Tl2O from First Principles. _J. Appl. Phys._ 2020, _128_ , 185111, DOI: 10.1063/5.0030588
  * Ding et al. 2021 Ding, J.; Liu, C.; Xi, L.; Xi, J.; Yang, J. Thermoelectric Transport Properties in Chalcogenides ZnX (X=S, Se): From the Role of Electron-Phonon Couplings. _J. Materiomics_ 2021, _7_ , 310–319, DOI: 10.1016/j.jmat.2020.10.007
  * Wu et al. 2020 Wu, Y.; Hou, B.; Chen, Y.; Cao, J.; Ma, C.; Shao, H.; Zhang, Y.; Lu, Z.; Zhu, H.; Fang, Z.; Zhang, R.; Zhang, H. Towards Accurate Predictions of Carrier Mobilities and Thermoelectric Performances in 2D Materials. 2020
  * Ge et al. 2020 Ge, Y.; Wan, W.; Ren, Y.; Li, F.; Liu, Y. Phonon-Limited Electronic Transport of Two-Dimensional Ultrawide Bandgap Material h-BeO. _Appl. Phys. Lett._ 2020, _117_ , 123101, DOI: 10.1063/5.0022426
  * Kane 1957 Kane, E. O. Band Structure of Indium Antimonide. _J. Phys. Chem. Solids_ 1957, _1_ , 249–261, DOI: 10.1016/0022-3697(57)90013-6
  * Tu 2011 Tu, L. W. _An Introduction to Manifolds_ , 2nd ed.; Universitext; Springer: New York, 2011
  * Ganose et al. 2021 Ganose, A. M.; Park, J.; Faghaninia, A.; Woods-Robinson, R.; Persson, K. A.; Jain, A. Efficient Calculation of Carrier Scattering Rates from First Principles. _Nat. Commun._ 2021, _12_ , 2222, DOI: 10.1038/s41467-021-22440-5
  * Park et al. 2021 Park, J.; Xia, Y.; Ozoliņš, V.; Jain, A. Optimal Band Structure for Thermoelectrics with Realistic Scattering and Bands. _npj Comp. Mater._ 2021, _7_ , 1–13, DOI: 10.1038/s41524-021-00512-w
  * Rode and Knight 1971 Rode, D. L.; Knight, S. Electron Transport in GaAs. _Phys. Rev. B_ 1971, _3_ , 2534–2541, DOI: 10.1103/PhysRevB.3.2534
  * Zhou and Bernardi 2016 Zhou, J.-J.; Bernardi, M. _Ab Initio_ Electron Mobility and Polar Phonon Scattering in GaAs. _Phys. Rev. B_ 2016, _94_ , 201201, DOI: 10.1103/PhysRevB.94.201201
  * Rode 1975 Rode, D. _Semiconductors and Semimetals_ ; Elsevier, 1975; Vol. 10; pp 1–89, DOI: 10.1016/S0080-8784(08)60331-2
  * Jain et al. 2013 Jain, A.; Ong, S. P.; Hautier, G.; Chen, W.; Richards, W. D.; Dacek, S.; Cholia, S.; Gunter, D.; Skinner, D.; Ceder, G.; Persson, K. A. Commentary: The Materials Project: A Materials Genome Approach to Accelerating Materials Innovation. _APL Mater._ 2013, _1_ , 011002, DOI: 10.1063/1.4812323
  * Perdew et al. 1996 Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. _Phys. Rev. Lett._ 1996, _77_ , 3865–3868, DOI: 10.1103/PhysRevLett.77.3865
  * Kresse and Furthmüller 1996 Kresse, G.; Furthmüller, J. Efficiency of Ab-Initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set. _Comput. Mater. Sci._ 1996, _6_ , 15–50, DOI: 10.1016/0927-0256(96)00008-0
  * Kresse and Furthmüller 1996 Kresse, G.; Furthmüller, J. Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. _Phys. Rev. B_ 1996, _54_ , 11169–11186, DOI: 10.1103/PhysRevB.54.11169
  * Jain et al. 2011 Jain, A.; Hautier, G.; Ong, S. P.; Moore, C. J.; Fischer, C. C.; Persson, K. A.; Ceder, G. Formation Enthalpies by Mixing GGA and GGA \+ U Calculations. _Phys. Rev. B_ 2011, _84_ , 045115, DOI: 10.1103/PhysRevB.84.045115
  * Blöchl 1994 Blöchl, P. E. Projector Augmented-Wave Method. _Phys. Rev. B_ 1994, _50_ , 17953–17979, DOI: 10.1103/PhysRevB.50.17953
  * Kresse and Joubert 1999 Kresse, G.; Joubert, D. From Ultrasoft Pseudopotentials to the Projector Augmented-Wave Method. _Phys. Rev. B_ 1999, _59_ , 1758–1775, DOI: 10.1103/PhysRevB.59.1758
  * Ong et al. 2013 Ong, S. P.; Richards, W. D.; Jain, A.; Hautier, G.; Kocher, M.; Cholia, S.; Gunter, D.; Chevrier, V. L.; Persson, K. A.; Ceder, G. Python Materials Genomics (Pymatgen): A Robust, Open-Source Python Library for Materials Analysis. _Comput. Mater. Sci._ 2013, _68_ , 314–319, DOI: 10.1016/j.commatsci.2012.10.028
  * de Jong et al. 2015 de Jong, M.; Chen, W.; Angsten, T.; Jain, A.; Notestine, R.; Gamst, A.; Sluiter, M.; Krishna Ande, C.; van der Zwaag, S.; Plata, J. J.; Toher, C.; Curtarolo, S.; Ceder, G.; Persson, K. A.; Asta, M. Charting the Complete Elastic Properties of Inorganic Crystalline Compounds. _Sci. Data_ 2015, _2_ , 150009, DOI: 10.1038/sdata.2015.9
  * Petousis et al. 2017 Petousis, I.; Mrdjenovich, D.; Ballouz, E.; Liu, M.; Winston, D.; Chen, W.; Graf, T.; Schladt, T. D.; Persson, K. A.; Prinz, F. B. High-Throughput Screening of Inorganic Compounds for the Discovery of Novel Dielectric and Optical Materials. _Sci. Data_ 2017, _4_ , 160134, DOI: 10.1038/sdata.2016.134
  * Jie et al. 2019 Jie, J.; Weng, M.; Li, S.; Chen, D.; Li, S.; Xiao, W.; Zheng, J.; Pan, F.; Wang, L. A New MaterialGo Database and Its Comparison with Other High-Throughput Electronic Structure Databases for Their Predicted Energy Band Gaps. _Sci. China Technol. Sci._ 2019, _62_ , 1423–1430, DOI: 10.1007/s11431-019-9514-5
  * Dunn et al. 2020 Dunn, A.; Wang, Q.; Ganose, A.; Dopp, D.; Jain, A. Benchmarking Materials Property Prediction Methods: The Matbench Test Set and Automatminer Reference Algorithm. _npj Comput Mater_ 2020, _6_ , 138, DOI: 10.1038/s41524-020-00406-3
  * Ward et al. 2018 Ward, L. et al. Matminer: An Open Source Toolkit for Materials Data Mining. _Comput. Mater. Sci._ 2018, _152_ , 60–69, DOI: 10.1016/j.commatsci.2018.05.018
  * Jain et al. 2015 Jain, A.; Ong, S. P.; Chen, W.; Medasani, B.; Qu, X.; Kocher, M.; Brafman, M.; Petretto, G.; Rignanese, G.-M.; Hautier, G., et al. FireWorks: a dynamic workflow system designed for high-throughput applications. _Concurr. Comput._ 2015, _27_ , 5037–5059, DOI: 10.1002/cpe.3505
  * Gibbs et al. 2017 Gibbs, Z. M.; Ricci, F.; Li, G.; Zhu, H.; Persson, K.; Ceder, G.; Hautier, G.; Jain, A.; Snyder, G. J. Effective Mass and Fermi Surface Complexity Factor from Ab Initio Band Structure Calculations. _Npj Comput. Mater._ 2017, _3_ , 1–6, DOI: 10.1038/s41524-017-0013-3


# Bayesian analysis of the spin distribution of LIGO/Virgo black holes

Juan García-Bellidoa<EMAIL_ADDRESS>José Francisco Nuño Silesa
<EMAIL_ADDRESS>Ester Ruiz Moralesb<EMAIL_ADDRESS>
ORCID: 0000-0002-0995-595X
aInstituto de Física Teórica UAM-CSIC, Universidad Autonóma de Madrid,
Cantoblanco, 28049 Madrid, Spain
bDepartamento de Física, ETSIDI, Universidad Politécnica de Madrid, 28012
Madrid, Spain

###### Abstract

Gravitational wave detection from binary black hole (BBH) inspirals has become
routine thanks to the LIGO/Virgo interferometers. The nature of these black
holes remains uncertain. We study here the spin distributions of LIGO/Virgo
black holes from the first catalogue GWTC-1 and the first four published BBH
events from run O3. We compute the Bayes evidence for several independent
priors: flat, isotropic, spin-aligned and anti-aligned. We find strong
evidence for low spins in all of the cases, and significant evidence for small
isotropic spins versus any other distribution. When considered as a
homogeneous population of black holes, these results give support to the idea
that LIGO/Virgo black holes are primordial.

gravitational waves, black hole binaries, spin

††preprint: IFT–UAM/CSIC–18–76, CERN–TH–2018–158

##  I Introduction

The regular detection with laser interferometers Abbott et al. (2019a); Abbott
et al. (2020a); Abbott et al. (2020b); Abbott et al. (2020c, d) of
gravitational wave (GW) events from BBH mergers has opened a new window into
the universe, and in particular to the exploration of the nature of black
holes.

Before the first BBH detections by LIGO, stellar black holes with masses in
the range $5-15~{}M_{\odot}$ had been detected as components of X-ray
binaries, only a few intermediate mass black holes (IMBH) were known with
masses above $100\,M_{\odot}$, while supermassive black holes (SMBH) were
known to exist at the centers of all galaxies. The origin of black holes in
such a large range of masses remains a mystery, and one fascinating
possibility is that part of these black holes are primordial in origin Clesse
and García-Bellido (2015); García-Bellido (2017). In fact, soon after the
first detection of black hole mergers by LIGO Abbott et al. (2016), there were
claims of their primordial nature Bird et al. (2016); Clesse and García-
Bellido (2017); Sasaki et al. (2016). Since then, the best scenario consistent
with all observational constraints so far Carr et al. (2021) is that of
spatially clustered and broad-mass distributed primordial black holes Clesse
and García-Bellido (2018); García-Bellido (2017).

As the number of binary black hole merger events detected with GW
interferometers increases, a new population of black holes is emerging with
unexpected properties in terms of their masses and spins. These properties
significantly differ from previous black holes detected through X-rays, via
stellar dynamics around SMBH at the center of galaxies, or via gravitational
lensing effects for IMBH. When the range of masses and distances accessible by
the GW interferometers improved, events like GW190425 and GW190521 appeared
with BH masses in the lower and upper mass gaps, challenging existing
astrophysical BH formation models. Moreover, events with small mass ratios
$q\ll 1$, like GW190814, are also difficult to generate in stellar binary
formation models, due to the expected mass transfer among the binary
components, see however Zevin et al. (2020).

One of the most striking features of this new population of black holes
detected in GW events is that they all seem to have a small spin. Although the
individual spin of each of the black holes in the binary is poorly
constrained,111It is worth pointing out, however, that the best determined
spin for a single black hole is that of the massive companion of GW190814,
with $S_{1}<0.07$ at 90% c.l., thanks to its particularly low mass ratio
$m_{2}/m_{1}=0.11$. a derived quantity called the effective spin, $\chi_{\rm
eff}$, can be well inferred from the GW waveform. The O1+O2 events Abbott et
al. (2019a) plus the four run O3 published events Abbott et al. (2020a);
Abbott et al. (2020b); Abbott et al. (2020c, d) show that, in almost all the
BBHs mergers events, the inferred $\chi_{\rm eff}$ posteriors cluster around
zero and are narrowly peaked. This single observation already constrains some
BBH formation models, because even if it is possible to generate a single
merger event with very small $\chi_{\rm eff}$ in almost all physical scenarios
of BBH formation, the fact that the whole population of BBH mergers have very
small $\chi_{\rm eff}$ cannot be explained by models involving high aligned
spins of the two BH’s in the binary.

Even before the first BBH merger detections, the effective spin had already
been considered an optimal variable to discriminate among formation channels.
By simulating populations of BBH mergers with different intrinsic spin values
and spin alignments, the authors in Vitale et al. (2017); Talbot and Thrane
(2017) showed that the effective spin could be relevant to distinguish the
astrophysical environment in which these binary systems formed, since binaries
evolved by dynamical interactions were expected to have their spins
isotropically oriented, while those coming from the evolution of an isolated
binary star system were more likely to have their spins highly aligned (see
however Belczynski et al. (2020) for updates on these models.).

The magnitude of the intrinsic spin of each black hole in the binary, though
much poorly constrained from the data than the effective spin, is important to
determine how these black holes were formed, and in particular to distinguish
a possible stellar versus primordial origin. Stellar black holes are expected
to have nonzero spin, due to conservation of angular momentum. On the
contrary, since the size of their Schwarzschild radius is essentially
identical to the size of the causal horizon at the moment of their formation,
primordial black holes are all expected to have zero or very near zero spin at
formation De Luca et al. (2019), although subsequent accretion could enhance
it slightly De Luca et al. (2020).

On the other hand, it has been argued that the more massive LIGO/Virgo BH
could originate from previous mergers Gerosa and Berti (2017); Fishbach et al.
(2017). However, it is known that the spin distribution of second generation
black holes is peaked very far from zero, near $S\sim 0.686$ Barausse and
Rezzolla (2009), and this is in disagreement with the fact that most of the
massive BH in LIGO/Virgo seem to have very low spin. In a careful analysis,
Ref Kimball et al. (2020) showed that GWTC-1 BBH catalog is consistent with
having no hierarchical mergers.

When the first detections of the O1 run became available, Farr et al. Farr et
al. (2017, 2018) analysed the discriminating power of the $\chi_{\rm eff}$
distribution with the first four GW LIGO/Virgo events, assuming equal masses
($q=m_{2}/m_{1}=1$) for all events and approximating the posterior $\chi_{\rm
eff}$ distributions by Gaussians. They compared the odds ratios for several
models of spin (modulus and orientations) of the underlying BH population,
finding a preference for either a population with an isotropic spin
distribution or with low intrinsic spin values of the merging black holes. A
more realistic analysis of the first six LIGO/Virgo mergers was made in Tiwari
et al. (2018), including Bayesian methods and taking into account the
$q-\chi_{\rm eff}$ correlations, reaching the general conclusion that highly
spinning black holes were disfavoured against low spins, see also Fernandez
and Profumo (2019).

With the publication of the GWTC-1 catalog, the LIGO/Virgo collaboration (LVC)
made a population analysis of the mass, redshift and spin properties of the
ten BBH mergers detected in O1+O2 runs Abbott et al. (2019b).222A hierarchical
Bayesian analysis of spin distributions was also performed in Miller et al.
(2020). Another Bayesian analysis, studying mass distributions and merger
rates but not spins, was done in Hall et al. (2020). Under different
assumptions for the parameters of the population models, they observed a
common trend of the inferred distribution for the BH spins to decrease with
increasing spin magnitude, but were not able to place strong constraints on
spin orientations, concluding that black hole spin measurements were not
informative enough at that moment to discern between isotropic and aligned
orientation distribution via $\chi_{\rm eff}$.

In this work, we improve such discriminating power by incorporating the
projection along the total angular momentum of the spin of the final BH formed
after the merger $a_{\rm f}$ as another variable of the BBH merger population
analysis. The effective spin and the final spin can be measured independently,
since $\chi_{\rm eff}$ is inferred from the inspiral part of the waveform
while $a_{\rm f}$ can be determined also from the ringdown part. We follow a
multivariate approach keeping all correlations among the variables
($q,\,\chi_{\rm eff},\,a_{\rm f}$) of each event, and show that this
significantly improves the Bayesian evidence for an underlying population of
black hole components with small spin magnitude and isotropic orientation.

There are still many uncertainties on the full spin (magnitude and
orientation) distribution of the different BBH formation channels, since the
computational models used to predict these properties depend on many
assumptions about poorly understood environmental conditions related to the
formation and evolution of the binaries. There are also many unknown aspects
related to the mass, spatial distribution and dynamics of BBH of primordial
origin. For this reason, instead of considering an specific astrophysical or
primordial model for the mass and spin distributions, we have chosen general
hypothesis on the underlying distributions of black hole spin magnitudes and
orientations and compare the different hypothesis in a full Bayesian analysis.

In Section 1, we describe the spin variables and parameters that will be
considered in our analysis of the BBH population. In Section 2, we describe
the method employed in calculating the Bayesian evidence from the published
parameter estimation samples, and define the population hypotheses together
with their priors for the parameters considered in our analysis. In Section 3
we use Bayesian methods to obtain the evidence for each hypothesis and the
corresponding Bayes ratios. Section 4 is devoted to the hierarchical modelling
method that allows us to infer a posterior distribution for the spin magnitude
of the underlying BH population. In Section 5 we give our conclusions.

##  II Spin observables from BBH events

In this section we will describe the observables relevant to the analysis of
spins in LVC black hole binaries. The main observables that we will consider
are the weighted-averaged effective projected spin, $\chi_{\rm eff}$, and the
the projection along the total angular momentum of the final spin $a_{\rm f}$.
We could have also chosen the effective precession spin parameter, $\chi_{p}$,
but this quantity is much worse measured in the LVC events published so far
and is mainly determined by the prior.

###  II.1 The effective spin

A derived physical quantity that can be measured very well by LIGO from the
waveform templates is the mass-averaged (so called effective) projected spin,

$\chi_{\rm
eff}=\frac{m_{1}\,\chi_{1}+m_{2}\,\chi_{2}}{m_{1}+m_{2}}=\frac{\chi_{1}+q\,\chi_{2}}{1+q}\,,\vspace{1mm}$
(1)

where $\chi_{i}\equiv(\vec{S}_{i}\cdot\vec{L})/(|\vec{L}|\,m_{i}^{2})$ is the
individual mass-weighted projection of each black hole spin onto the orbital
angular momentum $\vec{L}$ of the binary. Thus, $\chi_{\rm eff}$ gives some
information about the spin orientations of the inspiralling black holes w.r.t.
the orbital angular momentum. Here $q=m_{2}/m_{1}$ is the binary mass ratio,
in the range $q\in[0,1]$.

A derived quantity that will be useful is the so-called “chirp mass”,

${\cal
M}_{c}=\frac{(m_{1}\,m_{2})^{3/5}}{(m_{1}+m_{2})^{1/5}}=m_{1}\,\frac{q^{3/5}}{(1+q)^{1/5}}\,,$
(2)

which gives the mass ratio

$q=\frac{3\alpha\,(2/3)^{2/3}+(9\alpha+\sqrt{81\alpha^{2}-12\alpha^{3}})^{2/3}}{3\,(2/3)^{1/3}(9\alpha+\sqrt{81\alpha^{2}-12\alpha^{3}})^{1/3}}\,,$
(3)

where $\alpha=({\cal M}_{c}/m_{1})^{5}$.

Figure 1: The multidimensional LVC likelihoods projected on the planes
$(q,\,a_{\rm f})$ (top row), $(q,\,\chi_{\rm eff})$ (middle row) and
$(\chi_{\rm eff},\,a_{\rm f})$ (bottom row) for the ten GWTC-1 events of
LIGO/Virgo (left and middle column) plus the four runO3 events (right column).
The dotted lines correspond to the fixed points $\chi_{\rm eff}=0$ and $a_{\rm
f}=0.686$. The curves on the $(\chi_{\rm eff},\,a_{\rm f})$ plane correspond
to $q=1$ (black), $q=0.75$ (purple), $q=0.4$ (dark blue) and $q=0.11$ (light
blue), respectively.

###  II.2 The final spin

The second best-measured spin-related quantity in LIGO/Virgo binaries is the
projection along the total angular momentum of the spin of the final black
hole after merging, $a_{\rm f}$. We use the approximate expressions given in
Hofmann et al. (2016) for $a_{\rm f}$ in the particular case of spinning but
non-precessing black hole binaries. We first define some quantities and then
we assemble everything together.

First we need to define $a_{\rm tot}$ and $a_{\rm eff}$ in terms of the
projected spins, $a_{i}\equiv(\vec{S}_{i}\cdot\vec{L})/(|\vec{L}|\,m_{i})$,

$\displaystyle a_{\rm tot}$ $\displaystyle=$
$\displaystyle\frac{a_{1}+q^{2}\,a_{2}}{(1+q)^{2}}\,,\hskip 14.22636pta_{\rm
eff}=a_{\rm tot}+\xi\,\nu\,(a_{1}+a_{2})\,,$ (4) $\displaystyle{\rm with}$
$\displaystyle\hskip 8.53581pt\nu(q)=\frac{q}{(1+q)^{2}}\,,\hskip
28.45274pt\xi=0.474046\,.$

We then have to define the energy, angular momentum and size of the Innermost
Stable Circular Orbit (ISCO),

$\displaystyle E_{\rm ISCO}(a)=\sqrt{1-2/(3\,r_{\rm ISCO}(a))}\,,$ (5)
$\displaystyle L_{\rm ISCO}(a)=\frac{2}{3\sqrt{3}}\Big{(}1+2\sqrt{3\,r_{\rm
ISCO}(a)-2}\Big{)}\,,$ (6) $\displaystyle r_{\rm ISCO}(a)=3+Z_{2}-{\rm
sign}(a)\sqrt{(3-Z_{1})(3+Z_{1}+2Z_{2})}\,,$ $\displaystyle
Z_{1}(a)=1+(1-a^{2})^{1/3}\left((1+a)^{1/3}+(1-a)^{1/3}\right)\,,$ (7)
$\displaystyle Z_{2}(a)=\sqrt{3a^{2}+Z_{1}^{2}}\,.$ (8)

Then the final spin is given by (the coefficients $k_{ij}$ can be found in
Table 1)

$\displaystyle a_{\rm f}$ $\displaystyle=$ $\displaystyle a_{\rm
tot}+\nu\Big{(}L_{\rm ISCO}(a_{\rm eff})-2a_{\rm tot}(E_{\rm ISCO}(a_{\rm
eff})-1)\Big{)}$ $\displaystyle+$
$\displaystyle\nu^{2}\left(k_{00}+k_{01}\,a_{\rm eff}+k_{02}\,a_{\rm
eff}^{2}+k_{03}\,a_{\rm eff}^{3}+k_{04}\,a_{\rm eff}^{4}\right)$
$\displaystyle+$ $\displaystyle\nu^{3}\left(k_{10}+k_{11}\,a_{\rm
eff}+k_{12}\,a_{\rm eff}^{2}+k_{13}\,a_{\rm eff}^{3}+k_{14}\,a_{\rm
eff}^{4}\right)$ $\displaystyle+$
$\displaystyle\nu^{4}\left(k_{20}+k_{21}\,a_{\rm eff}+k_{22}\,a_{\rm
eff}^{2}+k_{23}\,a_{\rm eff}^{3}+k_{24}\,a_{\rm eff}^{4}\right)$
$\displaystyle+$ $\displaystyle\nu^{5}\left(k_{30}+k_{31}\,a_{\rm
eff}+k_{32}\,a_{\rm eff}^{2}+k_{33}\,a_{\rm eff}^{3}+k_{34}\,a_{\rm
eff}^{4}\right)\,.$

$-5.97723$ | $3.39221$ | $4.48865$ | $-5.77101$ | $-13.0459$
---|---|---|---|---
$35.1278$ | $-72.9336$ | $-86.0036$ | $93.7371$ | $200.975$
$-146.822$ | $387.184$ | $447.009$ | $-467.383$ | $-884.339$
$223.911$ | $-648.502$ | $-697.177$ | $753.738$ | $1166.89$
Table 1:  The coefficients $k_{ij}$ for $i\in[0,3]$ and $j\in[0,4]$.

Figure 2: The priors in the planes $(q,\,a_{\rm f})$ (top), $(q,\,\chi_{\rm
eff})$ (middle) and $(\chi_{\rm eff},\,a_{\rm f})$ (bottom) for the three
hypothesis $H_{i=2,3,4}$: isotropic spin (left), aligned spin (middle) and
anti-aligned spin (right), for spins centered at $\mu=0,\,0.5$ and 1 (blue,
green and red resp.).

##  III Bayesian Population Analysis

The Bayes theorem relates the likelihood ${\cal L}(d_{j}|\theta)$ of the data
$d_{j}$ for a given set of parameters $\theta$, with the posterior probability
of the parameters given the data, ${\cal P}(\theta|d_{j})$, via the prior
knowledge about the parameters of the population model $i$, $\Pi_{i}(\theta)$,

$\displaystyle{\cal P}_{i}(\theta|d_{j})$ $\displaystyle=$
$\displaystyle\frac{{\cal L}(d_{j}|\theta)\,\Pi_{i}(\theta)}{E(d_{j})}\,,$
$\displaystyle E_{ij}(d)$ $\displaystyle=$ $\displaystyle\int d\theta\,{\cal
L}(d_{j}|\theta)\,\Pi_{i}(\theta)\,,$ (9)

where $E_{ij}(d)$ is the Bayesian evidence for the data $d_{j}$ and the
population model $i$. Here we will consider different priors,
$\Pi_{i}(\theta)$, for the distribution of the parameters in each of the spin
population models characterized by the hypothesis $H_{i}$. We first compute
the multidimensional likelihood for each BBH event from the LIGO/Virgo
published samples, and then calculate the priors associated with each
population hypothesis to derive the evidence (9).

###  III.1 Likelihoods for LIGO/Virgo BBH events

We calculate the multidimensional likelihoods (i.e. including correlations)
from the posteriors and priors provided by the LIGO public documentation
page.333https://dcc.ligo.org/cgi-bin/DocDB/DocumentDatabase We use the
parameter estimation samples for each event, which are given in terms of the
fundamental parameters
$(m_{1},m_{2},s_{1},s_{2},\cos\theta_{1},\cos\theta_{2})$, and using the
expressions of $\chi_{\rm eff}$ in (1) and $a_{\rm f}$ in (II.2) for each
sample point, we construct the full multidimensional posterior and prior
distributions for the derived parameters $\theta=\\{q,\,\chi_{\rm
eff},\,a_{\rm f}\\}$, marginalizing over the other parameters. Finally we
generate the multidimensional likelihoods ${\cal L}_{\rm LVC}(d_{j}|\theta)$
dividing the posteriors ${\cal P}_{\rm LVC}(d_{j}|\theta)$ by their
corresponding priors $\Pi_{\rm LVC}(\theta)$. We have plotted the LVC
likelihoods for the three parameters $\theta=\\{q,\,\chi_{\rm eff},\,a_{\rm
f}\\}$ in Fig. 1.

When using these likelihoods to infer population properties in our Bayesian
analysis, we will not include possible selection effects in the ($\chi_{\rm
eff},\,a_{\rm f},\,q$) variables, like the observational bias towards positive
$\chi_{\rm eff}$ values described in Ng et al. (2018), or possible effects
from the fact that the waveform bank only considers aligned spins. We expect
these effects to be small given the spin variables used, but should eventually
be considered in a future analysis.

###  III.2 Priors on the BBH population

Once we have chosen to use the subset of parameters $\theta=\\{q,\,\chi_{\rm
eff},\,a_{\rm f}\\}$ for the Bayesian inference, we need to specify the
different hypotheses on the spin distribution of the merging BBH population.
We have chosen 5 different hypotheses $H_{i}$ on the spin magnitude and
orientation, which in turn depend on “hyperparameters”
$\Lambda=\\{\sigma,\,\mu\\}$ that determine the width and the mean of the
distribution of the spin magnitude.

These five basic Hypotheses are:

  * •

$H_{0}$: isotropic spin orientation, $\cos\theta_{LS_{i}}\\!\in\\![-1,1]$,
flat prior on spin magnitude in the range $S_{i}\in[0,1]$, flat prior on mass
ratio, $q\in[0,1]$. This is the ”null” hypothesis.

  * •

$H_{1}$: isotropic spin orientation, $\cos\theta_{LS_{i}}\\!\in\\![-1,1]$,
Gaussian prior on spin magnitude with $\mu=0$ and $\sigma\in[0,1]$, flat prior
on mass ratio, $q\in[0,1]$.

  * •

$H_{2}$: isotropic spin orientation, $\cos\theta_{LS_{i}}\\!\in\\![-1,1]$,
Gaussian prior on spin magnitude with $\sigma=0.2$ and $\mu\in[0,1]$, flat
prior on mass ratio, $q\in[0,1]$.

  * •

$H_{3}$: aligned spin orientation, $\cos\theta_{LS_{i}}$ sampled from Gaussian
centered at $+1$ and width $0.05$, Gaussian prior on spin magnitude for
$\sigma=0.2$ and $\mu\in[0,1]$, flat prior on mass ratio, $q\in[0,1]$.

  * •

$H_{4}$: anti-aligned spin orientation, $\cos\theta_{LS_{i}}$ sampled from
Gaussian centered at $\pm 1$ and width $0.05$, Gaussian prior on spin
magnitude for $\sigma=0.2$ and $\mu\in[0,1]$, flat prior on mass ratio,
$q\in[0,1]$.

We generate the multidimensional prior probability distributions for each
hypothesis $\Pi_{i}(\theta)$ from $10^{5}$ random realizations in
($q,\cos\theta_{i},\mu,\sigma$), giving rise to the corresponding prior
distributions in the $\theta=\\{q,\,\chi_{\rm eff},\,a_{\rm f}\\}$ parameter
space. To compute these distributions we have used the semianalytic expression
of $a_{f}$ given in equation (II.2) in terms of the fundamental parameters
$(m_{1},m_{2},s_{1},s_{2},\cos\theta_{1},\cos\theta_{2})$.

Fig. 2 shows the projections on the planes $(q,\,a_{\rm f})$, $(q,\,\chi_{\rm
eff})$ and $(\chi_{\rm eff},\,a_{\rm f})$ of the prior probability
distributions $\Pi_{i}(q,\,\chi_{\rm eff},\,a_{\rm f})$ for the last three
hypotheses $H_{i=2}$ (isotropic spins, left column), $H_{i=3}$ (aligned spins,
central column) and $H_{i=4}$ (anti-aligned spins, right column). In each of
these plots, we have chosen for illustrative purposes three different choices
for the spin magnitude: Gaussian distributions centered at $\mu=0,\,0.5,\,1$
(blue, green and red respectively) with a common width $\sigma=0.2$. The black
curves in the $(q,\,a_{\rm f})$ and $(q,\,\chi_{\rm eff})$ planes corresponds
to the limit case of zero spin of the underlying black hole population, while
the black curve in the $(\chi_{\rm eff},\,a_{\rm f})$ plane corresponds to the
limiting case of $q=1$.

These plots show that the different hypotheses on the underlying spin
configurations populate very different areas of the $(q,\,\chi_{\rm
eff},\,a_{\rm f})$ parameter space, and therefore, taken as priors and
integrated over the whole parameter space with the likelihood of each event
shown in Fig.1, we expect very different Bayesian evidence and therefore
significantly informative Bayes ratios for the different hypotheses.

Comparing the distributions on the different planes we see that while the
$\chi_{\rm eff}$ distributions are flat with respect to $q$, a well known
degeneracy for random spins, the $a_{\rm f}$ distributions have a strong
dependence on $q$. Models with high and aligned or anti-aligned spins simply
cannot produce events with low $a_{\rm f}$ and $q$ values, as found for event
GW190814. Only populations with close to zero spin or high but isotropic spins
can populate the lower left corner of the $(q,a_{\rm f})$ plane. If the
forthcoming O3a catalog includes more events of low $q$, we expect a much
better determination of the spin of the underlying black hole population.

These five hypotheses give rise to multidimensional priors in the space of
derived parameters $\theta=\\{q,\,\chi_{\rm eff},\,a_{\rm f}\\}$, which can
then be used to integrate them, together with the LVC likelihoods, to obtain
the Bayesian evidence for each prior hypothesis $H_{i}$, which depend on
$\Lambda=\\{\sigma,\,\mu\\},$

$E_{ij}(\Lambda)=\\!\\!\int\\!d^{3}\theta\,{\cal
L}(h_{j}|\theta)\,\Pi_{i}(\theta|\Lambda)\,.$ (10)

Once we have obtained the LVC likelihoods from their parameter estimation
samples, and computed the individual priors for our hypothesis, we can perform
the 3D integration with Mathematica.

Figure 3: The log of the Bayes ratios in the variables $(q,\,\chi_{\rm
eff},\,a_{\rm f})$ for the four hypothesis: $H_{1}$ = isotropic spin centered
at zero and width $\sigma\in[0,\,1]$ (top left) and $H_{2,3,4}$ = isotropic
spin (top right), aligned spin (bottom left) and anti-aligned spin (bottom
right), for spin distributions centered at $\mu\in[0,\,1]$, for $\sigma=0.2$.
We also study the sensitivity of the Bayes ratios to the evidence in the 2D
planes $(q,\,\chi_{\rm eff})$, $(q,\,a_{\rm f})$ and $(\chi_{\rm eff},\,a_{\rm
f})$, with increasing values of the Bayes ratios. We find strong evidence
$|\ln B_{12}|\geq 6$ at small values of the spin, for all four Hypotheses
$H_{i}$ versus the (null) “All Flat” Hypothesis $H_{0}$, attaining the maximum
value of $|\ln B_{12}|\geq 8$ for $H_{4}$ at $\mu\leq 0.01$. The events
contributing to these Bayes factors are the GWTC-1 (10-event catalog) plus the
four published run-O3 events from LIGO/Virgo Collaboration.

##  IV Bayes Factors

In order to evaluate the global significance of a prior hypothesis from the
full BBH catalog, we compute the global Bayes factor for the whole BH
population. For this, we will multiply the individual Bayes ratios, assuming
that all events in the catalog are independent,

$\ln\,B^{i}_{12}(\Lambda)=\sum_{j=1}^{N}\ln\frac{E_{ij}}{E_{0}}(\Lambda)\,,$
(11)

for each hypothesis $H_{i}(\Lambda=\\{\sigma,\,\mu\\})$.

In practice what we do is a reweighting of the priors, using our population
model hypothesis, $\Pi_{i}(\theta)$, versus the published LIGO/Virgo priors,
$\Pi_{\rm LVC}(\theta)$, which depend on the assumptions of the experiment on
each event,

$\displaystyle E_{ij}(\Lambda)$ $\displaystyle=$ $\displaystyle\int
d^{3}\theta\,{\cal L}_{\rm LVC}(d_{j}|\theta)\,\Pi_{i}(\theta|\Lambda)$ (12)
$\displaystyle=$ $\displaystyle E_{\rm LVC}\,\int d^{3}\theta\,{\cal P}_{\rm
LVC}(\theta|d_{j})\,\frac{\Pi_{i}(\theta|\Lambda)}{\Pi_{\rm LVC}(\theta)}\,.$

Since $E_{\rm LVC}$ is the same for all hypotheses, it factors out in the
Bayes ratio (11). We can then perform the integral (12) in the full 3D
parameter space $\theta=\\{q,\,\chi_{\rm eff},\,a_{\rm f}\\}$ maintaining all
the correlations in the multivariate priors and posteriors.

We show in Fig. 3 the Bayes ratios of the four different hypotheses
$H_{i=1,2,3,4}$ with respect to the null hypothesis $H_{0}$ of flat
distributions for mass ratios, spin magnitudes and random spin orientations.
According to Jeffrey’s scale Jeffreys (1939), when $\ln\,B_{12}\geq 5$,
hypothesis 1 is significantly more likely that hypothesis 2, see however
Nesseris and García-Bellido (2013).

It is clear from figure 3 that low spin magnitudes are significantly preferred
for the BBH population of LIGO/Virgo events, with Bayes ratios above $\ln
B_{12}=5$ for $\sigma,\mu<0.2$, therefore we find strong evidence for small
values of the spin, for all four Hypotheses $H_{i}$ versus the ”null” (All
Flat) Hypothesis $H_{0}$. In the case of Hypothesis $H_{1}$, whith zero spin
and allowing for variable width, we find strong evidence ($\ln B_{12}=6.2$)
for relatively narrow spin distributions $\sigma=0.15$. This low and narrow
spin hypothesis represents the isotropic spin distribution that one would
expect from an underlying population of primordial black holes. On the other
hand, when we vary the spin magnitude as a Gaussian centered at $\mu$ for
fixed width $\sigma=0.2$, we find, in all three orientation Hypothesis
$H_{2,3,4}$, that the maximum evidence occurs again for $\mu=0$. Therefore we
conclude that, whatever the orientation, there is very strong evidence for low
spins in LIGO/Virgo BH.

Note that the aligned spin hypothesis $H_{3}$ has Bayes factors $\ln
B_{12}<-5$ for spins $\mu>0.25$, reaching large Bayes ratios $\ln B_{12}<-100$
for large spins $\mu>0.8$, therefore our analysis of LVC catalog strongly
disfavoures a population of BBH with aligned spins and magnitudes greater than
$0.25$. Note also that hypothesis $H_{1}$ and $H_{4}$, for isotropic and anti-
aligned spins, have a similar behaviour in their Bayes factors at low spins,
$\mu<0.5$, reaching values above $\ln B_{12}\simeq 6$ for $\mu<0.06$ in the
isotropic case, and above $\ln B_{12}\simeq 8$ for $\mu<0.02$ in the anti-
aligned case, giving slightly higher significance for anti-aligned versus
isotropic spins in the range $\mu<0.5$. However, for large spins, $\mu>0.5$,
the anti-aligned hypothesis is much more strongly disfavoured than the
isotropic one.

We also compute the Bayes factors for the three pairs ($q,\,\chi_{\rm eff}$),
($q,\,a_{\rm f}$) and ($\chi_{\rm eff},\,a_{\rm f}$) finding less evidence for
low spins compared with the full analysis. The least significant is the Bayes
ratio in the plane ($q,\,\chi_{\rm eff}$), giving only mild preference for low
spins. As we include information coming from $a_{\rm f}$, as in ($\chi_{\rm
eff},\,a_{\rm f}$), ($q,\,a_{\rm f}$) and ($q,\,\chi_{\rm eff},\,a_{\rm f}$),
the evidence for low spin rises, in some cases up to $\ln B_{12}\sim 6$ or
above, which is very strong evidence in favor of that hypothesis.

Historically, the emphasis has been focused on the importance of
distinguishing between aligned and anti-aligned astrophysical models of BBH
formation. Now that the evidence for small spins is so strong, the orientation
of the spins becomes less relevant. It is much more interesting to
characterize the posterior distributions of the whole population and to
quantify the deviation from zero spins according to the various hypotheses
$H_{i}$, which take into account orientation. This will be useful in order to
characterize the possible origin of LVC BBH events from primordial black hole
populations.

Figure 4: The posterior PDFs $P(a)$ for the spin magnitude in the case of
isotropic (left), aligned (center) and anti-aligned (right) hypothesis. We
have computed the PDFs for three different spin widths, $\sigma=0.05$ (blue),
$\sigma=0.1$ (red) and $\sigma=0.25$ (green). The PDFs do not seem to depend
much on the width $\sigma$ of the distribution. In all cases, the PDFs peak
close to zero spin, specially for $\sigma=0.25$. The vertical gray lines
correspond to the 90% c.l. limit, which is almost the same for all PDFs.

##  V Bayesian hierarchical modeling

Given that the random isotropic spin hypothesis, $H_{1}$, gives a value of
$\ln\,B_{12}\sim 6$, strongly favouring low spin magnitudes $a<0.2$ for the
combined LVC BBH events, we would like to explore in detail the black hole
population distribution of spin magnitudes.

In order to go forward in the information content with respect to the Bayesian
approach of the previous section, we will use Bayesian hierarchical modeling,
see e.g. Talbot and Thrane (2017), as a way to estimate the posterior
distribution for spin magnitudes of the LVC black hole population. Here the
spin magnitude $a$ enters indirectly through the dependence of other measured
quantities, like the effective spin $\chi_{\rm eff}$ or the final spin $a_{\rm
f}$, on it. Then the posterior PDF for $a$ is computed as

$P_{i}(a)=\Pi(a)\,\prod_{j=1}^{N}E_{ij}(a)\,,$ (13)

where $E_{ij}(a)$ are the Bayesian evidences (10) computed from the LVC
likelihoods with the spin-dependent priors $\Pi_{i}(q,\,\chi_{\rm
eff},\,a_{\rm f}|\Lambda)$, with $\Lambda=\\{\sigma=0.05,0.1,0.25,\;\mu=a\\}$,
and $\Pi(a)$ is assumed here to be flat in spin magnitude $a$, not to give any
prior preference for any spin magnitude. We have plotted in Fig. 4 the
posterior distributions for the three spin Hypothesis (isotropic, aligned and
anti-aligned). In all cases, the preferred spin magnitude is below $a<0.12$,
within 50% c.l. For example, for isotropic spins, currently favoured by the
data, we find $a<0.25$ at 90% c.l. For aligned spins the PDF is even more
strongly peaked around zero spin, with $a<0.1$ at 90% c.l. All this suggests
that the inclusion of the four published run O3 events, with increased
sensitivity of the detectors, has revealed a property of the population of
LIGO/Virgo black holes that was not present in previous studies Abbott et al.
(2019b).

##  VI Conclusions

The routine detection of BBH inspirals by LVC has opened the door to a
detailed exploration of the nature of black holes and their populations. The
online availability of the strain time-streams and the parameter estimation
samples for each event allows for an independent analysis, opening the
possibility to compute multivariate likelihoods for combinations of derived
parameters like the mass ratio, or the effective spin.

It is then a matter of personal choice which parameters to use in order to
infer properties of the populations of black holes detected by LVC. In this
paper we have concentrated on just three parameters ($q,\,\chi_{\rm
eff},\,a_{\rm f}$) that we believe capture the essence of the spin nature of
the population of LVC black holes.

We have put forward four different spin-magnitudes and spin-orientation prior
hypothesis, $H_{i=1,...,4}$ (i.e. isotropic, aligned and anti-aligned), to
compare with the (null) all-flat prior hypothesis, $H_{0}$, and conducted a
Bayesian analysis study to determine the goodness of a given spin-distribution
hypothesis for the whole population of LIGO/Virgo black holes.

We find that all spin-orientation hypothesis have a larger Bayes factor for
low spins. In some cases the log of the Bayes factor reaches values well above
five, thus signalling a strong evidence in favour of low spins (below
magnitude $a=0.2$). Moreover, the largest Bayes factors are obtained by spin
distributions peaked at zero spin (Bayes factors above 5 for widths
$\mu<0.2$), with very small width.

We also note that aligned spins are strongly disfavoured, specially for large
spin magnitudes (with the log of the Bayes factors as low as $-100$), as would
be expected from astrophysical black holes from isolated binaries. On the
other hand, we find that LIGO/Virgo black hole population has a preference
towards low spins with isotropic orientations, consistent with what one would
expect from primordial black holes in clusters Clesse and García-Bellido
(2017).

We have then computed the posterior PDF for the spin magnitude in the case of
the three alternative spin hypothesis, for different spin widths
($\sigma=0.05,\,0.10,\,0.25$). In all cases the spin distribution is peaked at
very low values ($a\simeq 0.1,\,0.05,\,0.1$ for isotropic, aligned and anti-
aligned cases, respectively), which clearly indicates a preference of the
whole LIGO/Virgo BH population for low spins, irrespective of orientation.

Therefore, we conclude that only using the BH spin as a discriminator between
the astrophysical versus primordial nature of LIGO/Virgo black holes, when
considered as a homogeneous population of black holes, our analysis seems to
suggest very strongly that they are consistent with being primordial.

## Acknowledgements

JGB and ERM thank the Theory Division for their generous hospitality during
their Sabbatical at CERN in 2018, when this work originated. The authors
acknowledge support from the Research Project PGC2018-094773-B-C32 (MINECO-
FEDER) and the Centro de Excelencia Severo Ochoa Program SEV-2016-0597.

## References

  * Abbott et al. (2019a) B. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. X 9, 031040 (2019a), eprint 1811.12907.
  * Abbott et al. (2020a) B. P. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J. Lett. 892, L3 (2020a), eprint 2001.01761.
  * Abbott et al. (2020b) R. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. D102, 043015 (2020b), eprint 2004.08342.
  * Abbott et al. (2020c) R. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J. Lett. 896, L44 (2020c), eprint 2006.12611.
  * Abbott et al. (2020d) R. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 125, 101102 (2020d), eprint 2009.01075.
  * Clesse and García-Bellido (2015) S. Clesse and J. García-Bellido, Phys. Rev. D92, 023524 (2015), eprint 1501.07565.
  * García-Bellido (2017) J. García-Bellido, J. Phys. Conf. Ser. 840, 012032 (2017), eprint 1702.08275.
  * Abbott et al. (2016) B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016), eprint 1602.03837.
  * Bird et al. (2016) S. Bird, I. Cholis, J. B. Muñoz, Y. Ali-Haimoud, M. Kamionkowski, E. D. Kovetz, A. Raccanelli, and A. G. Riess, Phys. Rev. Lett. 116, 201301 (2016), eprint 1603.00464.
  * Clesse and García-Bellido (2017) S. Clesse and J. García-Bellido, Phys. Dark Univ. 15, 142 (2017), eprint 1603.05234.
  * Sasaki et al. (2016) M. Sasaki, T. Suyama, T. Tanaka, and S. Yokoyama, Phys. Rev. Lett. 117, 061101 (2016), eprint 1603.08338.
  * Carr et al. (2021) B. Carr, S. Clesse, J. García-Bellido, and F. Kuhnel, Phys. Dark Univ. 31, 100755 (2021), eprint 1906.08217.
  * Clesse and García-Bellido (2018) S. Clesse and J. García-Bellido, Phys. Dark Univ. 22, 137 (2018), eprint 1711.10458.
  * Zevin et al. (2020) M. Zevin, M. Spera, C. P. L. Berry, and V. Kalogera, Astrophys. J. Lett. 899, L1 (2020), eprint 2006.14573.
  * Vitale et al. (2017) S. Vitale, R. Lynch, R. Sturani, and P. Graff, Class. Quant. Grav. 34, 03LT01 (2017), eprint 1503.04307.
  * Talbot and Thrane (2017) C. Talbot and E. Thrane, Phys. Rev. D 96, 023012 (2017), eprint 1704.08370.
  * Belczynski et al. (2020) K. Belczynski et al., Astron. Astrophys. 636, A104 (2020), eprint 1706.07053.
  * De Luca et al. (2019) V. De Luca, V. Desjacques, G. Franciolini, A. Malhotra, and A. Riotto, JCAP 1905, 018 (2019), eprint 1903.01179.
  * De Luca et al. (2020) V. De Luca, G. Franciolini, P. Pani, and A. Riotto, JCAP 2004, 052 (2020), eprint 2003.02778.
  * Gerosa and Berti (2017) D. Gerosa and E. Berti, Phys. Rev. D95, 124046 (2017), eprint 1703.06223.
  * Fishbach et al. (2017) M. Fishbach, D. E. Holz, and B. Farr, Astrophys. J. Lett. 840, L24 (2017), eprint 1703.06869.
  * Barausse and Rezzolla (2009) E. Barausse and L. Rezzolla, Astrophys. J. Lett. 704, L40 (2009), eprint 0904.2577.
  * Kimball et al. (2020) C. Kimball, C. Talbot, C. P. L. Berry, M. Carney, M. Zevin, E. Thrane, and V. Kalogera, Astrophys. J. 900, 177 (2020), eprint 2005.00023.
  * Farr et al. (2017) W. M. Farr, S. Stevenson, M. Coleman Miller, I. Mandel, B. Farr, and A. Vecchio, Nature 548, 426 (2017), eprint 1706.01385.
  * Farr et al. (2018) B. Farr, D. E. Holz, and W. M. Farr, Astrophys. J. Lett. 854, L9 (2018), eprint 1709.07896.
  * Tiwari et al. (2018) V. Tiwari, S. Fairhurst, and M. Hannam, Astrophys. J. 868, 140 (2018), eprint 1809.01401.
  * Fernandez and Profumo (2019) N. Fernandez and S. Profumo, JCAP 08, 022 (2019), eprint 1905.13019.
  * Abbott et al. (2019b) B. Abbott et al. (LIGO Scientific, Virgo), Astrophys. J. Lett. 882, L24 (2019b), eprint 1811.12940.
  * Miller et al. (2020) S. Miller, T. A. Callister, and W. Farr, Astrophys. J. 895, 128 (2020), eprint 2001.06051.
  * Hall et al. (2020) A. Hall, A. D. Gow, and C. T. Byrnes (2020), eprint 2008.13704.
  * Hofmann et al. (2016) F. Hofmann, E. Barausse, and L. Rezzolla, Astrophys. J. Lett. 825, L19 (2016), eprint 1605.01938.
  * Ng et al. (2018) K. K. Y. Ng, S. Vitale, A. Zimmerman, K. Chatziioannou, D. Gerosa, and C.-J. Haster, Phys. Rev. D98, 083007 (2018), eprint 1805.03046.
  * Jeffreys (1939) H. Jeffreys, _The Theory of Probability_ , Oxford U.P. (1939).
  * Nesseris and García-Bellido (2013) S. Nesseris and J. García-Bellido, JCAP 1308, 036 (2013), eprint 1210.7652.