Constraining black-hole spins with gravitational wave observations

Vaibhav Tiwari, Stephen Fairhurst, Mark Hannam

I Introduction

Gravitational waves (GW) emitted by merging black holes are identified in the LIGO and Virgo data through the use of search analysis pipelines, which use the known waveform morphology to identify weak signals in the data (Abbott et al. 2016a; Abbott et al. 2017a; Abbott et al. 2017b; Abbott et al. 2017c). The observations are followed by parameter-estimation analyses that extract posterior probability distributions for the parameters of the binary — the masses and spins of the component black holes as well as the distance, sky location and orientation of the binary (Cutler & Flanagan 1994; Abbott et al. 2016b; Veitch et al. 2015). While some parameters are extracted with good precision, others cannot be accurately measured, and several sets of parameters are strongly correlated, for example the distance with binary orientation, and mass ratio with black hole spins. Nonetheless, the observed parameters from several observations can be combined to obtain the underlying astrophysical distributions of black-hole masses and spins. In this paper we use publicly available information of the measurements from the first six GW signals observed from merging black holes (GW150914, LVT151012, GW151226 and GW170104, GW170608 and GW170814)(Abbott et al. 2016a; Abbott et al. 2017a; Abbott et al. 2017b; Abbott et al. 2017c) to draw inferences about the underlying spin distribution of black holes We generate parameter distribution using confidence intervals reported in observation papers – please see section on method used in this paper..

The effective spin used in LIGO-Virgo analyses is related to the individual spins of the two black holes in the binary by (Damour 2001; Ajith et al. 2011)

where m1m_{1} and m2m_{2} are the component masses of the binary and χ1\chi_{1} and χ2\chi_{2} are the components of the (dimensionless) spins aligned with the angular momentum defined as χ=cS⋅L^/(Gm2)≤1\chi=c\mathbf{S\cdot\hat{L}}/(Gm^{2})\leq 1.

We are interested in inferring the spin distribution of the merging black holes binaries in the Universe by comparing the observed distribution of effective spins with those predicted by astrophysical models. At present, we have only a limited number of observations and therefore we focus on a discrete set of models. We follow (Farr et al. 2017) and introduce six possible spin distributions for comparison. We use three distributions for spin magnitude:

I.2 Two Primary Effects

Before we proceed to model selection we must consider two primary effects black hole spins have on the waveform.

Ideally, a combined analysis of masses and spins will naturally account for these effects. A flexible non-parametric prior, that maximizes the overall probability of observing all the gravitational wave signals, can be used to obtain the parameter distributions. Such an analysis will require hundreds of events (See also Wysocki et al. 2018 which constructs a phenomenological distribution with limited number of gravitational wave observations).

At present, we have only a limited number of observations and therefore we account for these effects by imposing an astrophysically motivated mass distribution on the universe: p(m1)∝m1−2.3p(m_{1})\propto m_{1}^{-2.3} with m2m_{2} uniformly distributed between 5 M⊙ and m1m_{1}. The choice is based on astrophysical observations that support the stellar initial mass function to follow a power-law distribution. Moreover, the power-law model provides a binary mass-ratio distribution supported by the population synthesis models Dominik et al. 2012; Rodriguez et al. 2016. Furthermore, independent of the assumed spin distribution, the GW measurement also support power-law model Abbott et al. 2017a. In summary, power-law is among simple models that are supported by the data.

II Method

Using the observed measurements of the effective spin from the six BBH mergers considered here, we use Bayesian model selection to calculate the odds ratio between the different models. While model selection is quite standard, care must be taken to ensure that the selection effects and mass priors are correctly incorporated; see also Loredo 2004; Mandel et al. 2016. We are interested in calculating

where {d}\{\bm{d}\} denotes the set of observations, p(λ)p(\lambda) is the prior on the model λ\lambda and p({d})p(\{\bm{d}\}) is formally given as the integral over λ\lambda of the numerator:

Since the model, λ\lambda, gives a distribution for the parameters of the signal, θ\bm{\theta}, we can express the probability of obtaining a given data set d\bm{d} corresponding to a single observation as

where the distribution of d\bm{d} given θ\bm{\theta} is calculated from the Gaussian likelihood as

where dXd_{X} denotes the data in detector ‘X’ and and hX(θ)h_{X}(\bm{\theta}) is the gravitational waveform expected in detector ‘X’ from a binary with parameters θ\bm{\theta}. The product is over detectors in the network, ⟨a∣b⟩\langle a|b\rangle is the noise weighted inner product, defined in the frequency domain as,

and S(f)S(f) is the power spectrum of the detector noise (Cutler & Flanagan 1994).

We must also take into account the fact that there is a separate threshold on the search. This arises in the normalization of the probability density for d\bm{d} above. When there is no threshold, the probability distribution in Equation 5 is correctly normalized. However, when we impose a threshold, we must take into account that not all sets of parameters θ\bm{\theta} are equally likely to lead to the identification of a gravitational wave signal. Thus, to normalize the probability, we must integrate the probability over all realizations of the data, d\bm{d}, which produces an event above the threshold ρ⋆\rho_{\star}.

is the total volume. Sensitive volume is a primary ingredient in accounting for the selection effects and can be estimated using semi-analytical (Abbott et al. 2016c) or numerical methods (Tiwari 2018).

We can now express the distribution for the data d\bm{d}, corresponding to a single observation, given the model λ\lambda as

The expression generalizes in a straightforward manner to a population of observed events, as we assume that the parameters of the signals are independent, such that

Finally, we can use Equation 3 to obtain the probability of a model λ\lambda given the set of observations {d}\{\bm{d}\} as:

This can then be used in a straightforward manner to perform model selection between two models λ1\lambda_{1} and λ2\lambda_{2} as

The three terms in the odds ratio are easily understood. The final term is simply the ratio of the priors for the two models. In this paper, when comparing models, we take an equal prior probability for the models so this term is equal to unity. The middle term is the probability for observing the data di\bm{d}_{i} given the model λ\lambda and the first term arises as a normalization due to the threshold in the identification of signals in the data. We note that the overall prior on the data p({d})p(\{\bm{d}\}) cancels as it appears in the same way for both the models.

When performing parameter estimation, we obtain a set of posterior samples θj\theta^{j} that describe the posterior distribution for θ\bm{\theta}. Explicitly, we can approximate an integral over the parameter space as

Thus, the integral in Equation 15 can be well approximated by a sum over the (appropriately weighted) posterior samples

The first term gives the ratio of the sensitive volumes for the two models, and favours the model with the lower sensitive volume. The second term sums over the re-weighted posterior samples, where the re-weighting factor is simply the ratio of the desired prior to the one used in the parameter estimation. The final term is the ratio of the priors of the two models. Let us now look in detail at the impact of the three factors when calculating the odds ratios. As discussed previously, we will always assume an equal prior between models, so the final term is unity.

In order to estimate the sensitive volumes of the different population models, defined in Equation 10, we perform Monte Carlo integration. To do so, we sample from the astrophysically expected distributions of parameters and determine the fraction of sources that would be observed. The six populations used in the analysis follow the same mass distribution but different spin distributions. Random samples of the masses and spins are drawn from the population and are assigned randomly chosen orientation and sky location. Samples are distributed in redshift as determined by standard cosmology. The expected signal-to-noise ratio (SNR) of the signals produced by these binaries at the detectors are calculated and signals that cross a certain SNR threshold are labeled as recovered. Since all of the events apart from GW170814 were observed by only the LIGO detectors, for simplicity we estimate the sensitive volume for the LIGO Hanford – LIGO Livingston detector network. We choose a fixed power spectrum for the detector noise operating close to the sensitivity of the LIGO detectors during the first observing run. While the sensitivity of the detectors varies over the runs, and between the first and second observing runs, the ratio of sensitive volumes for the different population models is relatively insensitive to changes in the detector sensitivity. To optimize the calculation, we estimate the SNR for face-on signals on a fiducial grid of binary masses and spins. The expected SNR of a binary with arbitrary masses, spins, location and orientation is calculated by linear interpolation in the mass and spin space and incorporating the loss in SNR due to the arbitrary orientation Schutz 2011. To test the efficacy of the procedure, we compare our result with the results reported in reference (Tiwari 2018). Our results are within 10% of the reported values.

The Monte Carlo equivalent of Equation 10 is given by

Next, let us consider the re-weighting of the posterior samples. In particular, the prior, π(θ)\pi(\bm{\theta}) is usually taken to be flat in m1m_{1} and m2m_{2}, subject to the condition that m1>m2m_{1}>m_{2} and flat in the z-components of the spin (Abbott et al. 2016b; Abbott et al. 2016a; Abbott et al. 2017a). In computing the probabilities for the various models under consideration we must vary the spin prior to match one of the six distributions under consideration. In addition, we would like to use a different mass prior, which is astrophysically motivated (Fishbach & Holz 2017; Talbot & Thrane 2018). In particular, we select p(m1)∝m1−αp(m_{1})\propto m_{1}^{-\alpha} and p(m2)p(m_{2}) uniform in m2m_{2} between 5 M⊙ and m1m_{1} (Abbott et al. 2016c; Abbott et al. 2016a). In Figure 1 we show the prior distribution for the mass ratio, given the flat prior π(θ)\pi(\bm{\theta}). In addition, we show the distribution that is obtained with the power law prior, with α=2.3\alpha=2.3, the value used to obtain the results, as well as α=0.9\alpha=0.9 and α=3.3\alpha=3.3. These values, with mean at α=2.3\alpha=2.3, cover one-sigma confidence interval of the possible values of α\alpha that are consistent with the observations (Abbott et al. 2017d).

The mass-ratio distribution based on an astrophysical model is significantly different from the one obtained with a flat prior on the component masses. In particular, close to equal mass ratio systems are significantly more likely while binaries with mass ratios greater than 5:1 are down-weighted by a factor between 1.5 and 30, depending upon the value of α\alpha. The higher the value of the power-law slope α\alpha, the more the distribution is skewed towards equal masses. Since there is a degeneracy between mass-ratio and aligned spin, a preference for close to equal mass binaries will provide tighter posteriors on the spin. So, applying an astrophysically motivated mass prior leads to a preference for lower spin values. This has a significant impact in down-weighting the high-spin distributions when summing over the re-weighted posterior distributions in Equation 18.

Finally, we discuss the generation of the posterior samples. The full LIGO-Virgo analysis produces thousands of posterior samples which are used to obtain the parameter distributions presented in the results papers. In principle, those could be used directly in the calculation of the odds ratio in Equation 18. However, at present only samples corresponding to LIGO’s first scientific run are publicly available (Abbott et al. 2018) so we instead generate distributions that mimic the full parameter estimation results, based upon the parameter values and uncertainties presented in the results papers. As we are interested in only the masses and spins, we do not perform a full parameter estimation analysis, but rather use the fact that the measurement of masses and spins is largely independent of the sky location and binary orientation (Singer & Price 2016).

Table 1 compares the credible intervals for the announced GW observations (Abbott et al. 2016a; Abbott et al. 2017a; Abbott et al. 2017b; Abbott et al. 2017c) with the credible intervals of the samples obtained using the method described above. The intervals obtained for the masses and effective spins are comparable. We have made several approximations which we would expect to give differences at the observed levels. In addition, the reported gravitational wave results make use the average results from both spin-aligned and full-spin precessing models, while in the posteriors that we generate we consider only aligned spins. The mass-ratios of the approximated posteriors are nearly equal to the reported values. We have verified that our results are insensitive to shifts in the posterior distributions at the levels reported in Table 1.

III Results and Discussion

Impact of selection effects are shown in Figure 3 that plots the ratio of the sensitive volume of the spin models and the sensitive volume of the low isotropic model. Binaries with a higher spin magnitude can be observed at a greater distance, so we expect that the low isotropic model, which leads to the population with the smallest spin magnitudes, to have the lowest sensitive volume. Thus, all other things being equal, for each event observed the model with the lower sensitive volume is preferred. This has a significant impact for the aligned spin models, but the volume ratio between models with isotropic spin distributions is close to unity. Thus, this doesn’t have a significant impact for a small number of events, but does give a factor of 5 contribution to the odds ratio with 50 events.

The results of the analysis are shown in Figure 4. This shows the odds ratio between each of the six models discussed above and the low isotropic model. The low isotropic distribution is preferred. All models with aligned spins are disfavoured at greater than 22000:1 or, equivalently, 4σ4\sigma. (The flat- and high-aligned distributions are disfavoured at >>5σ\sigma.) This provides strong evidence against aligned spins. The improvement arise from three different factors. The inclusion of two additional events (GW170608 and GW170814) increases the odds ratio by a factor of six, while an accurate treatment of the selection effects and the use of an approximately good mass ratio distribution to handle the mass ratio–spin degeneracy increase the odds by a factor of sixty. Incorporation of selection effects increases the odds ratio by a factor of four while accounting for correlations between mass ratio and spin increases the odds ratio by a factor of sixteen. Thus, these three effects combined explain the factor of around four hundred improvement in our ability to exclude aligned spin models in favour of isotropic distribution of spins.

However, as noted in this paper, systematic and selection effects will affect the inference made by a mixture model as it will likely introduce a bias towards LA models e.g. the odds-ratio for GW151226 reduces for the LA model on using re-weighted posterior samples instead of the standard posterior samples and odds-ratio for GW170608 reduces to close to unity.

The result is only mildly dependent upon the value of the slope, α\alpha. As we increase the value of α\alpha, we favour lower mass black hole binaries in the population while lower values of α\alpha have a larger fraction of high mass binaries. The impact of spin on visible volume is more significant for higher masses so, consequently, small values of α\alpha lead to larger difference in sensitive volumes between different spin distributions. However, for larger α\alpha, the mass ratio and aligned spin are more tightly constrained. This increases the importance of re-weighting the posterior samples. Thus, overall, changing the value of α\alpha has a limited impact on the results. We also note that choice of the minimum and the maximum masses in the power-law model has only mild effect on the results. Table 2 lists odds-ratios of HI and LA spin models in reference to LI model for different values of α\alpha and the maximum mass of the primary component of the binary.

We have shown that the first six observations of black hole binary mergers can be used to place a limit on the magnitude of black hole spins based on gravitational wave observations. The data show strong evidence for isotropic, rather than aligned, spins. Furthermore, there is emerging evidence that small low spin magnitudes are preferred to high spin magnitudes (see Figure 8 of (Farr et al. 2018) as well as of (Wysocki et al. 2018)). This contrasts the nominal spin magnitudes inferred from x-ray binary observations, which is more consistent with high spins.

We emphasise that these conclusions depend on our choices of possible spin distributions. If, for example, our low-spin distribution was restricted to much lower values of χeff\chi_{\rm eff}, then aligned-spin configurations would not be so strongly disfavoured. However, this would only strengthen our main conclusion, which is the preference for low spin magnitudes.

Distributions of black hole spins will be further refined through future gravitational wave observations. In the third advanced LIGO-Virgo observing run, there is an expectation of observing tens of black hole binary mergers. To get a sense of what we might expect, we simulated 30 observations from the low isotropic spin distribution and combined them with the six observations discussed above. With this set of observations, we would be able to exclude a population of black hole binaries where 20% have aligned spins and 80% with isotropic spins with a confidence of at least 4σ4\sigma. Furthermore, we would obtain an odds ratio in favour of low isotropic spins over high isotropic spins of around 130,000:1 (4.5σ4.5\sigma) and flat isotropic of around 300:1 (3.0σ3.0\sigma). Finally, we note that we have not made use of the precessing component of the spin. A clear observation of precession will give irrefutable evidence of spin misalignment, while observations of χp\chi_{p} consistent with zero will provide further evidence against high spin magnitudes.

Acknowledgments

We would thank the following for interesting and useful discussions: Will Farr, Frank Ohme, Richard O’Shaughnessy, Simon Stevenson, Eric Thrane and Vivien Raymond. The authors were funded by the Science and Technology Facilities Council (STFC) grants ST/L000962/1 and ST/N005430/1. MDH and SF was supported by the European Research Council Consolidator Grant 647839.

References