Searching for ultralight bosons within spin measurements of a population of binary black hole mergers

TL;DR

This paper proposes a hierarchical Bayesian method to combine multiple gravitational-wave black hole spin measurements, enabling stronger constraints on ultralight boson masses within a specific range, despite individual measurement uncertainties.

Contribution

It introduces a novel hierarchical Bayesian inference approach to improve constraints on ultralight bosons using multiple black hole spin measurements from gravitational-wave data.

Findings

Hundreds of high SNR detections can exclude certain boson mass ranges.

Current LIGO/Virgo data are insufficient for conclusive constraints.

The method's effectiveness depends on black hole spin distribution and boson mass.

Abstract

Ultralight bosons can form clouds around rotating black holes if their Compton wavelength is comparable to the black hole size. The boson cloud spins down the black hole through a process called superradiance, lowering the black hole spin to a characteristic value. It has been suggested that spin measurements of the black holes detected by ground-based gravitational-wave detectors can be used to constrain the mass of ultralight bosons. Unfortunately, a measurement of the individual black hole spins is often uncertain, resulting in inconclusive results. Instead, we use hierarchical Bayesian inference to combine information from multiple gravitational-wave sources and obtain stronger constraints. We show that hundreds of high signal-to-noise ratio gravitational-wave detections are enough to exclude (confirm) the existence of non-interacting bosons in the mass range $\left[10^{-13},3\times 10^{-12}\right]~\rm{eV}$ $\left([10^{-13},10^{-12}]~\rm{eV}\right)$. The precise number depends on the distribution of black hole spins at formation and the mass of the boson. From the few uninformative spin measurements of binary black hole mergers detected by LIGO and Virgo in their first two observing runs, we cannot draw statistically significant conclusions.

Tables

Table 1. Table 1: The estimated numbers of high SNR detections to rule out or confirm bosons for different combinations of formation spin distribution and boson.

Population models	Flat spin	Low spin
Astrophysical 777The statistical requirement is ${ℬ}_A^B=0.01$ to rule out bosons within $[{10}^{-13},3\times {10}^{-12}]$ eV.	${30}_{-25}^{+135}$	${140}_{-105}^{+120}$
${\mu }_s={10}^{-13}$ eV 888The statistical requirement is ${ℬ}_A^B=100$ to confirm the bosons.	${25}_{-15}^{+95}$	${80}_{-70}^{+210}$
${\mu }_s={10}^{-12}$ eV. 8	${65}_{-55}^{+165}$	${155}_{-145}^{+345}$ 999The upper bound is only an approximation since even using all the simulated signals we do not reach the desired threshold ${ℬ}_A^B=100$.

Equations

$$\begin{array}{rl}{\Gamma }_{[nlm]}^{\mathrm{inst}}=& {\mu }_s({\mu }_{s}M{)}^{4l+4}(m\chi -2{\mu }_s{r}_{+})\\ & \times \frac{2^{4l+2}(2l+1+n)!}{(l+1+n{)}^{2l+4}n!}{\left[\frac{l!}{(2l)!(2l+1)!}\right]}^2\\ & \times \prod _{k=1}^l[k^2(1-{\chi }^2)+(m\chi -2{\mu }_s{r}_{+}{)}^2],\end{array}$$

$$\frac{\alpha }{m}<\frac{\chi }{2\left(1+\sqrt{1-{\chi }^2}\right)}<\frac{1}{2}.$$

$${\chi }_{[nlm]}^{\infty }=\frac{4\alpha m}{4{\alpha }^2+m^2}.$$

$${\mathcal{H}}_{\mathrm{B}}:{\chi }_M=\{\begin{array}{ll}{\chi }_{[nlm]},& \text{if }{\tau }_{[nlm]}^{\mathrm{grow}}<{\tau }_s<{\tau }_{[n(l+1)(m+1)]}^{\mathrm{grow}}\\ {\chi }_F,& \text{otherwise}\end{array}$$

$${\mathcal{H}}_{\mathrm{A}}:{\chi }_M={\chi }_F.$$

$$p(\Lambda |\mathbf{d})\propto \pi (\Lambda )\prod _k^N\left[\frac{1}{\Sigma (\Lambda )}\int p(\theta |\Lambda )p(d_k|\theta )d\theta \right]$$

$$\Sigma (\Lambda )=\int p(\theta |\Lambda )p_{\det }(\theta )d\theta ,$$

$$p\left({\Lambda }_{{\mathcal{H}}_{\mathrm{B}}}|\mathbf{d},{\mathcal{H}}_{\mathrm{B}}\right)\propto \pi ({\Lambda }_{{\mathcal{H}}_{\mathrm{B}}})\prod _k^N\left\{\frac{1}{\Sigma ({\Lambda }_{{\mathcal{H}}_{\mathrm{B}}})}\int p(d_k|{\theta }_{{\mathcal{H}}_{\mathrm{B}}})\pi (M_1,M_2,{\tau }_s)\prod _{i=1}^2\left[p({\chi }_{M,i}|{\Lambda }_{{\mathcal{H}}_{\mathrm{B}}},M_i,{\tau }_s)d{M}_{i}d{\chi }_{M,i}\right]d{\tau }_s\right\},$$

$$\begin{array}{rl}& p({\chi }_{M,i}|\alpha ,\beta ,{\mu }_s,M_i,{\tau }_s)\\ & ={\int }_0^{1}p({\chi }_{M,i}|{\mu }_s,M_i,{\tau }_s,{\chi }_{F,i})p({\chi }_{F,i}|\alpha ,\beta )d{\chi }_{F,i},\end{array}$$

$$\begin{array}{rl}& p({\chi }_{M,i}|{\mu }_s,M_i,{\tau }_s,{\chi }_{F,i})\\ & =\delta \left[{\chi }_{M,i}-{\chi }_{[nlm]}({\mu }_s,M_i,{\tau }_s,{\chi }_{F,i})\right]\Theta \left[{\chi }_{F,i}-{\chi }_{[nlm]}\right]\\ & +\delta ({\chi }_{M,i}-{\chi }_{F,i})\Theta \left[{\chi }_{[nlm]}-{\chi }_{F,i}\right],\end{array}$$

$$\begin{array}{rl}p({\chi }_{M,i}|\alpha ,\beta ,& {\mu }_s,M_i,{\tau }_s)\approx f_{\mathrm{SR}}\delta ({\chi }_{M,i}-{\chi }_{[nlm]})\\ & +p({\chi }_{M,i}|\alpha ,\beta )\Theta ({\chi }_{[nlm]}-{\chi }_{M,i}),\end{array}$$

$$\begin{array}{rl}& f_{\mathrm{SR}}({\mu }_s,M_i,{\tau }_s,\alpha ,\beta )\\ \equiv & {\int }_0^{1}p({\chi }_{F,i}|\alpha ,\beta )\Theta \left[{\chi }_{F,i}-{\chi }_{[nlm]}({\mu }_s,M_i,{\tau }_s)\right]d{\chi }_{F,i}.\end{array}$$

Full text

I Introduction

Ultralight bosons with masses $\lesssim 10^{-11}\,\mathrm{eV}$, including axionlike particles [1, 2, 3], dilatons and, moduli [4, 5, 6] and fuzzy dark matter [7, 8, 9, 10, 11], have been proposed as a potential solution to various problems ranging from fundamental physics to cosmology [12, 2, 3, 13, 14, 15, 6, 16, 17, 18, 19, 20, 21]. Efforts are underway to search for these ultralight bosons using table-top experiments or astronomical observations [22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53]. When the Compton wavelength of the hypothetical boson is comparable to the size of a black hole (BH), i.e., $\alpha\equiv GM\mu_{s}/\hbar c^{3}\sim 1$, where $\alpha$ is the “gravitational fine-structure constant”, $M$ is the BH mass and $\mu_{s}$ is the boson mass, then a classical wave amplification process (superradiance) forms a bosonic cloud around the BH [54, 55, 56, 57, 58, 6, 17, 59, 60]. The formation of this cloud extracts rotational energy from the BH until the BH reaches a critical spin set by the Compton frequency of the boson and the mass of the BH [58, 6, 17, 59, 60]. This results on a critical spin curve on the BH mass-spin plane (e.g. Fig. 3 of Ref [17]): BHs born with spins above the critical spin curve will rapidly spin down until their final masses and spins lie on the curve. The net result of the superradiance process is thus to carve a region of the mass-spin plane, above the critical spin curve, where BHs are unlikely to be observed (“exclusion region”). Since the exact position and extent of the exclusion region depend on the boson mass (again, Fig. 3 of Ref [17]), one can use mass and spin measurements for a population of BHs to search and characterize ultralight bosons (e.g. [6, 17, 34, 35, 36, 46]).

Spin measurements of BHs in X-ray binaries (see e.g. References [61, 62]) could be used to search for bosons with their masses commensurate with BHs in the mass range $5~{}\mathrm{M}_{\odot}<M\lesssim 20~{}\mathrm{M}_{\odot}$ [61, 63]. Gravitational-wave (GW) signals emitted by binary black holes (BBHs) provide another avenue, as they encode the properties of their sources, including the masses and spins of the two component BHs. Since ground-based GW detectors such as LIGO and Virgo [64, 65] can detect heavier BHs ($M$ up to $\sim 100~{}\rm{M}_{\odot})$ [66, 67] than those found in X-ray binaries, the spin measurements inferred from GWs probe a lighter range of boson mass. This probe is complicated by two facts: (i) the measurements of individual BH spins with ground-based GW detectors are usually challenging [68, 69, 70, 71] and (ii) what is measured is the distribution of spins at merger, which is not only impacted by an eventual interaction with the bosons, but also on the distribution of spins at birth [72, 73, 74, 75].

In this paper we perform hierarchical Bayesian inference on a population of BBHs, to simultaneously infer the boson mass and the BH spin distribution at formation [76, 77, 78, 79, 80, 81, 82, 83, 84]. First, we show that the existence of ultralight bosons in the mass range between $10^{-13}$ eV and $3\times 10^{-12}~{}\mathrm{eV}$ can be ruled out with $\mathcal{O}(100)$ high signal-to-noise ratio (SNR) BBH detections. 111In the paper, we define “high SNR” as SNR$\geq 30$. Second, we illustrate how our method can confirm the existence of bosons with two examples of boson masses, $\mu_{s}=10^{-12}$ and $10^{-13}$ eV.

II Critical spin arising from superradiant instability

GW measurements yield the masses and spins of BHs at the time of merger. Therefore, one needs to account for the impact of superradiance on the evolution of spins, which we quickly review here. Superradiance causes the growth of a boson cloud in a time scale $\tau^{\rm inst}_{[nlm]}$, called instability time scale, and eventually spins down the host BH to a characteristic critical spin $\chi^{\infty}_{[nlm]}$ (which will be defined below) [58, 6, 17, 60]. The indices $[nlm]$ are the analog of the hydrogen atom quantum numbers, i.e., a set of radial, orbital azimuthal, and magnetic quantum numbers, of the cloud’s bound states. 222We follow Dolan’s notation, in which the ordinary principal quantum number $\tilde{n}=n+l+m$, such that $n=0$ corresponds to the dominant (nodeless) modes [58]. For any given $[nlm]$ mode, the instability time scale is representative of the exponential growth of the occupation number of the bosons in that mode, such that modes with small $\tau^{\rm inst}_{[nlm]}$ are populated quickly.

One often introduces the inverse of the instability time scale $\tau^{\mathrm{inst}}_{[nlm]}$, called the superradiant rate $\Gamma^{\mathrm{inst}}_{[nlm]}$, which can be analytically calculated for $\alpha\ll 1$ [85, 6, 17, 59, 60] (hereafter, we set $G=c=\hbar=1$),

[TABLE]

where $\chi$ is the dimensionless spin of the BH and $r_{+}\equiv M\left(1+\sqrt{1-\chi^{2}}\right)$ is the radial coordinate of the BH outer horizon. For the first few $l$’s, the fastest growing mode is the one with $l=m$, such that $\Gamma^{\mathrm{inst}}_{[nlm]}$ is the largest [17]. For any given $M$ and $\mu_{s}$, superradiance can happen as long as $\Gamma_{[nlm]}^{\mathrm{inst}}(\mu_{s},M,\chi)>0$; hence,

[TABLE]

which gives a condition on the BH spin $\chi$ such that superradiance can happen for the mode $[nlm]$. If the BH-cloud system had an infinite amount of time to evolve without disturbances, the spin angular momentum of the BH would eventually be lowered to the point where Eq. (2) cannot be satisfied or, equivalently, where $\Gamma_{[nlm]}^{\mathrm{inst}}=0$ (saturation of the superradiance). The critical spin for the $[nlm]$ mode is thus defined as the spin below which the superradiant growth of the cloud of the $[nlm]$ mode is forbidden,

[TABLE]

Besides the spin angular momentum of the BH, a small fraction $(\lesssim 10\%)$ of BH mass is also extracted and contributes to the cloud’s mass [17, 59, 86]. This is smaller or at most comparable to the BH mass uncertainty from GW measurements [70]. Therefore, we neglect the BH mass loss due to superradiance and assume that the BH masses at merger are the same as the masses at formation.

In this study, we do not allow for boson self-interaction, which would lead to additional phenomenology [17, 34, 87]. We refer to Ref. [46] for the analysis on self-interacting bosons using spin measurements of x-ray binaries.

III Postsuperradiant spin of astrophysical black holes

Astrophysical BHs in binaries have finite lifetimes, $\tau_{s}$, which implies their spins at merger will not reach the critical spin $\chi^{\infty}_{[nlm]}$. Therefore, we need to calculate the postsuperradiant BH spin for the $[nlm]$ mode, $\chi_{[nlm]}$ 333It is larger than $\chi^{\infty}_{[nlm]}$. by solving $\Gamma^{\rm inst}_{[nlm]}(\mu_{s},M,\chi_{[nlm]})=1/\tau_{s}$ for $\chi_{[nlm]}$, i.e., we truncate the spin evolution when the instability time scale decreases to the BH lifetime.

Naturally, if the lifetime of a BH is too short compared to the time required for the boson cloud to spin down its host BH, the effect might not even be measurable. Refs. [17, 34, 35] find that the bosonic field in the cloud should increase by $\sim 180$ $e$-foldings for the cloud to store a large amount of the BH angular momentum. This $e$-folding requirement translates to a “growth time scale” for the boson cloud such that it can significantly spin down the BH, $\tau^{\rm grow}_{[nlm]}\approx 180~{}\tau^{\rm inst}_{[nlm]}(\mu_{s},M,\chi_{I[nlm]})$, where $\chi_{I[nlm]}$ is the dimensionless spin of the BH at the onset of superradiance of the $[nlm]$ mode. Therefore, a BH that is born with a spin at formation $\chi_{F}$ and merges in $\tau_{s}$ can be spun down to the postsuperradiant spin $\chi_{[nlm]}$ (which is given by the solution of $\Gamma^{\rm inst}_{[nlm]}(\mu_{s},M,\chi_{[nlm]})=1/\tau_{s}$) only if $\chi_{F}>\chi_{[nlm]}$ and $\tau_{s}>\tau^{\rm growth}_{[nlm]}$.

Multiple clouds with different modes can be excited within the lifetime of a BH. The highest mode of the cloud that can be populated is given by the condition $\tau^{\rm grow}_{[nlm]}~{}{<}~{}\tau_{s}~{}{<}~{}\tau^{\rm grow}_{[n(l+1)(m+1)]}$. One can then estimate the BH spin at merger $\chi_{M}$ to be the postsuperradiant spin of this highest mode $\chi_{[nlm]}$, given the BH spin at formation and its lifetime. For example, if a BH is born with $\chi_{F}>\chi_{[011]}$, then the formation of the initial cloud (with the time scale $\tau^{\rm grow}_{[011]}$) slows down the BH spin to $\chi_{[011]}$. After its formation, the cloud dissipates away emitting nearly monochromatic gravitational waves [17, 35, 36, 37, 88]. Next, a second (“higher mode”) cloud is formed with the time scale $\tau^{\rm grow}_{[022]}$, and the BH spins down further to the postsuperradiant spin of this next mode $\chi_{[022]}$, if $\chi_{[011]}>\chi_{[022]}$. This cycle repeats until the BBH merges at time $\tau_{s}$. In Fig. 1, we show a schematic picture of a system for which both the $l=m=1$ and the $l=m=2$ clouds have enough time to form before merger. We note that only the first few growing modes $n=0$ and $l=m\leq 3$ are relevant for the typical astrophysical time scales $\tau_{s}\lesssim 10$ Gyr.

In principle, the BBH merger could cause the cloud to fall back to its host through level mixing [17, 86, 89, 90, 91]. While one could think that this results in a transfer of the cloud’s angular momentum back to the host BH, which would then be spun-up, most of the in-falling modes have nonpositive angular momentum (i.e. $m\leq 0$) due to selection rules. [89, 90, 91]. This is why recent studies have suggested that the in-falling cloud instead spins down the host BH, or might decrease the orbital angular momentum in the binary system [86, 89, 90, 92, 93, 91]. An eventual decrease in the BH spin by the fallback would further increase the size of the exclusion region on the BH mass-spin plane, making it easier to verify the existence of bosons with our method. To be conservative, we ignore this binary effect on the postsuperradiant spins.

IV Testing the Ultralight-boson Hypothesis using Hierarchical Bayesian inference

An astrophysical distribution of BH spins at birth, which produces mainly small spins in the absence of superradiance, is partially degenerate with the postsuperradiant spin distribution that originates from a moderate (or high) spin distribution at birth in the presence of superradiance. 444Although in the latter case one would expect a characteristic peak at around the postsuperradiant spin curve. Hence, we need to simultaneously infer the spin distribution at birth and the boson mass to properly account for this degeneracy.

We use hierarchical Bayesian inference [79, 78, 81, 82, 83], and consider two competing models: (i) in the “boson model,” $\mathcal{H}_{\mathrm{B}}$, we assume that a boson exists such that BHs can spin down to the corresponding postsuperradiant spins $\chi_{[nlm]}$ through superradiance (Sec. III); (ii) in the “astrophysical model,” $\mathcal{H}_{\rm{A}}$, ultralight bosons do not exist, and the spins of BHs merging in binaries are entirely determined by their astrophysical evolution.

The two hypotheses are distinguishable through the resulting distribution of the BH spins at merger. Specifically, for $\mathcal{H}_{\mathrm{B}}$, we assume

[TABLE]

where $\chi_{F}$ and $\chi_{M}$ are the individual BH spins at formation and at merger, respectively, and $\tau_{s}$ is the BH lifetime. The condition $\tau^{\rm grow}_{[nlm]}<\tau_{s}<\tau^{\rm grow}_{[n(l+1)(m+1)]}$ implies that $\chi_{[nlm]}$ is obtained as the postsuperradiance spin of the highest mode that can be populated within the binary lifetime. If superradiance does not happen for the specific BH (either because of its parameters at formation or because of the time scale to merger) then the spin at merger is the same as the spin at formation, i.e. $\chi_{M}=\chi_{F}$. On the other hand, there is no superradiance and therefore no spin evolution for $\mathcal{H}_{\rm{A}}$. Thus, for all black holes, the spin remains unchanged from formation to merger,

[TABLE]

For both models, we parametrize the distribution of BH spins at formation with a beta distribution, controlled by two unknown shape parameters $\alpha>0$ and $\beta>0$: $p(\chi_{F}|\alpha,\beta)\propto\chi_{F}^{\alpha-1}(1-\chi_{F})^{\beta-1}$. This is a generic functional form that can capture multiple different formation pathways [94, 80]. The boson model thus depends on three hyperparameters $\Lambda_{\mathcal{H}_{\mathrm{B}}}=(\alpha,\beta,\mu_{s})$, while the astrophysical model only has two hyperparameters $\Lambda_{\mathcal{H}_{\rm{A}}}=(\alpha,\beta)$. We aim at measuring the hyperparameters $\Lambda$, given a set of $N$ GW observations $\mathbf{d}=\left\{d_{k}\right\}$, whose morphology depends on a set of parameters $\theta$, such as BH masses, spins and distance [95]. The distribution of measured $\Lambda$, known as hyperposterior, can be written as [79, 78, 81, 82, 83]:

[TABLE]

where $p(\theta|\Lambda)$ is the expected distribution of the individual events parameters, given the hyperparameters; $\pi(\Lambda)$ are the priors of the hyperparameters; $p(d_{k}|\theta)$ is the likelihood of the $k$th GW source; and $\Sigma(\Lambda)$ is the normalization factor given by

[TABLE]

where $p_{\mathrm{det}}(\theta)$ is the detection probability for a BBH with parameters $\theta$. The normalization factor $\Sigma(\Lambda)$ can thus be interpreted as the fraction of detectable BBHs.

When working with the boson model, $\Lambda_{\mathcal{H}_{\mathrm{B}}}=(\alpha,\beta,\mu_{s})$ and $\theta_{\mathcal{H}_{\mathrm{B}}}=(M_{1},M_{2},\chi_{M,1},\chi_{M,2},\tau_{s})$, where $M_{i}$ and $\chi_{M,i}$ are the masses and spins (at merger) of the two compact objects in the binary. One thus has:

[TABLE]

In this expression, $\pi(M_{1},M_{2},\tau_{s})$ is the prior on the component masses and the merger time of BBHs, $\pi(\Lambda_{\mathcal{H}_{\mathrm{B}}})$ is the prior on the hyperparameters of the model $\mathcal{H}_{\mathrm{B}}$, and $p(\chi_{M,i}|\Lambda_{\mathcal{H}_{\mathrm{B}}},M_{i},\tau_{s})$ is the distribution of the spin magnitude at merger of $\mathcal{H}_{\mathrm{B}}$. It can be derived from the spin-magnitude distribution at formation as follows:

[TABLE]

in which we define the conditional probability

[TABLE]

where $\delta[...]$ is a Dirac delta that maps from the spin at formation to the spin at merger, and $\Theta[...]$ is a Heaviside step function that enforces the superradiance condition $\chi_{F,i}>\chi_{[nlm]}$. Since the growth time scale $\tau^{\rm grow}_{[nlm]}$ depends mildly on the initial spin at the onset of the superradiance, the postsuperradiant spin $\chi_{[nlm]}$ depends on four parameters $(\mu_{s},M,\tau_{s},\chi_{F})$ in principle. One needs to calculate $\chi_{[nlm]}$ for each $\chi_{F,i}$ to precisely evaluate the integral in Eq. (9). To facilitate the evaluation of this integral during the sampling process, we simplify the dependence of $\chi_{F,i}$ in $\chi_{[nlm]}$ with the following approximations to fix the values of $\chi_{I[nlm]}$. For the dominant $[011]$ mode, we set the spin at the onset to the midpoint between the minimum spin value required for superradiance and the maximal Kerr spin: $\chi_{I[011]}\approx(1+\chi_{[011]})/2$. This is justified because the growth time scale varies by only around 1 order of magnitude for spins in the range $\chi_{I[011]}\in[\chi_{[011]},1]$. After spin-down, the BH settles on the postsuperradiant spin of the given mode, $\chi_{[nlm]}$. Therefore, we approximate the initial spin of each subsequent mode $\chi_{[n(l+1)(m+1)]}$ by the preceding mode postsuperradiant spin, i.e., $\chi_{I[n(l+1)(m+1)]}\approx\chi_{[nlm]}$. With the above simplifications, $\chi_{[nlm]}$ only depends on $(\mu_{s},M_{i},\tau_{s})$ and we can approximate Eq. (9) as

[TABLE]

with $f_{\rm SR}$ being the differential fraction of BHs that undergo superradiance at the BH mass $M_{i}$:

[TABLE]

In the astrophysical model, $\mathcal{H}_{\rm{A}}$, one obtains a similar expression for $p(\Lambda_{\mathcal{H}_{\rm{A}}}|\mathbf{d},\mathcal{H}_{\rm{A}})$ by replacing $\mathcal{H}_{\mathrm{B}}\rightarrow\mathcal{H}_{\rm{A}}$ everywhere, and removing all references to $\tau_{s}$, which is not a relevant parameter of $\mathcal{H}_{\rm{A}}$.

While calculating Eq. (8), we use a power law prior $\pi(M_{1})\propto M_{1}^{-2.35}$ for the primary mass [96], and uniform prior for mass ratio $q=M_{2}/M_{1}$ with $M_{2}\leq M_{1}$. We also assume that $\tau_{s}$ is known (for the $\mathcal{H}_{\mathrm{B}}$ model) and fixed at $10\,\rm Myr$: $\pi(\tau_{s})=\delta(\tau_{s}-10\,{\rm Myr})$, which is toward the lower limit of the typical inspiral time scales, according to numerical simulations ($\sim 10\,{\rm Myr}-10\,{\rm Gyr})$ [97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107]. From one side, this choice is conservative as it allows for the least time for BHs to spin-down due to superradiance, making it harder to find evidence for bosons. From the other side, restricting the merger time prior overestimates the prior information thus overestimates the evidence for the boson hypothesis $\mathcal{H}_{B}$ in Eq. (8). Nevertheless, the additional parameter space in $\tau_{s}$ is expected to contribute modestly to the Bayes factor because the corrections to the available mass-spin parameter space are mostly smaller than the mass-spin measurement uncertainties. Therefore, while more realistic models for the merger time prior $\pi(\tau_{s})$ could be used, our choice is sufficient for advanced GW detectors, given their limited precision in the measurement of component masses and spins.

Since we assume the BH mass distribution is known, we do not need to consider its selection effect while calculating $\Sigma(\Lambda)$. We can also ignore selection effects due to BH spins as the expected number of observations only varies by $\lesssim 10\%$ for different spin models [76, 94, 108]. Furthermore, for $\Lambda_{B}$, the fraction of detectable BBHs does not depend on $\mu_{s}$. Based on the above arguments, we therefore assume $\Sigma(\Lambda_{\mathcal{H}_{\mathrm{B}}})$ and $\Sigma(\Lambda_{\mathcal{H}_{\rm{A}}})$ are constants in the evaluation of the hyperposteriors. In generating the simulations, however, we fully account for all selection effects so that the number of sources can be interpreted as the expected number of detections in future observations.

Integrating Eq. (8), and the equivalent expression for $\mathcal{H}_{\rm{A}}$, over the whole hyperparameters space yields evidences $Z_{\mathcal{H}_{\mathrm{B}}}$ and $Z_{\mathcal{H}_{\rm{A}}}$ that can be used to calculate the Bayes factor between the boson and astrophysical hypothesis: $\mathcal{B}_{\rm{A}}^{\rm{B}}={Z_{\mathcal{H}_{\mathrm{B}}}}/{Z_{\mathcal{H}_{\rm{A}}}}$. We also perform a Monte Carlo simulation of 50 different sets of sources in every simulated universe. This allows us to estimate the probability distribution of the Bayes factors due to Poisson fluctuation for each number of detections $N$.

V Mock data analysis

The method described above can be applied to both simulated and real detections. We first demonstrate its use on three different simulated “universes”: (i) one with a boson scalar field with $\mu_{s}=10^{-13}\,\rm eV$, (ii) one with a boson with $\mu_{s}=10^{-12}\,\rm eV$, and (iii) one where no boson exists (“astrophysical population”). To create the mock populations, we generate BBHs with component masses $M_{1,2}$ following the same prior in the model: $\pi(M_{1})\propto M_{1}^{-2.35}$, uniform distribution for $q$ in $[0.1,1]$ and require both $\{M_{1},M_{2}\}\in[5,50]\rm M_{\odot}$, consistently with Ref. [80]. The BBHs are distributed uniformly in the source-frame comoving volume, as well as the sky positions, orbital orientations and polarization angles in the unit sphere. The astrophysical processes that set the initial spin magnitude and orientation are still to be fully understood [72, 73, 74, 75]. For each of the three universes, we consider two distributions of formation spin magnitudes $\chi_{F,i}$: (a) uniform in $[0,1)$ (“flat spin”) and (b) $p(\chi_{F,i})\propto(1-\chi_{F,i})$ (“low spin”), with an isotropic spin orientation in both cases. The true shape parameters of the beta distribution are $\alpha=\beta=1$ and $\alpha=1,\beta=2$ for the “flat spin” and “low spin” populations, respectively.

When simulating the universes where bosons exist, we need to evolve the BH spins at formation to the spins at merger using Eq. (4). We assume all BBHs have a short merger time scale $\tau_{s}=10$ Myr, which minimizes the effect of superradiance and is thus a conservative choice. To keep the computational cost of the analysis reasonable, of all the sources we generate, we only analyze those for which SNR$>30$. These are the only sources that will contribute to the test since individual spins are hard to measure for low or medium SNR BBHs [68, 70]. The populations of synthetic BBH sources are thus added into simulated noise of the LIGO and Virgo detectors at design sensitivity [109, 110]. We use the LALInference [111, 112] algorithm with the IMRPhenomPv2 waveform family [113] to obtain posterior and likelihood distributions for the compact binary parameters of the simulated sources, which can be used to infer the population hyperparameters as described in the previous section. For all of the hyperparameters, we use uniform-in-log priors, with ranges $[0.01,10]$ for $\alpha$ and $\beta$, as well as $[10^{-13},3\times 10^{-12}]$ eV for $\mu_{s}$, which is the range of $\mu_{s}$ that can be realistically probed with ground-based GW detectors [17, 34, 35, 88].

In Fig. 2, we show the evolution of the log Bayes factor boson versus astrophysical model, $\log_{10}\mathcal{B}^{\rm B}_{\rm A}$ as more events are used for the test. The bottom $x$-axes show the numbers of loud events, while the top ones show the numbers of total events. 555Since the distribution of SNRs for BBH detected by advanced detectors is known analytically and goes as $P(\rho)\propto\rho^{-4}$ [114], one can calculate that there is one event with ${\rm SNR}>30$ for each 16 events with SNR$\geq 12$ on average.

All curves show that the underlying hypothesis is correctly preferred by the method, given enough number of observations. In Table 1, we report the expected numbers of observations required to significantly 666We follow Ref. [115] and strongly prefer the boson (astrophysical) hypothesis if $\mathcal{B}_{\rm{A}}^{\rm{B}}\geq 100$ ($\leq 0.01$). prefer a hypothesis for all pairs of the spin distribution and boson mass. In general, we would expect that fewer sources are required to disprove the boson hypothesis (first row of Table 1) than to confirm it. This is because even one highly spinning BH measurement can contradict $\mathcal{H}_{\mathrm{B}}$, whereas multiple BHs that match the predicted postsuperradiant spins are necessary to favor $\mathcal{H}_{\mathrm{B}}$. On the other hand, for some values of the boson mass, the morphology of the postsuperradiant spin distribution, including its dependence of the BH mass, may be very different from the astrophysical spin model, making it easier to prefer the boson hypothesis over the astrophysical hypothesis. This is, for example, the case for $\mu_{s}\sim 10^{-13}\,\rm eV$, whose exclusion region does not constrain BH masses below $\sim 15~{}\mathrm{M}_{\odot}$ (see Fig. 2 of  [53]). Unless the astrophysical spin distribution is significantly correlated with the BH mass (which does not seem to be the case based on the latest LVK results [116]), it is harder for the astrophysical model to match the expected postsuperradiant spin distribution, and hence is easier to verify the boson hypothesis for that boson mass. For heavier bosons, however, the resulting postsuperradiant spin distribution is similar to an astrophysical model (in the absence of bosons) with low BH spins at birth, which makes the two models harder to distinguish. This explains why smaller numbers of sources are required on average to confirm the existence of a boson with mass $\mu_{s}=10^{-13}\,\rm eV$ than to rule out bosons in the mass range $10^{-13}\text{ eV}\leq\mu_{s}\leq 3\times 10^{-12}\text{ eV}$, in the flat spin and low spin scenarios, Table 1.

We also notice that more sources are required for the test if BHs generally have low formation spins (“low spin” population) than for the “flat spin” population. This is expected since it is harder to prove the existence of a dearth of highly spinning BHs due to superradiant spin-down given a population with small formation spins.

Next, we look at the estimation of the individual hyperparameters. As an example, we take 300 high SNR detections drawn from the simulated Universe with “flat spin” and $\mu_{s}=10^{-12}$ eV. Figure 3 shows the corner plot of $(\alpha,\beta,\mu_{s})$ in $\mathrm{log}_{10}$ space, assuming $\mathcal{H}_{\mathrm{B}}$ (blue) or $\mathcal{H}_{\rm{A}}$ (orange). First, we notice the model $\mathcal{H}_{\mathrm{B}}$ results in a bimodal hyperposterior for $\mu_{s}$. This is because the exclusion region generated by the first superradiant mode of a boson with mass $\mu_{s}=10^{-12}$ eV is similar to the one generated by the second mode of a boson with roughly twice the mass. In turn, this implies at least a partial degeneracy between the two configurations. We note that the true $\mu_{s}$ is found at the primary peak, and the secondary peak becomes less prominent as the number of detections increases.

Second, using the astrophysical model $\mathcal{H}_{\rm{A}}$, we recover heavily biased values of $(\alpha,\beta)$, which control the shape of the spins at formation. To better visualize this bias, we recast the $(\alpha,\beta)$ hyperposteriors of both models into $p(\chi_{F}|\alpha,\beta)$, as shown in Fig. 4. The model $\mathcal{H}_{\rm{A}}$ (orange band) is indeed more consistent with the postsuperradiant spin distribution at merger $p(\chi_{M}|\alpha=1,\beta=1,\mu_{s}=10^{-12}\,\mathrm{eV})$ (black dashed line), instead of the spin distribution at formation $p(\chi_{F}|\alpha=1,\beta=1)$ (black solid line). This is not surprising, since $\mathcal{H}_{\rm{A}}$ cannot account for the superradiant spin loss and simply treats the spin at merger as if it were the spin at formation, i.e., $p(\chi_{M})$ as $p(\chi_{F})$. On the other hand, $\mathcal{H}_{\mathrm{B}}$ (blue band) can “undo” the superradiance and reconstruct $p(\chi_{F})$ much closer to the “true” distribution at formation (black solid line) in our simulation. Hence, both the hyperposterior and the $\chi_{F}$ distribution, inferred by the model $\mathcal{H}_{\mathrm{B}}$, are unbiased.

VI Discussion

In this paper, we have illustrated the use of hierarchical Bayesian inference to simultaneously measure the BH spin distribution at formation and to search for ultralight bosons, by combining mass and spin measurements from a population of BBHs detected by GW observatories. Our method relies on the morphology of the expected exclusion region in the BH mass-spin plane, which depends on the properties of the boson mass, to search for and characterize the boson. Applying this method on a mock data set, we have shown that BBHs discovered by ground-based GW detectors can be used to rule out the existence of ultralight bosons in the mass range $[10^{-13},3\times 10^{-12}]$ eV. Our method is also capable to reveal the existence of a ultralight boson, and measure its mass, as we have explicitly shown with two example of boson masses, $\mu_{s}=10^{-13}$ eV and $10^{-12}$ eV. In order to investigate the impact of spin distribution at formation on the statistical power of our method, we have generated populations of simulated BBHs with either a uniform distribution or a linearly decreasing distribution of BH spins at formation. We found that in both cases combining $\lesssim 300$ high SNR events will be enough to rule out or confirm the existence of an ultralight boson within the mass range $[10^{-13},10^{-12}]$ eV.

While we only consider scalar bosons in this study, the method we developed is applicable to vector or tensor boson fields, which have much shorter instability and GW emission time scales [36, 38, 117]. Our analysis of simulated BBHs has made a few simplifying assumptions which make it conservative. First, we have assumed that all BBHs merge in 10 Myr, which is toward the lower limit of what is usually obtained in numerical simulations [97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107]. Since most of the BBHs in our parameters space of interest would have undergone superradiance within 10 Myr except for the very low mass systems $\sim 5M_{\odot}$, assuming longer merger times does not significantly improve the searching efficiency. Second, we have assumed that only sources with SNR$>30$ will contribute to this test, as their spins are easier to measure. In reality, while the component spins of weaker events are harder to measure, they will still contribute to the test.

In the analysis, we ignored the $\lesssim 10\%$ BH mass loss during superradiance. In order to assess the impact of this choice, we repeated the analysis with all BH mass posteriors shifted by $5\%$ toward the light side, hence mimicking the effect of mass loss. This translates to a systematic overestimation of $\sim 5\%$ boson mass, which is still within the statistical uncertainty, in our inference with $\mathcal{O}(100)$ high SNR sources. However, we expect this systematic error will dominate when the number of events grow to $\mathcal{O}(10^{4})$ in the era of next-generation detectors [118, 119, 120, 121, 122, 123, 124, 125].

The true distribution of spins at formation plays the most important role: the number of events needed to perform this test will be larger if the astrophysical distribution of spins at formation is such that small spins are preferred. Conversely, if many highly spinning BHs are formed, potentially with significant misalignments between spin and angular momenta (both of which make spins easier to measure), then fewer sources will be necessary. Given the measured BBH merger rate, ground-based interferometers will detect hundreds of BBHs per year at design sensitivity [126, 127, 128, 129, 130]. In this large-number observations regime, one will want to use more sophisticated models which also capture eventual correlations between the masses and spins of astrophysical BHs [72, 73, 75]. This may boost or suppress the statistical power of testing boson hypothesis in a different boson mass range, depending on the actual joint distribution of BH mass and spin at formation. We will leave investigating these systematics to future work.

Within the assumptions made in this study, it seems feasible to rule out the existence of ultralight bosons everywhere in the mass range of $[10^{-13},3\times 10^{-12}]$ eV with a few years of advanced detectors data. However, we note that fewer GW events would be required to rule out the bosons in a narrower boson mass range. Statistically proving the existence of these bosons will take longer, as more sources are required: the planned upgrades of LIGO and Virgo to their “plus” configurations might yield thousands of BBH events per year, which will make it more plausible to gather the evidence for the existence of ultralight bosons in the mass range of $[10^{-13},3\times 10^{-12}]$ eV [119, 121].

VII Acknowledgements

We thank the anonymous referees for their suggestions which significantly improved this paper. We also thank Emanuele Berti, Richard Brito, Will Farr, Carl Haster, Max Isi, and Kaze Wong for the valuable discussions and suggestions. K. K. Y. N. and S. V. acknowledge the support of the National Science Foundation (NSF) through the NSF Grant No. PHY-1836814. LIGO was constructed by the California Institute of Technology and Massachusetts Institute of Technology with funding from the National Science Foundation and operates under cooperative agreement PHY-1764464. The work of O. A. H. was supported by the Hong Kong Ph.D. Fellowship Scheme issued by the Research Grants Council of Hong Kong before resubmission and was supported by the research program of the Netherlands Organization for Scientific Research. The work of T. G. F. L. was partially supported by grants from the Research Grants Council of Hong Kong (Projects No. CUHK14306218, No. CUHK14310816, and No. CUHK24304317), Research Committee of the Chinese University of Hong Kong, and the Croucher Foundation in Hong Kong. This research made use of data, software, and/or web tools obtained from the Gravitational Wave Open Science Center [131], a service of LIGO Laboratory, the LIGO Scientific Collaboration, and the Virgo Collaboration. The authors are grateful for computational resources provided by the LIGO Lab and supported by the NSF Grants No. PHY-0757058 and No. PHY-0823459.