A first search for a stochastic gravitational-wave background from ultralight bosons

TL;DR

This paper develops a Bayesian analysis framework to search for a stochastic gravitational-wave background from ultralight bosons using LIGO and Virgo data, aiming to constrain boson masses and distinguish signals from other sources.

Contribution

It introduces an improved Bayesian pipeline incorporating binary merger remnants to analyze gravitational-wave data for ultralight bosons, and applies it to real LIGO data.

Findings

No evidence of the bosonic cloud signal was found.

Boson masses between 2.0e-13 eV and 3.8e-13 eV are excluded under optimistic assumptions.

The analysis highlights degeneracies affecting robust exclusion of boson presence.

Abstract

In this work, we develop a Bayesian data analysis framework to study the SGWB from bosonic clouds using data from Advanced LIGO and Advanced Virgo, building on previous work by Brito et.al (2017). We further improve this model by adding a BH population of binary merger remnants. To assess the performance of our pipeline, we quantify the range of boson masses that can be constrained by Advanced LIGO and Advanced Virgo measurements at design sensitivity. Furthermore, we explore our capability to distinguish an ultralight boson SGWB from a stochastic signal due to distant compact binary coalescences (CBC). Finally, we present results of a search for the SGWB from bosonic clouds using data from Advanced LIGO's first observing run. We find no evidence of such a signal. Due to degeneracies between the boson mass and unknown astrophysical quantities such as the distribution of isolated BH spins, our analysis cannot robustly exclude the presence of a bosonic field at any mass. Nevertheless, we show that under optimistic assumptions about the BH formation rate and spin distribution, boson masses in the range $ \SI{2.0e-13}{eV}\leq m_\mathrm{b}\leq\SI{3.8e-13}{eV} $ are excluded at 95% credibility, although with less optimistic spin distributions, no masses can be excluded. The framework established here can be used to learn about the nature of fundamental bosonic fields with future gravitational wave observations.

Tables

Table 1. Table 1: Parameters in each recovered background model.

Models	CBC	SI$+$CBC
Parameters	${\Omega }_0,\alpha$	$m_b,{\chi }_{\mathrm{ul}},{\Omega }_0,\alpha$

Equations

$$m_b{c}^2$$

$$0<\omega <m{\Omega }_H,$$

$${\Omega }_H=\frac{\chi }{2r_{+}}.$$

$$\mu \equiv \frac{m_b}{\hslash },$$

$$\begin{array}{rl}{\omega }_{lmn}& \equiv {\omega }_R+i{\omega }_I\\ & \simeq \mu +i2\gamma \mu r_{+}(m{\Omega }_H-\mu )(M_i\mu {)}^{4l+4},\end{array}$$

$$\displaystyle\tau_{\mathrm{inst}}\sim 0.07\chi^{-1}\left(\frac{M_{i}}{$10\text{\,}\textit{M}_{\odot}$}\right)\left(\frac{0.1}{M_{i}\mu}\right)^{9}\mathrm{yr}.$$

$$\frac{\mathrm{d}{E}_{\mathrm{GW}}}{\mathrm{d}t}\approx \frac{484+9{\pi }^2}{23040}{\left(\frac{M_S}{M_f}\right)}^2(M_f\mu {)}^{14},$$

$$\displaystyle\tau_{\mathrm{GW}}\sim 6\times 10^{4}\chi^{-1}\left(\frac{M_{f}}{$10\text{\,}\textit{M}_{\odot}$}\right)\left(\frac{0.1}{M_{f}\mu}\right)^{15}\mathrm{yr}.$$

$$J_f=J_i-\frac{m}{{\omega }_R}(M_i-M_f).$$

$$M_f=\frac{m^3-\sqrt{m^6-16m^2{\omega }_R^2(m{M}_i-{\omega }_R{J}_i{)}^2}}{8{\omega }_R^2(m{M}_i-{\omega }_R{J}_i)}.$$

$$M_S^{\max }=M_i-M_f\approx 0.1M_i\chi \left(\frac{M_i\mu }{0.1}\right).$$

$$E_{\mathrm{GW}}$$

$${\Omega }_{\mathrm{GW}}(f)\equiv \frac{1}{{\rho }_c}\frac{\mathrm{d}{\rho }_{\mathrm{GW}}}{\mathrm{d}\ln (f)},$$

$${\Omega }_{\mathrm{GW}}(f)=\frac{f}{{\rho }_c}\int \mathrm{d}z\frac{\mathrm{d}t}{\mathrm{d}z}\int \mathrm{d}\mathbf{\theta }p(\mathbf{\theta })R(z;\mathbf{\theta })\frac{\mathrm{d}{E}_s}{\mathrm{d}{f}_s}(\mathbf{\theta };f(1+z)),$$

$$\frac{\mathrm{d}t}{\mathrm{d}z}=\frac{1}{(1+z)H_0\sqrt{{\Omega }_M(1+z{)}^3+{\Omega }_{\Lambda }}},$$

$$\frac{\mathrm{d}{E}_{\mathrm{s}}}{\mathrm{d}{f}_s}\approx E_{\mathrm{GW}}\delta (f(1+z)-f_0),$$

$${\Omega }_{\mathrm{GW}}^{\mathrm{iso}}(f)=\frac{f}{{\rho }_c}\int \mathrm{d}z\frac{\mathrm{d}t}{\mathrm{d}z}\int \mathrm{d}M\mathrm{d}\chi p(\chi )\frac{\mathrm{d}\dot{n}}{\mathrm{d}M}\frac{\mathrm{d}{E}_s}{\mathrm{d}{f}_s},$$

$$p(\chi )=\{\begin{array}{ll}0& (\chi <{\chi }_{\mathrm{ll}},{\chi }_{\mathrm{ul}}<\chi )\\ \frac{1}{{\chi }_{\mathrm{ul}}-{\chi }_{\mathrm{ll}}}& ({\chi }_{\mathrm{ll}}\le \chi \le {\chi }_{\mathrm{ul}}),\end{array}$$

$$\frac{\mathrm{d}\dot{n}}{\mathrm{d}M}$$

$$=\psi (z_f)\frac{\xi ({\mathcal{M}}_{\ast })}{{\mathcal{M}}_{\ast }}\frac{\mathrm{d}{\mathcal{M}}_{\ast }}{\mathrm{d}M}.$$

$$\begin{array}{r}\psi (z)=\nu \frac{a\exp (b(z-z_m))}{a-b+b\exp (a(z-z_m))},\\ a=2.37,b=1.80,\nu =0.178,z_m=\mathrm{2.00.}\end{array}$$

$$\xi ({\mathcal{M}}_{\ast })=\frac{{\mathcal{M}}_{\ast }^{-1.35}}{{\int }_{0.1M_{\odot }}^{100M_{\odot }}{\mathcal{M}}_{\ast }^{-1.35}\mathrm{d}{\mathcal{M}}_{\ast }}.$$

$$\frac{\mathrm{d}{\mathcal{M}}_{\ast }}{\mathrm{d}M}={\left(\frac{\mathrm{d}g({\mathcal{M}}_{\ast },Z)}{\mathrm{d}{\mathcal{M}}_{\ast }}\right)}^{-1}.$$

$$\begin{array}{rl}{\Omega }_{\mathrm{GW}}^{\mathrm{rem}}(f)=\frac{f}{{\rho }_c}& \int \mathrm{d}z\frac{\mathrm{d}t}{\mathrm{d}z}\\ \times & \int \mathrm{d}\mathbf{m}\mathrm{d}\chi p(\mathbf{m})R_m(z;\mathbf{m})p(\chi )\frac{\mathrm{d}{E}_s}{\mathrm{d}{f}_s}.\end{array}$$

$$p(\chi )=\delta (\chi -0.7).$$

$$p(m_1,m_2)\propto \frac{m_1^{-2.35}}{m_1-5M_{\odot }},$$

$$M\approx m_1+m_2-5.7\times 10^{-2}\frac{m_1{m}_2}{m_1+m_2},$$

$$R_m(z;\mathbf{m})={\int }_{t_{\min }}^{t_{\max }}{R}_f(z_f;\mathbf{m})p(t_d)\mathrm{d}{t}_d,$$

$$R_f(z_f;\mathbf{m})\propto \{\begin{array}{ll}\psi (z_f)& \mathrm{if}{m}_1,m_2<30M_{\odot }\\ e(z_f)\psi (z_f)& \mathrm{otherwise}.\end{array}$$

$$\int p(\mathbf{m})R_m(z=0;\mathbf{m})\mathrm{d}\mathbf{m}=103{\mathrm{Gpc}}^{-3}{\mathrm{yr}}^{-1}.$$

$${\Omega }_{\mathrm{GW}}(f)={\Omega }_{\mathrm{GW}}^{\mathrm{iso}}(f)+{\Omega }_{\mathrm{GW}}^{\mathrm{rem}}(f),$$

Full text

gravitational waves, stochastic backgrounds, superradiant instability, Bayesian inference

I Introduction

The first detection of gravitational waves (GWs) from binary black hole (BBH) Abbott et al. (2016a) and neutron star (BNS) Abbott et al. (2017a) coalescences by the Advanced Laser Interometer Gravitational-wave Observatory (LIGO) Aasi et al. (2015a) and Advanced Virgo Acernese et al. (2015) represented a historical breakthrough, creating an alternative window through which to view the Universe. Furthermore, subsequent observations of BBH coalescence events have firmly established GW astronomy Abbott et al. (2016b, 2017b, 2017c, 2017d, 2018a), shedding new light on fundamental physics Abbott et al. (2016c) and the nature of the stellar mass black hole (BH) population Abbott et al. (2016d, 2018b). In the near future, Advanced LIGO and Advanced Virgo will be joined by additional detectors, like KAGRA Aso et al. (2013) and LIGO-India Iyer et al. (2011). A major target for the advanced detector network is the stochastic gravitational-wave background (SGWB), a superposition of many sources too faint to resolve individually Abbott et al. (2016e, 2017e).

Ultralight bosons around spinning BHs have been proposed as a possible source of the SGWB Brito et al. (2017); Brito (2017). Superradiant instabilities induced by the bosonic fields co-rotating around BHs result in angular momentum extraction analogous to the so-called Penrose process Penrose (1969). This instability not only spins down the BH, but also triggers an exponential growth of the bosonic field. Subsequently the resulting bosonic cloud exhibits a time-dependent quadrupole moment, leading to GW radiation. We could obeserve these GWs with current ground-based GW detectors in two regimes – the “resolvable” regime in which nearby sources within $\mathcal{O}(10\,{\rm Mpc})$ can be directly detected Brito et al. (2017); Palomba et al. (2012); D’Antonio et al. (2018); Isi et al. (2018) and the “unresolvable” regime where a superposition of all other sources in the Universe will contribute to a SGWB Brito et al. (2017); Brito (2017).

Ground-based detectors are sensitive to bosons with a given mass scale that can be determined through dimensional-analysis Brito et al. (2017); Brito (2017). We expect a bosonic field with the mass $m_{b}$ to couple strongly to BHs whose Schwarzschild radius $r_{H}$ is comparable to the boson’s Compton wavelength $\lambda_{C}=\hbar/m_{b}c$, where $\hbar$ is the reduced Planck constant and $c$ is the speed of light. For BHs of around $10M_{\odot}$, this implies

[TABLE]

where $G$ is the gravitational constant and $M$ is the mass of a BH. In other words, only BHs within a relatively narrow window will significantly couple to bosonic field with a certain mass scale. This can be physically interpreted as follows. For a given boson mass, if BH mass is too large or small, the timescale of superradiant instability increases and its effect is exponentially suppressed Brito et al. (2017, 2015a); Dolan (2007). As we will demonstrate, Advanced LIGO could detect superradiant instabilities surrounding $\sim 10M_{\odot}$ BHs, and therefore can probe boson masses near $m_{b}\approx${10}^{-13}\text{,}\mathrm{eV}$$.

Light bosonic fields, such as QCD axions or axion-like particles have attracted attention, as a promising candidate of dark matter Arvanitaki (2010); Essig et al. (2013). It is challenging to use traditional particle physics experiments to probe the existence of light particles that do not couple strongly to ordinary matter, but we can use GW observations to study particle physics beyond the Standard Model Arvanitaki (2015, 2017).

Since angular momentum extraction is followed by characteristic GW emission, the detection of this class of GWs could also provide an explanation for the possible abundance of low spin BHs Brito et al. (2015b), which is consistent with the current BH spin measurement by Advanced LIGO Abbott et al. (2016f, 2017c, 2017d, 2018b). Additionally, the BH spin distribution is considered to be a crucial indicator of different BH formation channels Belczynski et al. (2016); Tutukov and YungelSon (1993); Sigurdsson and Hernquist (1993); Rodriguez et al. (2016). The discovery of ultralight bosons regulating BH spins would have significant implications for understanding BH formation scenarios.

Previous work has studied the SGWB from bosonic clouds Fan (2018); Brito (2017). Ref. Fan (2018) examined the detectability of a stochastic background from BHs spinning down, assuming the current ground-based detector network operates at design sensitivity. They assumed two different signal spectra, referred to as the Gaussian model and the quasi-normal-mode model, and using a Fisher analysis they also assessed the capability of extracting the component of their signal models from the total background including the CBC background. Ref. Brito (2017) computed the background spectrum expected from the superradiant instability assuming only an isolated BH population model. They also compared the predicted spectrum with the power-law integrated curves Thrane (2013) of Advanced LIGO to study its detectability. In this work, we will develop a Bayesian analysis framework to search for the SGWB from the superradiant instability and apply it to data from Advanced LIGO’s first observing run. We construct our signal model based on Brito et al. (2017); Brito (2017), additionally incorporating the expected contribution from remnant BHs produced by compact binary mergers to obtain a more accurate prediction.

This paper is structured as follows. In Section II we review the model of superradiant instabilities as previously computed in Brito et al. (2017); Brito (2017); Brito et al. (2015a). In Section III, we present the predicted SGWB signal, including the contribution from population of binary merger remnants. In Section IV, we present a Bayesian framework with which to search for this background. In Section V, we discuss the sensitive range in the boson mass and the model selection capability between this background and the projected SGWB that arises due to unresolved CBCs. Furthermore, in Section VI we show the range of excluded boson masses using data from Advanced LIGO’s first observing run. Finally, Section VII summarizes the results and future work.

In what follows, all quantities are described in units $G=c=1$.

II Superradiant instability

Before describing the data analysis methods, first we briefly review the theoretical basis of the bosonic cloud model as described in Brito et al. (2017); Brito (2017). For a review of supperradiance, see Brito et al. (2015b). Hereafter, we will restrict ourselves to scalar (spin-0) fields. Although vector (spin-1) East (2017); Pani (2012); Baryakhtar et al. (2017) and tensor (spin-2) Brito (2013) fields can additionally induce superradiant instabilities, the scalar case has been investigated most extensively Arvanitaki (2015, 2017); Dolan (2007); Yoshino and Kodama (2014); Brito et al. (2015b).

Incident radiation can by amplified by extracting energy and angular momentum from a BH. As a result, the mass and angular momentum of the BH transitions from the initial states $M_{i},J_{i}$ to the final states $M_{f},J_{f}$, respectively. Denoting the frequency of the radiation by $\omega$ and letting $l,m$ be the quantum numbers for angular momentum, the condition for superradiance to occur is given by Brito et al. (2015b)

[TABLE]

where $\Omega_{H}$ is the horizon angular velocity of the BH, given by

[TABLE]

Here, $\chi$ is the dimensionless spin parameter and $r_{+}$ is the outer event horizon of a Kerr BH.

Using black hole perturbation theory, in the limit that the scalar field’s Compton wavelength is much larger than the radius of the BH ($M_{i}\mu\ll 1$), where

[TABLE]

an approximate expression can be found for the eigenfrequencies of the scalar field Dolan (2007); Brito et al. (2017):

[TABLE]

where $n,l,m$ are the quantum numbers for energy and angular momentum for a wave function in a spherical potential, $\gamma$ is a positive numerical factor depending on $n,l,m$, and $\omega_{R/I}$ are real quantities. The imaginary part of the eigenfrequency, $\omega_{I}$, will be positive, indicating that the mode is (linearly) unstable, whenever Eq. (2) is satisfied111Note that we have used the relation $\omega_{R}=\mu$, which is valid to leading order in $M_{i}\mu$.. This instability leads to the formation of a scalar “cloud” around the BH. We define the instability timescale $\tau_{\mathrm{inst}}\equiv\omega_{I}^{-1}$ as the characteristic timescale over which the cloud forms. We will focus only on the $n=0$, $l=m=1$ mode hereafter since this mode has the shortest timescale and dominates the gravitational-wave radiation. In the large Compton wavelength limit, the instability timescale is given by Brito (2017); Brito et al. (2015a, 2017)

[TABLE]

As illustrated in Figure 1, after the scalar cloud has been formed, the scalar field profile will source a stress-energy tensor with a time-dependent quadrupole moment. The energy stored in the rotating scalar field configuration will emit gravitational-waves, with a frequency of $\omega_{\rm GW}=2\omega_{R}$ and a luminosity given by Brito (2017); Brito et al. (2015a, 2017)

[TABLE]

where $M_{S}$ is the total mass (energy) of the amplified scalar field configuration. Eventually, the scalar cloud will lose energy due to gravitational radiation. This process will happen over the course of the emission timescale $\tau_{\rm GW}$, given by Brito (2017); Brito et al. (2015a, 2017)

[TABLE]

Crucially, there is a significant difference between the instability and emission timescales, $\tau_{\mathrm{GW}}/\tau_{\mathrm{inst}}\sim 10^{6}$. Thus, we make the approximation that the scalar cloud forms more or less instantaneously and then starts to emit GWs. Neglecting the emission of GWs during the growth of the scalar cloud, conservation of energy and angular momentum implies that Brito (2017)

[TABLE]

The amplification of the cloud stops when the condition $\omega_{R}=m\Omega_{H}$ is satisfied. Combining this with Eq. (9), we obtain the final mass of the BH

[TABLE]

For the lowest mode $m=1$, in the long Compton wavelength limit $M_{i}\mu\ll 1$, the maximum mass of the scalar cloud reads

[TABLE]

Thus, for typical parameters $M_{i}\mu\approx 0.1$, about 10% of the energy of the initial BH is stored in the scalar cloud configuration.

Finally, we can evaluate the total GW energy emitted between the superradiance saturation and the present time as follows

[TABLE]

where $\tau_{\mathrm{GW}}$ is the GW emission timescale. Here, the signal duration $\Delta t$ is reasonably defined as $\Delta t\equiv\min(\tau_{\mathrm{GW}},t_{0})$, in which $t_{0}$ is the age of the Universe $\approx$13.8\text{,}\mathrm{G}\mathrm{y}\mathrm{r}$$. The SGWB will be calculated by summing the energy emitted by each source over the population of rotating BHs.

III Modeling the SGWB

Now we turn to modelling the SGWB produced by scalar clouds which form around BHs. We consider isolated BHs, following the model of Brito et al. (2017). We also derive the contribution due to remnant BHs formed from compact binary coalescences, and verify that it is smaller than the contribution of isolated BHs for most scalar masses of interest.

III.1 Framework to compute the SGWB

The SGWB is described by its energy density spectrum Allen and Romano (1999); Romano and Cornish (2017), defined as

[TABLE]

where $\rho_{\mathrm{GW}}$ is the energy density of GWs existing in the Universe and $\rho_{c}$ is the critical energy density required to have a spatially flat Universe. The astrophysical SGWB has been studied for a wide array of sources Kowalska-Leszczynska et al. (2015); Rosado (2011); Marassi (2011); Wu (2012); Ölmez (2010); Jenkins and Sakellariadou (2018) and is generally given by

[TABLE]

where $\,\mathrm{d}E_{s}/\,\mathrm{d}f_{s}$ is the source-frame energy spectrum of an individual astrophysical event, $R(z;\bm{\theta})$ is the number of sources per unit comoving volume per unit source frame time, and $p(\bm{\theta})$ is the multivariate probability distribution of the source parameters $\bm{\theta}$. Note that the source frame frequency $f_{s}$ is related to the observed frequency $f$ by a factor of the redshift: $f_{s}=(1+z)f$. The function ${\rm d}t/{\rm d}z$ is determined by standard cosmology

[TABLE]

where $\Omega_{M}$ and $\Omega_{\Lambda}$ are the dimensionless matter density and the dimensionless cosmological constant density, respectively. We use the cosmological parameters inferred from Planck Ade et al. (2016), i.e. $H_{0}=$68\text{,}\mathrm{k}\mathrm{m}~{}\mathrm{s}^{-1}\mathrm{M}\mathrm{p}\mathrm{c}^{-1}$$ and $\Omega_{M}=1-\Omega_{\Lambda}=0.308$ .

When constructing a model for the SGWB from a superradiant instability, we take the source parameters $\bm{\theta}$ to be BH masses and their dimensionless initial spin parameter $\chi$. This $\chi$ represents the BH spin before spinning down due to the instability. We note that the GWs emitted from an individual system actually have a slight frequency drift caused by the change in the cloud’s binding energy during the GW emission Isi et al. (2018). However, its overall drift is well within the width of each frequency bin that we use for our analysis. This leads to the following simple spectrum for a single source

[TABLE]

where $f_{0}=\omega_{R}/\pi\approx\mu/\pi$ (see Eq. (5)) and $E_{\mathrm{GW}}$ is given by Eq. (12).

In order to compute the background due to unresolved sources, we only integrate over the range of parameters $\theta$ that produces a signal-to-noise ratio (SNR) of less than 8 in a typical search for the astrophysical system in question D’Antonio et al. (2018); Isi et al. (2018). Nevertheless, the specifics of this cutoff do not change the overall shape of the predicted spectra significantly.

III.2 BH population models

To search for and constrain the stochastic background from a superradiant instability, we will need to assume a specific source number density and the mass and spin distributions of BHs. Below, we will consider two possible BH populations: (i) isolated BHs formed by core-collapse supernovae (CCSNe) and (ii) BBH merger remnants.

III.2.1 Isolated black holes

Here we assume a population of isolated BHs born from CCSNe. Eq. (14) can be reorganized as

[TABLE]

where ${\,\mathrm{d}\dot{n}}/{\,\mathrm{d}M}$ is the source-frame BH formation rate per comoving volume per BH mass $M$. For the probability density of $\chi$, we assume a uniform distribution

[TABLE]

where $\chi_{\mathrm{ll}},\chi_{\mathrm{ul}}$ are the lower and upper limit of the distribution. Inspired by Brito (2017), we adopt two different parameterizations of $p(\chi)$, either (a) varying the lower limit $\chi_{\mathrm{ll}}$ and fixing $\chi_{\mathrm{ul}}=1$, or (b) varying the upper limit $\chi_{\mathrm{ul}}$ and fixing $\chi_{\mathrm{ll}}=0$. Note that the first case is more optimistic than the second case, as it ensures a population of high-spin BHs that readily yield superradiant instabilities. Although these models for $p(\chi)$ admittedly quite simple, the true spin distribution of isolated BHs is extremely uncertain. As we will show, different $\chi_{\mathrm{ll,ul}}$ values crucially affect the background spectrum.

The BH formation rate $\,\mathrm{d}\dot{n}/\,\mathrm{d}M$ reads

[TABLE]

$\mathcal{M}_{*}$ is the mass of the BH’s progenitor star, whose population properties follow the cosmic star formation rate (SFR) $\psi(z)$ and initial mass function (IMF) $\xi(\mathcal{M}_{*})$. We adopt the SFR model proposed in Vangioni et al. (2015),

[TABLE]

$z_{f}$ is the redshift at the progenitor’s time of birth, that is, $t_{f}=t-\tau(\mathcal{M}_{*})$, where $\tau(\mathcal{M}_{*})$ is the lifetime of a progenitor based on Schaerer, D. (2002). Note that $\xi(\mathcal{M})$ is the IMF defined in terms of stellar mass fraction. This implies that $\xi(\mathcal{M}_{*})\,\mathrm{d}\mathcal{M}_{*}$ yields the ratio of the total mass of stars whose mass lies between $\mathcal{M}_{*}\sim\mathcal{M}_{*}+\,\mathrm{d}\mathcal{M}_{*}$ to the whole stellar mass. Since the Salpeter IMF is chosen as $\xi(\mathcal{M}_{*})$ in the mass range $\mathcal{M}_{*}\in[0.1\sim 100]~{}M_{\odot}$, its expression is222One should not confuse it with the other definition, that is, the number fraction of stars whose mass lies between $\mathcal{M}_{*}\sim\mathcal{M}_{*}+\,\mathrm{d}\mathcal{M}_{*}$. Let this be $\phi(\mathcal{M}_{*})$. The Salpeter function in this definition follows $\phi(\mathcal{M}_{*})\propto\mathcal{M}_{*}^{-2.35}$.

[TABLE]

The BH mass is related to the mass and metalliticity of the progenitor star via $M=g(\mathcal{M}_{*},Z)$. In this work we use a numerical fit for $g(\mathcal{M}_{*},Z)$ given by Fryer et al. (2012). It follows that

[TABLE]

Note that $g(\mathcal{M},Z)$ implicitly depends on redshift via the stellar metallicity $Z$, whose evolution over the cosmic history is evaluated in Ma et al. (2016). The lower BH mass cutoff of $g(\mathcal{M},Z)$ is set as $3\text{\,}\textit{M}_{\odot}$.

III.2.2 Binary black hole merger remnants

We also consider remnant BHs formed by BBH mergers, whose background spectrum is evaluated such that

[TABLE]

$\bm{m}$ represents a set of component masses of BBHs, i.e. the primary mass $m_{1}$ and the secondary mass $m_{2}$ ($m_{1}>m_{2}$ by definition) and $\,\mathrm{d}\bm{m}\equiv\,\mathrm{d}m_{1}\,\mathrm{d}m_{2}$. $R_{m}(z;\bm{m})$ is the BBH merger rate density for a given $\bm{m}$ and $z$. Unlike the spin distribution of isolated BHs, here $p(\chi)$ can be motivated by both the spin measurement of final remnant BHs by Advanced LIGO and Virgo Abbott et al. (2016f, 2017b, 2017d, 2017c) and several numerical simulations for mergers of similar-mass BHs González (2007); Berti (2007); Rezzolla et al. (2008). Both suggest that the spin magnitude of final remnant BHs is around $0.7$. Therefore, we assume all the remnant BHs initially have $\chi=0.7$, namely

[TABLE]

Following Abbott et al. (2017e, f, 2016f), we adopt the BH mass distribution

[TABLE]

with the constraints that $5M_{\odot}<m_{1},m_{2}<95M_{\odot}$ and $m_{1}+m_{2}<100M_{\odot}$. We approximate the mass of remnant BHs as

[TABLE]

where the binding energy of the inner most stable orbit in a binary system is subtracted from the total mass Maggiore (2007).

We evaluate the merger rate as described in Dominik et al. (2013); Abbott et al. (2016e, 2017e, 2017f),

[TABLE]

where $t_{d}$ is the time delay between a binary formation and its merger and $p(t_{d})$ is the time delay distribution. $z_{f}$ is the redshift at the binary formation time $t_{f}=t(z)-t_{d}$, in which $t(z)$ is the cosmic time at merger. Following Abbott et al. (2016e, 2017e), we assume that: (i) the distribution $p(t_{d})$ is taken as $p(t_{d})\propto 1/t_{d}$ in the range $t_{\mathrm{min}}<t_{d}<t_{\mathrm{max}}$, where $t_{\mathrm{min}}$ is $50~{}\mathrm{Myr}$ Abbott et al. (2016e); Dominik et al. (2013) and $t_{\mathrm{max}}$ is the Hubble time Nakar (2007); O’Shaughnessy et al. (2008); Dominik et al. (2012, 2013), (ii) the formation rate of BBHs, both of which are lighter than $30M_{\odot}$, evolves proportionally to the SFR $\psi(z)$, (iii) the SFR of Vangioni et al. (2015) is consistent with the one used for the isolated BHs model (see Eq. (21)), (iv) massive binaries in which at least one of the two component masses is above $30M_{\odot}$ cannot be formed in a high metallicity environment where $Z>Z_{\odot}/2$ ($Z_{\odot}$ is the solar metallicity). For the assumption (iv), we instead multiply the SFR by a weighting factor $e(z)$ to account for the fraction of star formation in an environment where $Z\leq Z_{\odot}/2$ Madau and Dickinson (2014). In other words,

[TABLE]

After marginalizing Eq. (28) over the component mass distribution $p(\bm{m})$, the calibration at $z=0$ is performed with the published estimate of local merger rate with the power-law mass distribution333At the time of writing the paper, Advanced LIGO and Virgo collaboration published a new paper on the inferred BBH population property using their first and second observing runs and the local merger rate estimate has been updated. However, we stick to the local merger rate, previously estimated in Abbott et al. (2017b), to demonstrate the pipeline performance. The normalization with the new merger rate estimate will be considered in future work., $103~{}\mathrm{Gpc}^{-3}\mathrm{yr}^{-1}$ Abbott et al. (2017b) such that

[TABLE]

Although we have assumed this rate is known exactly, we note that the true rate is uncertain – the current 90% credible bounds on the BBH merger rate span [9.7, 101]$~{}\mathrm{Gpc}^{-3}\mathrm{yr}^{-1}$. This uncertainty will correspondingly introduce a systematic uncertainty in our prediction of the remnant spectrum $\Omega^{\mathrm{rem}}_{\textsc{gw}}(f)$, whose effect should be considered in future work.

III.3 Total background model

The actual background spectrum that one would observe is the sum of the contributions from these two BH populations, that is,

[TABLE]

where the superscripts represent each of the isolated BH and BBH merger remnant population given by Eq. (17) and Eq. (24), respectively. For the comparison between these components, Fig. 2 shows background spectra computed from each BH population model, plotted with power-law integrated curves Thrane (2013) of Advanced LIGO at different phases and the CBC background approximated as a power-law spectrum Abbott et al. (2017e).

The solid curves are the energy density spectra contributed from the isolated BH population for different scalar masses under an assumption of a uniform distribution $p(\chi)$ over $0\text{\,}1\text{\,}$. The spin distribution skewed towards high $\chi$ would simply scale the overall spectrum. See Fig. 2 of Brito (2017) for its dependence on different spin distributions. The dotted curves, meanwhile, show the contribution from BBH merger remnants. In general, the abundance of isolated BHs exceeds the merger rate by four orders of magnitude. The isolated BH channel therefore dominates the SGWB for scalar masses $m_{b}\geq 10^{-12.5}\mathrm{eV}$. However, if the $m_{b}\sim 10^{-13}\mathrm{eV}$, the SGWB from isolated BHs is significantly suppressed due to the lack of BHs heavier than $50\,M_{\odot}$. In other words, although scalar fields with $m_{\mathrm{b}}\sim${10}^{-13}\text{,}\mathrm{e}\mathrm{V}$$ is sensitively coupled to these heavy BHs (see Fig. 1 of Brito et al. (2017)), such massive BHs are only rarely produced by CCSNe Fryer et al. (2012). In contrast, heavy BHs can readily be produced by the merger remnant channel, and so remnant BHs dominate the SGWB when $m_{b}\sim 10^{-13}\mathrm{eV}$.

It is worth noting that the SGWB may dominate the projected CBC background, which is approximated as

[TABLE]

as inferred from Abbott et al. (2017e). We will discuss the capability to distinguish between these two signal models based on the Bayesian framework in Section V.4.

IV Search Method

This section provides the overview of a statistical methodology that can be used to claim a detection or to place constraints related to the model of the superradiant instability. Here, we assume the GW background to be; (a) isotropic, (b) unpolarized, (c) stationary and (d) Gaussian.

IV.1 Setup

The SGWB is analyzed by taking the cross-correlation between outputs from a pair of detectors. In this work we will focus on the case of two detectors for simplicity, but the formalism can be extended to handle a larger network of detectors Allen and Romano (1999); Romano and Cornish (2017). Following the notation of Callister et al. (2017), we define a cross-correlation estimator for each frequency bin Aasi et al. (2015b); Romano and Cornish (2017) as,

[TABLE]

Here, $\tilde{s}_{i}(f)$ is the Fourier transform of time series output of the $i$-th detector and $T$ is the total observation time. This is normalized such that

[TABLE]

where $\gamma(f)$ is the overlap reduction function to encode the geometry and separation between a pair of detectors Flanagan (1993). In the low signal-to-noise limit, the variance is approximately given by

[TABLE]

where $\Delta f$ is frequency resolution and $P_{i}(f)$ is the power spectral density (PSD) of the $i$-th detector.

Then, the unbiased broadband estimator and its variance can be constructed from $\hat{C}(f)$ and $\sigma^{2}(f)$ as follows

[TABLE]

and

[TABLE]

where $w(f)$ is a linear filter formed by a given background model $\Omega_{\mathcal{A}}(f)$, so that

[TABLE]

Once the broadband estimator is obtained, the SNR can be computed in a straightforward manner.

[TABLE]

IV.2 Bayesian inference

We now develop a Bayesian formalism to perform parameter estimation and model selection for the superradiant instability model, following Mandic et al. (2012). In practice, the GW signal can be described by a set of unknown parameters $\bm{\theta}=\{\theta_{1},\theta_{2}...,\theta_{N}\}$ and one needs to evaluate their posterior probability given new data. Then, Bayes’ theorem states that

[TABLE]

where $\{\hat{C}\}$ is the cross-correlation estimator over frequency band computed from observed data and $\bm{\theta}_{\mathcal{A}}$ is a set of parameters characteristic of GW signals in the signal hypothesis $\mathcal{A}$. $p(\bm{\theta}_{\mathcal{A}}|\{\hat{C}\},\mathcal{A})$ is the posterior probability on the multi-dimensional space, $L(\{\hat{C}\}|\bm{\theta}_{\mathcal{A}},\mathcal{A})$ is the likelihood, $\pi(\bm{\theta}_{\mathcal{A}}|\mathcal{A})$ is the prior probability and $Z(\{\hat{C}\}|\mathcal{A})$ is the evidence. To construct a posterior in an efficient manner, we apply the PyMultiNest package Buchner, J. et al. (2014) to our search pipeline. PyMultiNest is a python interface to the nested sampling package MultiNest Feroz and Hobson (2008); Feroz et al. (2009, 2013), which produces a set of samples drawn from an estimated posterior.

IV.2.1 Likelihood

For one realization of $\hat{C}(f)$, the joint likelihood over frequency bins is given by Mandic et al. (2012)

[TABLE]

$L(\hat{C}(f)|\bm{\theta}_{\mathcal{A}},\mathcal{A})$ is the likelihood within a single frequency bin. We assume a Gaussian likelihood, such that

[TABLE]

Here, $\Omega_{\mathcal{A}}(f;\bm{\theta}_{\mathcal{A}})$ is a model energy-density spectrum for a given set of parameters $\bm{\theta}_{\mathcal{A}}$.

IV.2.2 Posterior

We choose a uniform prior on every parameter (e.g. scalar mass $m_{b}$ and spin limits $\chi_{\mathrm{ul}/\mathrm{ll}}$):

[TABLE]

Our posteriors are therefore proportional to the likelihood:

[TABLE]

We will often consider the marginalized posterior for a particular parameter, $\theta_{1}$, which is defined as

[TABLE]

The marginalized posterior will allow us to define credible intervals for the parameter $\theta_{1}$ in a standard way.

IV.2.3 Model selection

In order to perform model selection between different models of the SGWB, we need to compute the Bayesian evidence for each hypothesis. Let $\{\hat{C}\}$ be the cross-correlation estimator obtained from the data and suppose that one assesses which model, $\mathcal{A}$ or $\mathcal{B}$, is better supported by the data. It is straightforward to compute the odds ratio $\mathcal{O}^{\mathcal{A}}_{\mathcal{B}}$, which is defined as

[TABLE]

where $Z(\{\hat{C}\}|\mathcal{A})$ and $Z(\{\hat{C}\}|\mathcal{B})$ are the evidences for each model. The evidence $Z$ is calculated through

[TABLE]

This expression can also be interpreted as the fully-marginalized likelihood. In the case where no signal is present (the null hypothesis $H_{0}$), the evidence is obtained by fixing $\Omega_{\mathcal{A}}(f;\theta_{\mathcal{A}})$ to zero. Also $Z(\{\hat{C}\}|\mathcal{A})/Z(\{\hat{C}\}|\mathcal{B})$ in Eq. (46) is the so-called Bayes factor, while $\pi(\mathcal{A})/\pi(\mathcal{B})$ is a priori probability ratio for the two models, which we will set as unity. Since in our framework the odds ratio Eq. (46) is effectively equivalent to the Bayes factor, hereafter we will evaluate statistical significance in terms of the Bayes factor and follow the convention that a log Bayes factor of $\approx 8$ indicates a favor for a model over the other with great confidence Jeffreys (1961).

V Results

In this section, we first describe the parameterization of our signal model. Section V.2 shows a scheme for injecting simulated signals into Gaussian noise and recovering those signals. We then discuss the implications of performed analysis tests, including the construction of a “sensitivity window” of scalar masses to which Advanced LIGO might be sensitive. In Section V.4, we examine prospects for successfully discriminating between a CBC background and a SGWB due to supperradiant instability.

V.1 Parameters for the background model

In this work, the signal model is parameterized by the scalar mass $m_{b}$ and either $\chi_{\mathrm{ll}}$ or $\chi_{\mathrm{ul}}$ (refer to Sections III.2.1 for the details). We note that this model contains several sources of systematic uncertainty, such as the BH formation and merger rates as well as the underlying SFR. Nevertheless, we will assume a particular SFR (Eq. (21)) and binary formation rate (Eq. (30)). Our results are also conditional on our specific parameterization of the $p(\chi)$. A different choice for the BH spin distribution will generically yield different constraints on $m_{b}$. Also, posteriors will be constructed over the two-dimensional parameter space of $m_{b}$ and either $\chi_{\mathrm{ll}}$ or $\chi_{\mathrm{ul}}$.

V.2 Injection scheme and signal recovery

Here we describe the overview of the injection scheme adopted in this work. We assume only the LIGO detector pair, that is, the Hanford and Livingston sites. Injections are performed in the frequency domain. Given a PSD of the LIGO detectors, one constructs the variance of a cross-correlation estimator given by Eq. (35). A cross-correlation spectrum consistent with Gaussian noise is then constructed from this variance and the SGWB predicted from superradiant instability is added to the Gaussian noise. Therefore, the simulated cross-correlation spectrum is defined as

[TABLE]

where \gamma_{\mathrm{HL}}(f)\ is the overlap reduction function for the LIGO baseline Flanagan (1993), $\Omega_{\mathrm{inj}}(f;\bm{\theta}_{\mathrm{inj}})$ is the injected background constructed from the parameters $\bm{\theta}_{\mathrm{inj}}$, and $\hat{n}$ is a random variable drawn from a Gaussian distribution with zero mean and unit variance. From this simulated cross-correlation spectrum, our pipeline infers the parameters $\bm{\theta}_{\mathrm{inj}}$ by computing the likelihood and constructing a Bayesian posterior. We adopt a uniform prior for $m_{b}$ across the range ${10}^{-13}\text{\,}\mathrm{e}\mathrm{V}{10}^{-12}\text{\,}\mathrm{e}\mathrm{V}$ and a uniform prior for $\chi_{\mathrm{ul}}$ or $\chi_{\mathrm{ll}}$ between 0 and 1.

Fig. 3 shows an example cross-correlation spectrum $\hat{C}_{\mathrm{sim}}(f)$ as a function of frequency using two years of observation with Advanced LIGO’s design sensitivity. The red dashed line gives the injected signal. The injection has SNR $\sim$ 32 and corresponds to the superradiant instability model with $m_{\mathrm{b}}=$7\text{\times}{10}^{-13}\text{,}\mathrm{e}\mathrm{V}$$ and $\chi_{\mathrm{ll}}=0.5$. The magenta region represents the model evaluated with 100 sets of the parameters randomly drawn from the posterior. Fig. 4 shows the resulting posterior on $m_{b}$ and $\chi_{\mathrm{ll}}$ derived from this simulated observation. These figures both demonstrate that the pipeline can recover a loud injection with great confidence. A statistical discussion for a consistency test of the parameter estimation is described in Appendix A.

V.3 Sensitivity window

To study the scalar mass space we can probe through this model, which we refer to as the “sensitivity window”, we make a number of injections and compute their Bayes factors between the signal and noise hypotheses. The signal injection is performed in the same way as Eq. (48). For this test, we conservatively adopt the $\chi_{\mathrm{ul}}$ parameterization and the injected $\chi_{\mathrm{ul}}$ value is fixed as $0.8$. The $m_{b}$ value for each of the 500 injections is uniformly drawn from the range $0.0\text{\times}{10}^{-13}\text{\,}\mathrm{e}\mathrm{V}0.0\text{\times}{10}^{-12}\text{\,}\mathrm{e}\mathrm{V}$. In Fig. 5, we recognize the injections above a log Bayes factor of 8 as detectable. Therefore, the $m_{b}$ range in which a log Bayes factor is larger than 8 (presented as the red line in Fig. 5) can be interpreted as the sensitivity window, which is around $1.8\text{\times}{10}^{-13}\text{\,}\mathrm{e}\mathrm{V}7.5\text{\times}{10}^{-13}\text{\,}\mathrm{e}\mathrm{V}$ in Fig. 5. Note that the actual sensitivity window depends on spin parameter $\chi_{\mathrm{ul}}$. Since larger $\chi_{\mathrm{ul}}$ simply scales the overall spectrum by a constant factor, we can have a wider detectable range of scalar masses for a larger $\chi_{\mathrm{ul}}$.

V.4 Distinguishing superradiant instability and CBC backgrounds

As seen in Fig. 2, the SGWB signal from the superradiant instability model dominates over the projected CBC background for some choices of scalar mass and BH spins. Thus, the next question to address is “Given a detected signal, can we distinguish these two models from one another?” Here, we consider the case where both the superradiant instability and the projected CBC background are present, so that

[TABLE]

where $\Omega_{\mathrm{inj}}^{\mathrm{SI}}(f;\bm{\theta})$ is the background due to superradiant instabilities under the $\chi_{\mathrm{ul}}$ parametrization given by Eq. (31), and $\Omega_{\mathrm{inj}}^{\mathrm{CBC}}(f)$ is the fixed CBC background approximated as the power law spectrum shown in Eq. (32). Note that, after three years of observation with Advanced LIGO’s design sensitivity, $\Omega_{\mathrm{inj}}^{\mathrm{CBC}}(f)$ alone is detectable with a log Bayes factor of 8.8 between a CBC-only and a noise model, which corresponds to an SNR of 5.5.

Given a simulated measurement of the combined background shown in Eq. (49), we recover our measurements with different two models: a CBC-only model and a joint superradiant instabiltity and CBC model (SI$+$CBC). The CBC background model for recovery is parameterized by a power-law index and an amplitude at the reference frequency $25\mathrm{Hz}$, such that

[TABLE]

We refer to Eq. (50) as a CBC model, even though this could be generally called “power-law spectrum model”, because the degrees of freedom in this model can be used to recover the parameters of the CBC injection. Also, priors for $\Omega_{0}$ and $\alpha$ are taken as a uniform distribution in the range of ${10}^{-10}\text{\,}{10}^{-17}\text{\,}$ and $-5\text{\,}5\text{\,}$ respectively, while the priors for $m_{b}$ and $\chi_{\mathrm{ul}}$ are identical to those described in Section V.2. The parameters considered when evaluating the evidence for each background model are listed in Table 1.

We compute a log Bayes factor between these two hypotheses. The computation is repeated, changing injected $(m_{b},\chi_{\mathrm{ul}})$ values until we explore a grid over the entire prior space. Thus, we can probe the parameter space on which the superradiant instability signal can be discerned from the expected CBC background with statistical significance. Fig. 6 is a gray-scale map of a log Bayes factor, effectively equivalent to Eq. (46), with two contours of $\ln(\mathrm{BF})=8$ (magenta) and [math] (cyan). Since the injected CBC background is loud enough to be detected in this simulation, Fig. 6 implies that inside the magenta contour one can expect the superradiant instability background to be distinguished from the CBC’s.

VI Application to the first observing run of Advanced LIGO

We now apply our method to the cross-correlation spectrum measured in Advanced LIGO’s first observing run Abbott et al. (2017f) (the data products given in that paper are publicly available at LIGO-G (1602446)). No statistically significant signal is detected, so we attempt to exclude some of the two-dimensional space of $m_{b}$ and $\chi_{\mathrm{ul},\mathrm{ll}}$ through the search. Figs. 7 and 8 show the posteriors obtained under the $\chi_{\mathrm{ul}}$ and $\chi_{\mathrm{ll}}$ parametrizations, respectively. Fig. 7 indicates that, for the $\chi_{\mathrm{ul}}$ case, the posterior is almost uniformly distributed. Thus, no significant constraint can currently be placed on $m_{b}$ when $\chi_{\mathrm{ul}}$ is allowed to vary. This can be also verified by the fact that, according to Fig. 2, almost none of the possible spectra with the $\chi_{\mathrm{ul}}$ parameterization can reach the O1 power-law integrated curve. On the other hand, if we fix $\chi_{\mathrm{ul}}$ and allow the lower bound $\chi_{\mathrm{ll}}$ to vary, then Fig. 8 suggests that the mass range 2.0\text{\times}{10}^{-13}\text{\,}\mathrm{e}\mathrm{V}$\leq m_{\mathrm{b}}\leq$3.8\text{\times}{10}^{-13}\text{\,}\mathrm{e}\mathrm{V} is excluded. This can be understood by Fig. 2 as well, which shows the largest SGWB amplitude when $m_{b}\sim$e-12.5\text{,}\mathrm{eV}$\approx$3.2\text{\times}{10}^{-13}\text{,}\mathrm{eV}$$. As a forecast for the future, Appendix B shows the constraints that will be possible once Advanced LIGO reaches its design sensitivity.

VII Conclusion

This paper presents a first search for signs of supperradiant instability in the SGWB. We describe several performance tests of our search pipeline in Section IV. First, a detectable window in scalar mass is estimated from a Bayes factor of injections as shown in Fig. 5. Second, we study the capability to distinguish the model presented here from the fiducial CBC background model. The gray-scale map of log Bayes factors suggests that it is separable from the CBC model with statistical significance in some parameter space (see Fig. 6).

Finally, we present results obtained by analyzing data from Advanced LIGO’s first observing run. No signal is detected with either of our parameterizations of the BH spin distribution, and so we present constraints on possible boson masses and BH spin bounds. Using the $\chi_{\mathrm{ul}}$ parametrization, we cannot place any meaningful constraints on the scalar mass, while the $\chi_{\mathrm{ll}}$ case rules out the scalar mass range 2.0\text{\times}{10}^{-13}\text{\,}\mathrm{e}\mathrm{V}$\leq m_{\mathrm{b}}\leq$3.8\text{\times}{10}^{-13}\text{\,}\mathrm{e}\mathrm{V} with 95% percent credibility. Let us note that this constraint is still subject to our choice of the BH spin distribution as well as specific astrophysical models we adopt in this work.

Future work will generalize the model to take into account effects of astrophysical uncertainties such as BH formation rate and the initial spin distribution for BBH remnant BHs. Also, the updated local merger rate Abbott et al. (2018b) should be revisited. As the sensitivity increases with future observing runs of Advanced LIGO and Advanced Virgo, this framework can be used to place stronger constraints on the existence of ultralight bosons.

Acknowledgements.

We thank Richard Brito and Irina Dvorkin for their fruitful discussion and comments. L.T is supported by JSPS KAKENHI Grant Number JP18J21709. 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-0757058. A.M and P.M are supported by NSF grant PHY1505870. Parts of this research were conducted by the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav), through project number CE170100004. This paper carries LIGO Document Number LIGO-P1800232.

Appendix A Consistency test of the parameter estimation

As a statistical consistency check of the parameter estimation performed by our pipeline using PyMultiNest, we construct a probability-probability plot (p-p plot) by performing 500 injections into simulated Gaussian noise drawn from the observed $\sigma(f)$ spectrum from Advanced LIGO’s O1 run. For each injection, we first construct the marginalized posterior for $m_{b}$. We then compute the percentile (or p-value), which is defined by

[TABLE]

where $i$ denotes the label of a given injection, $p(m_{b})$ the marginalized posterior for $m_{b}$, $m^{*}_{b}$ the injected value and $m_{0}$ the lower limit of the $m_{b}$ prior space. We then plot the fraction of injections with $p$-values less then a given threshold $p_{i}$, as a function of the threshold. For correctly constructed posteriors we expect that for a given threshold $p_{i}$, the fraction of injections with p-values smaller than the threshold is $p_{i}$. If this is the case, then the injections will form a straight line with unit slope on the pp-plot. Fig. 9 shows the pp-plot we obtain with our injection campaign. The gray shadow is the $95\%$ credible error region predicted by statistical fluctuations from the ideal diagonal. Since the orange curve obtained by data lies within the error region, the pipeline consistently reproduces the parameters regardless of the loudness of those injections.

Appendix B Demonstration of possible constraints on scalar mass at the design sensitivity

We run the pipeline with synthesized noise data to demonstrate constraints on scalar mass expected from null result by three years of observation with the Advanced LIGO’s design sensitivity. Following the definition Eq. (33) and Eq. (35), the cross-correlation estimator $\hat{C}(f)$ and its variance $\sigma^{2}(f)$ for individual frequency bins are computed.

Figs. 10 and 11 each shows a two-dimensional posterior with different $p(\chi)$ parameterizations. When $\chi_{\mathrm{ll}}$ is left free, the posterior excludes all scalar masses above $m_{b}\geq$1.4\text{\times}{10}^{-13}\text{,}\mathrm{eV}, as shown in Fig. [10](#A2.F10). In Fig. [11](#A2.F11), on the other hand, more stringent constraints are placed in lighter scalar mass and the constraints become looser in the heavier boson mass regime. This is due to the strong dependence of the GW background amplitude on $\chi_{\mathrm{ul}}$. More interestingly, the projected design sensitivity would place quite different constraints depending on the $p(\chi)$ parameterization. While $\chi_{\mathrm{ll}}$ parameterization excludes upper part of prior space, $m_{b}\geq$1.4\text{\times}{10}^{-13}\text{\,}\mathrm{eV}, the $\chi_{\mathrm{ul}}$ parameterization could exclude the lower edge.