THE MASS DISTRIBUTION OF STELLAR-MASS BLACK HOLES

The most massive stars probably end their lives with a supernova explosion or a quiet core collapse, becoming stellar-mass black holes. The mass distribution of such black holes can provide important clues to the end stages of evolution of these stars. In addition, the mass distribution of stellar-mass black holes is an important input in calculations of rates of gravitational wave emission events from coalescing neutron star–black hole and black hole–black hole binaries in the LIGO gravitational wave observatory (Abadie et al. 2010).

Observations of X-ray binaries in both the optical and X-ray bands can provide a measurement of the mass of the compact object in these systems. The current sample of stellar-mass black holes with dynamically measured masses includes 15 systems with low-mass, Roche lobe overflowing donors and 5 wind-fed systems with high-mass donors. Hence, sophisticated statistical analyses of the black hole mass distribution in these systems are possible.

The first study of the mass distribution of stellar-mass black holes, in Bailyn et al. (1998), examined a sample of seven low-mass X-ray binaries thought to contain a black hole, concluding in a Bayesian analysis that the mass function was strongly peaked around seven solar masses.⁵ Bailyn et al. (1998) found evidence of a "gap" between the least massive black hole and a "safe" upper limit for neutron star masses of 3 M_☉ (e.g., Kalogera & Baym 1996). Such a gap is puzzling in light of theoretical studies that predict a continuous distribution of compact object supernova remnant masses with a smooth transition from neutron stars to black holes (Fryer & Kalogera 2001). (We note that Fryer & Kalogera (2001) considered binary evolution effects only heuristically and put forward some possible explanations for the gap from Bailyn et al. (1998) both in the context of selection effects or in connection to the energetics of supernova explosions.)

Toward the end of our analysis work, we became aware of a more recent study (Ozel et al. 2010), also in a Bayesian framework, analyzing the low-mass X-ray binary sample. Our results are largely consistent with those obtained by Ozel et al. (2010), who examined 16 low-mass X-ray binary systems containing black holes and found a strongly peaked distribution at 7.8 ± 1.2 M_☉. They used two models for the mass function: a Gaussian and a decaying exponential with a minimum "turn-on" mass (motivated by the analytical model of the black hole mass function in Fryer & Kalogera 2001). We note that Ozel et al. (2010) do not provide confidence limits for the minimum black hole mass, instead discussing only the model parameters at the peak of their posterior distributions. They also do not perform any model selection analysis; thus, they give the distribution of parameters within each of their models, but cannot say which model is more likely to correspond to the true distribution of black hole masses. Nevertheless, it appears that their analysis confirms the existence of a mass gap. Ozel et al. (2010) discuss possible selection effects that could lead to the appearance of a mass gap, but conclude that these effects could not produce the observed gap, which they therefore claim is a real property of the black hole mass distribution.

We use a Bayesian Markov Chain Monte Carlo (MCMC) analysis to quantitatively assess a wide range of models for the black hole mass function for both samples. We include both parametric models, such as a Gaussian, and non-parametric models where the mass function is represented by histograms with various numbers of bins. (Our set of models includes those of Ozel et al. 2010 and Bailyn et al. 1998.) After computing posterior distributions for the model parameters, we use model selection techniques (including a new technique for efficient reversible-jump MCMC; Farr & Mandel 2011) to compare the evidence for the various models from both samples.

We define the "minimum black hole mass" to be the 1% quantile, $M_{1\%}$, in the black hole mass distribution (see Section 5). In qualitative agreement with Ozel et al. (2010) and Bailyn et al. (1998), we find strong evidence for a mass gap among the best models for both samples. Our analysis gives distributions for $M_{1\%}$ implied by the data in the context of each of our models for the black hole mass distribution. In the context of the best model for the low-mass systems (a power law), the distribution for $M_{1\%}$ gives $M_{1\%} > 4.3\,M_\odot$ with 90% confidence; in the context of the best model for the combined sample of lower- and high-mass systems the distribution of $M_{1\%}$ has $M_{1\%} > 4.5\,M_\odot$ with 90% confidence. Further, in the context of models with lower evidence, most also have a mass gap, with 90% confidence bounds on $M_{1\%}$ significantly above a "safe" maximum neutron star mass of 3 M_☉ (Kalogera & Baym 1996).

We find that, for the low-mass X-ray binary sample, the theoretical model from Fryer & Kalogera (2001)—a decaying exponential—is strongly disfavored by our model selection. We find that the low-mass systems are best described by a power law, followed closely by a Gaussian (which is the second model considered by Ozel et al. 2010). On the other hand, we find that the theoretical model from Fryer & Kalogera (2001) is the preferred model for the combined sample of low- and high-mass X-ray binaries. A model with two separate Gaussian peaks also has relatively high evidence for the combined sample of systems. The difference in best-fitting model indicates that the low-mass subsample is not consistent with being drawn from the distribution of the combined population.

The structure of this paper is as follows. In Section 2, we discuss the 15 systems that comprise the low-mass X-ray binary black hole sample and the 5 additional high-mass, wind-fed systems that make up the combined sample. In Section 3, we discuss the Bayesian techniques we use to analyze the black hole mass distribution, the techniques we use for model selection, and the parametric and non-parametric models we will use for the black hole mass distribution. In Section 4, we discuss the results of our analysis and model selection. In Section 5, we discuss the distribution of the minimum black hole mass implied by the analysis of Section 4. In Section 6, we summarize our results and comment on the significance of the observed mass gap in the context of theoretical models. Appendix A describes MCMC techniques in some detail. Appendix B explains our novel algorithm for efficiently performing the reversible-jump MCMC computations used in the model comparisons of Section 4 (but see also Farr & Mandel 2011).

The 20 X-ray binary systems on which this study is based are listed in Table 1. We separate the systems into 15 low-mass systems in which the central black hole appears to be fed by Roche lobe overflow from the secondary, and 5 high-mass systems in which the black hole is fed via winds (these systems all have a secondary that appears to be more massive than the black hole). The low- and high-mass systems undoubtedly have different evolutionary tracks, and therefore it is reasonable that they would have different black hole mass distributions. We will first analyze the 15 low-mass systems alone (Section 4.1), and then the combined sample of 20 systems (Section 4.2).

In each of these systems, spectroscopic measurements of the secondary star provide an orbital period for the system and a semi-amplitude for the secondary's velocity curve. These measurements can be combined into the mass function,

$$\begin{equation} f(M) = \frac{{\it P K}^3}{2\pi G} = \frac{M \sin ^3 i}{\left(1 + q \right)^2}, \end{equation} \tag{ 1 }$$

where P is the orbital period, K is the secondary's velocity semi-amplitude, M is the black hole mass, i is the inclination of the system, and q ≡ M₂/M is the mass ratio of the system.

The mass function defines a lower limit on the mass: f(M) < M. To accurately determine the mass of the black hole, the inclination i and mass ratio q must be measured. Ideally, this can be accomplished by fitting ellipsoidal light curves and study of the rotational broadening of spectral lines from the secondary, but even in the most studied case (see, e.g., Cantrell et al. 2010 on A0620) this procedure is complicated. In particular, contributions from an accretion disk and hot spots in the disk can significantly distort the measured inclination and mass ratios. For some systems (e.g., GS 1354; Casares et al. 2009) strong variability completely prevents determination of the inclination from the light curve; in these cases an upper limit on the inclination often comes from the observed lack of eclipses in the light curve. In general, accurately determining q and i requires a careful system-by-system analysis.

For the purposes of this paper, we adopt the following simplified approach to the estimation of the black hole mass from the observed data. When an observable is well constrained, we assume that the true value is normally distributed about the measured value with a standard deviation equal to the quoted observational error. This is the case for the mass function in all the systems we use, and for many systems' mass ratios and inclinations. When a large range is quoted in the literature for an observable, we take the true value to be distributed uniformly (for the mass ratio) or isotropically (for the inclination) within the quoted range. Table 1 gives the assumed distribution for the observables in the 20 systems we use. We do not attempt to deal with the systematic biases in the observational determination of f, q, and i in any realistic way; we are currently investigating more realistic treatments of the errors (including observational biases that can shift the peak of the true mass distribution away from the "best-fit" mass in the observations). This treatment will appear in future work.

From these assumptions, we can generate probability distributions for the true mass of the black hole given the observations and errors via the Monte Carlo method: drawing samples of f, q, and i from the assumed distributions and computing the mass implied by Equation (1) gives samples of M from the distribution induced by the relationship in Equation (1). Mass distributions generated in this way for the systems used in this work are shown in Figures 1 and 2. Systems for which i is poorly constrained have broad "tails" on their mass distributions. These mass distributions constitute the "observational data" we will use in the remainder of this paper.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this section, we describe the statistical analysis we will apply to various models for the underlying mass distribution from which the low-mass sample and the combined sample of X-ray binary systems in Table 1 were drawn. The results of our analysis are presented in Section 4.

3.1. Bayesian Inference

The end result of our statistical analysis will be the probability distribution for the parameters of each model implied by the data from Section 2 in combination with our prior assumptions about the probability distribution for the parameters. Bayes' rule relates these quantities. For a model with parameters $\vec{\theta }$ in the presence of data d, Bayes' rule states

$$\begin{equation} p(\vec{\theta }| d) = \frac{p(d | \vec{\theta }) p(\vec{\theta })}{p(d)}. \end{equation} \tag{ 2 }$$

Here, $p(\vec{\theta }|d)$, called the posterior probability distribution function, is the probability distribution for the parameters $\vec{\theta }$ implied by the data d; $p(d|\vec{\theta })$, called the likelihood, is the probability of observing data d given that the model parameters are $\vec{\theta }$; $p(\vec{\theta })$, called the prior, reflects our estimate of the probability of the various model parameters in the absence of any data; and p(d), called the evidence, is an overall normalizing constant ensuring that

$$\begin{equation} \int d\theta \, p(\vec{\theta }|d) = 1, \end{equation} \tag{ 3 }$$

whence

$$\begin{equation} p(d) = \int d\vec{\theta }\, p(d|\vec{\theta }) p(\vec{\theta }). \end{equation} \tag{ 4 }$$

In our context, the data are the mass distributions given in Section 2: d = {p_i(M)|i = 1, 2, ..., 20}. We assume that the measurements in Section 2 are independent, so the complete likelihood is given by a product of the likelihoods for the individual measurements. For a model with parameters $\vec{\theta }$ that predicts a mass distribution $p(M|\vec{\theta })$ for black holes, we have

$$\begin{equation} p(d|\vec{\theta }) = \prod _i \int {\it dM}\, p_i(M) p(M|\vec{\theta }). \end{equation} \tag{ 5 }$$

That is, the likelihood of an observation is the average over the individual mass distribution implied by the observation, p_i(M), of the probability for a black hole of that mass to exist according to the model of the mass distribution, $p(M | \vec{\theta })$. We approximate the integrals as averages of $p(M|\vec{\theta })$ over the Monte Carlo mass samples drawn from the distributions in Table 1 (also see Figures 1 and 2):

$$\begin{equation} p(d|\vec{\theta }) \approx \prod _i \frac{1}{N_i} \sum _{j = 1}^{N_i} p(M_{ij} | \vec{\theta }), \end{equation} \tag{ 6 }$$

where M_ij is the jth sample (out of a total N_i) from the ith individual mass distribution.

Our calculation of the likelihood of each observation does not include any attempt to account for selection effects in the observations. We simply assume (almost certainly incorrectly) that any black hole drawn from the underlying mass distribution is equally likely to be observed. The results of Ozel et al. (2010) suggest that selection effects are unlikely to significantly bias our analysis.

For a mass distribution with several parameters, $p(\vec{\theta }| d)$ lives in a multi-dimensional space. Previous works (Ozel et al. 2010; Bailyn et al. 1998) have considered models with only two parameters; for such models evaluating $p(\vec{\theta }|d)$ on a grid may be a reliable method. Many of our models for the underlying mass distribution have three or more parameters. Exploring the entirety of these parameter spaces with a grid rapidly becomes prohibitive as the number of parameters increases. A more efficient way to explore the distribution $p(\vec{\theta }| d)$ is to use a MCMC method (see Appendix A). MCMC methods produce a chain (sequence) of parameter samples, $\lbrace \vec{\theta }_i \, | \, i = 1, \ldots \rbrace$, such that a particular parameter sample, $\vec{\theta }$, appears in the sequence with a frequency proportional to its posterior probability, $p(\vec{\theta }|d)$. In this way, regions of parameter space where $p(\vec{\theta }|d)$ is large are sampled densely while regions where $p(\vec{\theta }|d)$ is small are effectively ignored.

Once we have a chain of samples from $p(\vec{\theta }|d)$, the distribution for any quantity of interest can be computed by evaluating it on each sample in the chain and forming a histogram of these values. For example, to compute the one-dimensional distribution for a single parameter obtained by integrating over all other dimensions in parameter space, called the "marginalized" distribution, one plots the histogram of the values of that parameter appearing in the chain.

3.2. Priors

3.3. Parametric Models for the Black Hole Mass Distribution

Here, we discuss the various parametric models of the underlying black hole mass distribution considered in this paper.

3.3.1. Power-law Models

Many astrophysical distributions are power laws. Let us assume that the black hole mass distribution is given by

$$\begin{eqnarray} p(M|\vec{\theta }) &=& p(M|\lbrace M_{\rm{min}}, M_{\rm{max}}, \alpha \rbrace)\nonumber\\ &=& \left\lbrace \begin{array}{@{}l@{\quad }l@{}}A M^\alpha & M_{\rm{min}}\le m \le M_{\rm{max}}\\ \ms 0 & \textnormal {otherwise} \end{array}\right.. \end{eqnarray} \tag{ 7 }$$

The normalizing constant A is

$$\begin{equation} A = \frac{1+\alpha }{M_{\rm{max}}^{1+\alpha } - M_{\rm{min}}^{1+\alpha }}. \end{equation} \tag{ 8 }$$

We choose uniform priors on M_min and M_max ⩾ M_min between 0 and 40 M_☉, and uniform priors on the exponent α in a broad range between −15 and 13:

$$\begin{eqnarray} &&p(\vec{\theta }) = p(\lbrace M_{\rm{min}}, M_{\rm{max}}, \alpha \rbrace) \nonumber\\ &&= \left\lbrace \begin{array}{@{}l@{\quad }l@{}}2 \frac{1}{40^2}\frac{1}{28} & 0 \le M_{\rm{min}}\le M_{\rm{max}}\le 40, \quad -15 \le \alpha \le 13 \\ \ms 0 & \textnormal {otherwise} \end{array}\right..\nonumber\\ \end{eqnarray} \tag{ 9 }$$

Our MCMC analysis output is a list of {M_min, M_max, α} values distributed according to the posterior

$$\begin{eqnarray} p(\vec{\theta }|d) &=& p(\lbrace M_{\rm{min}}, M_{\rm{max}}, \alpha \rbrace |d)\nonumber\\ & \propto & p(d|\lbrace M_{\rm{min}}, M_{\rm{max}}, \alpha \rbrace) p(\lbrace M_{\rm{min}}, M_{\rm{max}}, \alpha \rbrace),\quad\quad \end{eqnarray} \tag{ 10 }$$

with the likelihood p(d|{M_min, M_max, α}) defined in Equation (5).

3.3.2. Decaying Exponential

Fryer & Kalogera (2001) studied the relation between progenitor and remnant mass in simulations of supernova explosions. Combining this with the mass function for supernova progenitors, they suggested that the black hole mass distribution may be well represented by a decaying exponential with a minimum mass:

$$\begin{eqnarray} p(M|\vec{\theta }) &=& p(M|\lbrace M_{\rm{min}}, M_0\rbrace)\nonumber\\ & =& \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\frac{e^{\frac{M_{\rm{min}}}{M_0}}}{M_0} \exp \left[ - \frac{M}{M_0} \right] & M \ge M_{\rm{min}}\\ \ms 0 & \textnormal {otherwise} \end{array}\right.. \end{eqnarray} \tag{ 11 }$$

We choose a prior for this model where M_min is uniform between 0 and 40 M_☉. For each M_min, we choose M₀ uniformly within a range ensuring that 40 M_☉ is at least two scale masses above the cutoff: 40 M_☉ ⩾ M_min + 2M₀. This ensures that the majority of the mass probability lies in the range 0 M_☉ ⩽ M ⩽ 40 M_☉. The resulting prior is

$$\begin{eqnarray} && p(\vec{\theta }) = p(\lbrace M_{\rm{min}}, M_0\rbrace)\nonumber\\ && = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\frac{4}{40^2} & 0 \le M_{\rm{min}}\le 40, \quad 0 < M_0, \quad M_{\rm{min}}+ 2 M_0 \le 40 \\ \ms 0 & \textnormal {otherwise} \end{array}\right..\nonumber\\ \end{eqnarray} \tag{ 12 }$$

3.3.3. Gaussian and Two-Gaussian Models

The mass distributions in Figure 1 all peak in a relatively narrow range near ∼10 M_☉. The prototypical single-peaked probability distribution is a Gaussian:

$$\begin{equation} p(M|\vec{\theta }) = p(M|\lbrace \mu, \sigma \rbrace) = \frac{1}{\sigma \sqrt{2\pi }} \exp \left[ - \left(\frac{M - \mu }{\sqrt{2} \sigma } \right)^2 \right]. \end{equation} \tag{ 13 }$$

We use a prior on the mean mass, μ, and the standard deviation, σ, that ensures that the majority of the mass distribution lies below 40 M_☉:

$$\begin{equation} p(\lbrace \mu,\sigma \rbrace) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\frac{8}{40^2} & 0 \le \mu \le 40, \quad \sigma \ge 0, \quad \mu + 2\sigma \le 40 \\ \ms 0 & \textnormal {otherwise} \end{array}\right., \end{equation} \tag{ 14 }$$

where both μ and σ are measured in solar masses.

Though we do not expect to find a second peak in the low-mass distribution, we may find evidence of one when exploring the combined low- and high-mass samples. To look for a second peak in the black hole mass distribution, we use a two-Gaussian model:

$$\begin{eqnarray} && p(M|\vec{\theta }) = p(M|\lbrace \mu _1, \mu _2, \sigma _1, \sigma _2, \alpha \rbrace) = \frac{\alpha }{\sigma _1 \sqrt{2\pi }}\nonumber\\&& \quad \times \exp \left[ - \left(\frac{M - \mu _1}{\sqrt{2}\sigma _1} \right)^2 \right] + \frac{1-\alpha }{\sigma _2 \sqrt{2\pi }} \exp \left[ - \left(\frac{M - \mu _2}{\sqrt{2}\sigma _2} \right)^2 \right].\nonumber\\ \end{eqnarray} \tag{ 15 }$$

The probability is a linear combination of two Gaussians with weights α and 1 − α. We restrict μ₁ < μ₂ and also impose combined conditions on μ_i and σ_i that ensure that most of the mass probability lies below 40 M_☉ with the prior

$$\begin{eqnarray} && p(\lbrace \mu _1, \mu _2, \sigma _1, \sigma _2, \alpha \rbrace)\nonumber\\ && = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}2 p(\lbrace \mu _1, \sigma _1\rbrace) p(\lbrace \mu _2, \sigma _2\rbrace) & \mu _1 \le \mu _2, \quad 0 \le \alpha \le 1 \\ \ms 0 & \textnormal {otherwise} \end{array}\right.,\quad\quad \end{eqnarray} \tag{ 16 }$$

where the single-Gaussian prior, p({μ_i, σ_i}), is defined in Equation (14).

3.3.4. Log Normal

Many of the mass distributions for the systems in Figure 1 rise rapidly to a peak and then fall off more slowly in a longer tail toward high masses. So far, none of the parameterized distributions we have discussed have this property. In this section, we consider a log-normal model for the underlying mass distribution; the log-normal distribution has a rise to a peak with a slower falloff in a long tail.

The log-normal distribution gives log M a Gaussian distribution with mean μ and standard deviation σ:

$$\begin{eqnarray} &&p(M|\vec{\theta }) = p(M|\lbrace \mu, \sigma \rbrace) = \frac{1}{ \sqrt{2\pi } M \sigma } \exp \left[ -\frac{\left(\log M - \mu \right)^2}{2 \sigma ^2} \right].\nonumber\\ \end{eqnarray} \tag{ 17 }$$

The parameters μ and σ are dimensionless; the mean mass 〈M〉 and mass standard deviation σ_M are related to μ and σ by

$$\begin{eqnarray} \langle M \rangle & = & \exp \left(\mu + \frac{1}{2} \sigma ^2 \right) \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} \sigma _M & = & \langle M \rangle \sqrt{\exp \left(\sigma ^2 \right) - 1}. \end{eqnarray} \tag{ 19 }$$

For a fair comparison with the other models, we impose a prior that is flat in 〈M〉 and σ_M. To ensure that most of the probability in this model occurs for masses below 40 M_☉, we require 〈M〉 + 2σ_M ⩽ 40, resulting in a prior

$$\begin{eqnarray} p(\vec{\theta }) &=& p(\lbrace \mu,\sigma \rbrace)\nonumber\\ & =& \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\frac{4}{40^2} \left| \frac{\partial \left(\langle M \rangle, \sigma _M \right)}{\partial \left(\mu, \sigma \right)} \right| & \sigma > 0, \quad \langle M \rangle + 2 \sigma _M \le 40 \\ \ms 0 & \textnormal {otherwise} \end{array}\right.,\quad\quad \end{eqnarray} \tag{ 20 }$$

where

$$\begin{equation} \left| \frac{\partial \left(\langle M \rangle, \sigma _M \right)}{\partial \left(\mu, \sigma \right)} \right| = \frac{\exp (2 (\mu + \sigma ^2)) \sigma }{\sqrt{\exp (\sigma ^2) - 1}} \end{equation} \tag{ 21 }$$

is the determinant of the Jacobian of the map in Equations (18) and (19).

3.4. Non-parametric Models for the Black Hole Mass Distribution

The previous subsection discussed models for the underlying black hole mass distribution that assumed particular parameterized shapes for the distribution. In this subsection, we will discuss models that do not assume a priori a shape for the black hole mass distribution. The fundamental non-parametric distribution in this section is a histogram with some number of bins, N_bin. Such a distribution is piecewise constant in M.

One choice for representing such a histogram would be to fix the bin locations, and allow the heights to vary. With this approach, one should be careful not to "split" features of the mass distribution across more than one bin in order to avoid diluting the sensitivity to such features; similarly, one should avoid including more than "one" feature in each bin. The locations of the bins, then, are crucial. An alternative representation of histogram mass distributions avoids this difficulty.

We choose to represent a histogram mass distribution with N_bin bins by allocating a fixed probability, 1/N_bin, to each bin. The lower and upper bounds for each bin are allowed to vary; when these are close to each other (i.e., the bin is narrow), the distribution will have a large value, and conversely when the bounds are far from each other. We assume that the non-zero region of the distribution is contiguous, so we can represent the boundaries of the bins as a non-decreasing array of masses, $w_0 \le w_1 \le \ldots \le w_{N_{\rm{bin}}}$, with w₀ the minimum and $w_{N_{\rm{bin}}}$ the maximum mass for which the distribution has support. This gives the distribution

$$\begin{eqnarray} p(M|\theta) &=& p(M|\lbrace w_0, \ldots, w_{N_{\rm{bin}}}\rbrace)\nonumber\\ & =& \left\lbrace \begin{array}{@{}l@{\quad }l@{}}0 & M < w_0 \textnormal { or } w_{N_{\rm{bin}}} \le M \\ \ms \frac{1}{N_{\rm{bin}}} \frac{1}{w_{i+1} - w_i} & w_i \le M < w_{i+1} \end{array}\right..\quad\quad \end{eqnarray} \tag{ 22 }$$

For priors on the histogram model with N_bin bins, we assume that the bin boundaries are uniformly distributed between 0 and 40 M_☉ subject only to the constraint that the boundaries are non-decreasing from w₀ to $w_{N_{\rm{bin}}}$:

$$\begin{eqnarray} && p(\lbrace w_0, \ldots, w_{N_{\rm{bin}}}\rbrace)\nonumber\\ && = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\frac{\left(N_{\rm{bin}}+1\right)!}{40^{N_{\rm{bin}}+1}} & 0 \le w_0 \le w_1 \le \ldots \le w_{N_{\rm{bin}}} \le 40 \\ \ms 0 & \textnormal {otherwise} \end{array}\right..\quad \end{eqnarray} \tag{ 23 }$$

We consider histograms with up to five bins in this work. We will see that the evidence for the histogram models (see Sections 3.5, 4.1.7, and 4.2.7) from both the low-mass and combined data sets is decreasing as the number of bins reaches five, indicating that increasing the number of bins beyond five would not sufficiently improve the fit to the mass distribution to compensate for the extra parameter space volume implied by the additional parameters.

3.5. Bayesian Model Selection

In Sections 3.3 and 3.4, we discussed a series of models for the underlying black hole mass distribution. Our MCMC analysis will provide the posterior distribution of the parameters within each model, but does not tell us which models are more likely to correspond to the actual distribution. This model selection problem is the topic of this section.

Consider a set of models, {M_i|i = 1, ...}, each with corresponding parameters $\vec{\theta }_i$. Re-writing Equation (2) to be explicit about the assumption of a particular model, we have

$$\begin{equation} p(\vec{\theta }_i | d, M_i) = \frac{p(d|\vec{\theta }_i, M_i) p(\vec{\theta }_i | M_i)}{p(d|M_i)}. \end{equation} \tag{ 24 }$$

This gives the posterior probability of the parameters $\vec{\theta }_i$ in the context of model M_i. But the model itself can be regarded as a discrete parameter in a larger "super-model" that encompasses all the M_i. The parameters for the super-model are $\lbrace M_i, \vec{\theta }_i\rbrace$: a choice of model and the corresponding parameter values within that model. Each point in the super-model parameter space is a statement that, e.g., "the underlying mass distribution is a Gaussian, with parameters μ and σ," or "the underlying mass distribution is a triple-bin histogram with parameters w₁, w₂, w₃, and w₄," or... The posterior probability of the super-model parameters is given by Bayes' rule:

$$\begin{equation} p(\vec{\theta }_i, M_i|d) = \frac{p(d|\vec{\theta }_i, M_i) p(\vec{\theta }_i |M_i) p(M_i)}{p(d)}, \end{equation} \tag{ 25 }$$

where we have introduced the model prior p(M_i), which represents our estimate on the probability that model M_i is correct in the absence of the data d. The normalizing evidence is now

$$\begin{eqnarray} p(d) &=& \sum _i \int d\vec{\theta }_i\, p(d|\vec{\theta }_i, M_i) p(\vec{\theta }_i |M_i) p(M_i)\nonumber\\ & =& \sum _i p(d|M_i) p(M_i), \end{eqnarray} \tag{ 26 }$$

writing the single-model evidence from Equation (4) as p(d|M_i) to be explicit about the dependence on the choice of model.

To compare the various models M_i, we are interested in the marginalized posterior probability of M_i:

$$\begin{equation} p(M_i|d) \equiv \int d\vec{\theta }_i\, p(\vec{\theta }_i, M_i|d). \end{equation} \tag{ 27 }$$

This is the integral of the posterior over the entire parameter space of model M_i. The marginalized posterior probability of model M_i can be rewritten in terms of the single-model evidence, p(d|M_i) (see Equations (25) and (4)):

$$\begin{eqnarray} p(M_i|d) &=& \int d\vec{\theta }_i\, p(\vec{\theta }_i, M_i|d) = \frac{p(M_i)}{p(d)}\nonumber\\ && \times \int d\vec{\theta }_i p(d|\vec{\theta }_i,M_i) p(\vec{\theta }_i|M_i) = \frac{p(d|M_i) p(M_i)}{p(d)}.\nonumber\\ \end{eqnarray} \tag{ 28 }$$

Here and throughout, we assume that any of the models in Section 3 are equally likely a priori, so the model priors are equal:

$$\begin{equation} p(M_i) = \frac{1}{N_{\rm{model}}}. \end{equation} \tag{ 29 }$$

A powerful technique⁶ for computing p(M_i|d) is the reversible-jump MCMC (Green 1995). Reversible-jump MCMC, discussed in more detail in Appendix B, is a standard MCMC analysis conducted in the super-model. The result of a reversible-jump MCMC is a chain of samples, $\lbrace M_i, \vec{\theta }_i\, | \, i = 1, \ldots \rbrace$, from the super-model parameter space. The integral in Equation (28) can be estimated by counting the number of times that a given model M_i appears in the reversible-jump MCMC chain:

$$\begin{equation} p(M_i|d) = \int d\vec{\theta }_i p(M_i, \vec{\theta }_i|d) \approx \frac{N_i}{N}, \end{equation} \tag{ 30 }$$

where N_i is the number of MCMC samples that have discrete parameter M_i and N is the total number of samples in the MCMC.

Naively implemented reversible-jump MCMCs can be very inefficient when the posteriors for a model or models are strongly peaked. In this circumstance, a proposed MCMC jump into one of the peaked models is unlikely to land on the peak by chance; since it is rare to propose a jump into the important regions of parameter space of the peaked model in a naive reversible-jump MCMC, the output chain must be very long to ensure that all models have been compared fairly. We describe a new algorithm in Appendix B that produces very efficient jump proposals for a reversible-jump MCMC by exploiting the information about the model posteriors we have from the single-model MCMC samples. (See also Farr & Mandel 2011.) With this algorithm, reasonable chain lengths can fairly compare all the models under consideration. We have used this algorithm to perform 10-way reversible-jump MCMCs to calculate the relative evidence for both the parametric and non-parametric models in this study. These results appear in Section 4.

In this section, we discuss the results of our MCMC analysis of the posterior distributions of parameters for the models in Sections 3.3 and 3.4. We also discuss model selection results. The results in Section 4.1 apply to the low-mass sample of systems, while those of Section 4.2 apply to the combined sample of systems.

4.1. Low-mass Systems

Table 2 gives quantiles of the marginalized parameter distributions of the parametric models implied by the low-mass data. Table 3 gives the quantiles of the histogram bin boundaries in the non-parametric analysis implied by the low-mass data.

Table 2. Quantiles of the Marginalized Distribution for Each of the Parameters in the Models Discussed in Section 3.3 Implied by the Low-mass Data

Model	Parameter	5%	15%	50%	85%	95%
Power law (Equation (7))	Mmin	1.2786	4.1831	6.1001	6.5011	6.6250
	Mmax	8.5578	8.9214	23.3274	36.0002	38.8113
	α	−12.4191	−10.1894	−6.3861	2.8476	5.6954
Exponential (Equation (11))	Mmin	5.0185	5.4439	6.0313	6.3785	6.5316
	M0	0.7796	0.9971	1.5516	2.4635	3.2518
Gaussian (Equation (13))	μ	6.6349	6.9130	7.3475	7.7845	8.0798
	σ	0.7478	0.9050	1.2500	1.7335	2.1134
Two Gaussian (Equation (15))	μ1	5.4506	6.3877	7.1514	7.6728	7.9803
	μ2	7.2355	7.7387	12.3986	25.2456	31.4216
	σ1	0.3758	0.7626	1.2104	1.7981	2.3065
	σ2	0.2048	0.6421	1.9182	5.2757	7.2625
	α	0.0983	0.3526	0.8871	0.9792	0.9936
Log normal (Equation (17))	〈M〉	6.7619	7.0122	7.4336	7.9159	8.2942
	σM	0.7292	0.8920	1.2704	1.8695	2.4069

Notes. We indicate the 5%, 15%, 50% (median), 85%, and 95% quantiles. The marginalized distribution can be misleading when there are strong correlations between variables. For example, while the marginalized distributions for the power-law parameters are quite broad, the distribution of mass distributions implied by the power-law MCMC samples is similar to the other models. This occurs in spite of the broad marginalized distributions because of the correlations between the slope and limits of the power law discussed in Section 3.3.1.

Download table as:  ASCIITypeset image

Table 3. The 5%, 15%, 50% (Median), 85%, and 95% Quantiles for the Bin Boundaries in the One-through Five-bin Histogram Models Discussed in Section 3.4

Bins	Boundary	5%	15%	50%	85%	95%
1	w0	3.94488	4.55603	5.43333	6.02557	6.29749
	w1	8.50844	8.69262	9.11784	9.83477	10.5128
2	w0	3.3426	4.2047	5.39132	6.18413	6.47553
	w1	6.41972	6.72605	7.43421	8.2489	8.52885
	w2	8.46161	8.65077	9.12694	10.1113	11.2595
3	w0	2.18176	3.54345	5.16094	6.16473	6.44697
	w1	5.68876	6.14223	6.68829	7.38725	8.04235
	w2	6.8297	7.22718	8.1451	8.7512	9.27296
	w3	8.44307	8.67362	9.25718	12.1688	21.92
4	w0	1.32131	2.7934	4.66156	5.78459	6.17946
	w1	5.20112	5.77331	6.42501	6.98427	7.44584
	w2	6.41805	6.73535	7.43826	8.32958	8.64212
	w3	7.40302	7.95608	8.58976	9.33897	10.3992
	w4	8.56724	8.8059	10.2451	24.3573	34.2423
5	w0	0.9392	2.28789	4.33389	5.7012	6.21166
	w1	4.69778	5.44302	6.26575	6.76407	7.14427
	w2	6.1388	6.47155	7.00606	7.97325	8.38259
	w3	6.82058	7.28677	8.22514	8.81555	9.41012
	w4	8.02335	8.36993	8.94879	11.3206	17.3349
	w5	8.7112	9.25208	16.2059	31.897	37.2738

Download table as:  ASCIITypeset image

Recall that each MCMC sample in our analysis gives the parameters for a model of the black hole mass distribution. The chain of samples of parameters for a particular model gives us a distribution of black hole mass distributions. Figure 3 gives a sense of the shape and range of the distributions of black hole mass distributions that result from our MCMC analysis. In Figure 3, we plot the median, 10%, and 90% values of the black hole mass distributions that result from the MCMC chains. Because the choice of parameters that gives, for example, the median distribution value at one mass need not give the median distribution at another mass, these curves do not necessarily look like the underlying model for the mass distribution. For the same reason, they are not necessarily normalized.

Zoom In Zoom Out Reset image size

4.1.1. Power Law

In Figure 4, we display a histogram of the resulting samples in each of the parameters M_min, M_max, and α for the power-law model (see Equation (7)); this represents the one-dimensional "marginalized" distribution

$$\begin{equation} p(\alpha |d) = \int dM_{\rm{min}}\, dM_{\rm{max}}\, p(\lbrace M_{\rm{min}}, M_{\rm{max}}, \alpha \rbrace |d), \end{equation} \tag{ 31 }$$

and similarly for M_min and M_max.

Zoom In Zoom Out Reset image size

The marginalized distribution for α is broad, with

$$\begin{equation} -11.8 < \alpha < 6.8 \end{equation} \tag{ 32 }$$

enclosing 90% of the probability (excluding 5% on each side). We have p(α < 0) = 0.6. The median value is α = −3.35. The broadness of the marginalized distribution for α comes from the need to match the relatively narrow range in mass of the low-mass systems. When α is negative, the resulting mass distribution slopes down; M_min is constrained to be near the lowest mass of the observed black holes, while M_max is essentially irrelevant. Conversely, when α is positive and the mass distribution slopes up, M_max must be close to the largest mass observed, while M_min is essentially irrelevant. Figure 5 illustrates this effect, showing the correlations between α and M_min and α and M_max. When we include the high-mass systems in the analysis, the long tail will eliminate this effect by bringing both M_min and M_max into play for all values of α.

Zoom In Zoom Out Reset image size

4.1.2. Decaying Exponential

Figure 6 displays the marginalized posterior distribution for the scale mass of the exponential, M₀, and the cutoff mass, M_min (see Equation 11). The median scale mass is M₀ = 1.55, and 0.78 ⩽ M₀ ⩽ 3.25 with 90% confidence. This model was one of those considered by Ozel et al. (2010), whose results (M₀ ∼ 1.5 and M_min ∼ 6.5) are broadly consistent with ours. Figure 7 displays the MCMC samples in the M_min, M₀ plane for this model. There is a small correlation between smaller M_min and larger M₀, which is driven by the need to widen the distribution to encompass the peak of the mass measurements in Figure 1 when the minimum mass is smaller.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.1.3. Gaussian

Figure 8 shows the resulting marginalized distributions for the parameters μ and σ. We constrain the peak of the Gaussian between 6.63 ⩽ μ ⩽ 8.08 with 90% confidence. This model also appeared in Ozel et al. (2010); they found μ ∼ 7.8 and σ ∼ 1.2, consistent with our results here.

Zoom In Zoom Out Reset image size

4.1.4. Two Gaussian

Figure 9 shows the marginalized distributions for the two-Gaussian model parameters from our MCMC runs. We find α > 0.8 with 62% probability, clearly favoring the Gaussian with smaller mean. The distributions for μ₁ and σ₁ are similar to those of the single Gaussian displayed in Figure 8, indicating that this Gaussian is centered around the peaks of the low-mass distributions. The second Gaussian's parameter distributions are much broader. The second Gaussian appears to be sampling the tail of the mass samples. In spite of the extra degrees of freedom in this model, we find that this model is strongly disfavored relative to the single-Gaussian model for this data set: $p(\textnormal {Gaussian}|d) / p(\textnormal {two Gaussian}|d) \simeq 4.7$ (see Sections 3.5 and 4.1.7 for discussion).

Zoom In Zoom Out Reset image size

4.1.5. Log Normal

The marginal distributions for 〈M〉 and σ_M appear in Figure 10. The distributions are similar to those for μ and σ from the Gaussian model in Section 3.3.3.

Zoom In Zoom Out Reset image size

4.1.6. Histogram Models

The median values of the histogram mass distributions that result from the MCMC samples of the posterior distribution for the w_i parameters for one-, two-, three-, four-, and five-bin histogram models are shown in Figure 3. Table 3 gives quantiles of the marginalized bin boundary distributions for the histogram models.

As the number of bins increases, the models are better able to capture features of the mass distribution, but we find that the one-bin histogram is the most probable of the histogram models for the low-mass data (see Section 4.1.7 for discussion). This occurs because the extra fitting power does not sufficiently improve the fit to compensate for the vastly larger parameter space of the models with more bins.

4.1.7. Model Selection for the Low-mass Sample

We have performed a suite of 500 independent reversible-jump MCMCs jumping between all the models (both parametric and non-parametric) described in Section 3 using the single-model MCMC samples to construct an efficient jump proposal for each model as described above (see Appendix B). The numbers of counts in each model are consistent across the MCMCs in the suite; Figure 11 displays the average probability for each model across the suite, along with the 1σ errors on the average inferred from the standard deviation of the model counts across the suite. Table 4 gives the numerical values of the average probability for each model across the suite of MCMCs.

Zoom In Zoom Out Reset image size

Table 4. Relative Probabilities of the Various Models from Section 3 Implied by the Low-mass Data

Model	Relative Evidence
Power law (Section 3.3.1)	0.331488
Gaussian (Section 3.3.3)	0.288129
Log normal (Section 3.3.4)	0.138435
Exponential (Section 3.3.2)	0.0916218
Two Gaussian (Section 3.3.3)	0.0662577
Histogram (1 Bin, Section 3.4)	0.0641941
Histogram (2 Bin, Section 3.4)	0.015184
Histogram (3 Bin, Section 3.4)	0.00332933
Histogram (4 Bin, Section 3.4)	0.000999976
Histogram (5 Bin, Section 3.4)	0.0003614

Notes. These probabilities have been computed from reversible-jump MCMC samples using the efficient jump proposal algorithm in Appendix B. (See also Figure 11.)

Download table as:  ASCIITypeset image

The most favored model is the power law from Section 3.3.1, followed by the Gaussian model from Section 3.3.3. Interestingly, the theoretical curve from Fryer & Kalogera (2001; the exponential model of Section 3.3.2) places fourth in the ranking of evidence.

Though the model probabilities presented in this section have small statistical error, they are subject to large "systematic" error. The source of this error is both the particular choice of model prior (uniform across models) and the choice of priors on the parameters within each model used for this work. For example, the theoretically preferred exponential model (Section 3.3.2) is only a factor of ∼3 away from the power-law model (Section 3.3.1), which does not have theoretical support. Is such support worth a factor of three in the model prior? Alternately, we may say we know (in advance of any mass measurements) that black holes must exist with mass ≲ 10 M_☉; then we could, for example, impose a prior on the minimum mass in the exponential model (M_min) that is uniform between 0 and 10 M_☉, which would reduce the prior volume available for the model by a factor of four without significantly reducing the posterior support for the model. This has the same effect as increasing the model prior by a factor of four, which would move this model from fourth to first place. Of course, we would then have to modify the prior support for the other models to take into account the restriction that there must be black holes with M ≲ 10 M_☉... Linder & Miquel (2008) discuss these issues in the context of cosmological model selection, concluding with a warning against overreliance on model selection probabilities.

Nevertheless, we believe that our model comparison is reasonably fair (see the discussion of priors in Section 3.2). It seems safe to conclude that "single-peaked" models (the power law and Gaussian) are preferred over "extended" models (the exponential or log normal), or those with "structure" (the many-bin histograms or two-Gaussian model). Previous studies have also supported the "single, narrow peak" mass distribution (Bailyn et al. 1998; Ozel et al. 2010). In this light, poor performance of the single-bin histogram is surprising.

4.2. Combined Sample

This section repeats the analysis of the models from Section 3, but including the high-mass, wind-fed systems from Table 1 (see also Figure 2) in the sample. Figure 12 displays bounds on the value of the underlying mass distribution for the various models in Section 3 applied to this data set; compare to Figure 3. The inclusion of the high-mass, wind-fed systems tends to widen the distribution toward the high-mass end and, in models that allow it, produces a second, high-mass peak in addition to the one in Figure 3.

Zoom In Zoom Out Reset image size

4.2.1. Power Law

Figure 13 presents the marginalized distribution for the three power-law parameters M_min, M_max, and α (Section 3.3.1) from an analysis including the high-mass systems. The distribution for M_max is quite broad because the best-fit power laws slope downward (α < 0), making this parameter less relevant. The range −5.05 ⩽ α ⩽ −1.77 encloses 90% of the probability; the median value of α is −3.23. The presence of the high-mass samples in the analysis produces a distinctive tail, eliminating the correlations discussed in Section 3.3.1 and displayed in Figure 5 for the low-mass subset of the observations.

Zoom In Zoom Out Reset image size

4.2.2. Decaying Exponential

Figure 14 displays the marginalized distributions for the exponential parameters M_min and M₀ (Section 3.3.2) from an analysis including the high-mass systems. The distribution for the scale mass, M₀, has moved to higher masses relative to Figure 6 to fit the tail of the mass distribution; the distribution for M_min is less affected, though it has broadened somewhat toward low masses.

Zoom In Zoom Out Reset image size

4.2.3. Gaussian

Figure 15 displays the marginalized distributions for the Gaussian parameters (Section 3.3.3) when the high-mass objects are included in the mass distribution. The mean mass, μ, and the mass standard deviation, σ, are both increased relative to Figure 8 to account for the broader distribution and high-mass tail.

Zoom In Zoom Out Reset image size

4.2.4. Two Gaussian

The analysis of the two-Gaussian model shows the largest change when the high-mass samples are included. Figure 16 shows the marginalized distributions for the two-Gaussian parameters (Section 3.3.3) when the high-mass samples are included in the analysis. In stark contrast to Figure 9, there are two well-defined, separated peaks; the low-mass peak reproduces the results from the low-mass samples, while the high-mass peak (13.5534 ⩽ μ₂ ⩽ 27.9481 with 90% confidence; median 20.3839) matches the new high-mass samples. The peak in α near 0.8 is consistent with approximately four-fifths the total probability being concentrated in the 15 low-mass samples.

Zoom In Zoom Out Reset image size

4.2.5. Log Normal

The marginalized distributions for the log-normal parameters (Section 3.3.4) when the high-mass samples are included in the analysis are displayed in Figure 17. The changes when the high-mass samples are included (compare to Figure 10) are similar to the changes in the Gaussian distribution: the mean mass moves to higher masses and the distribution broadens. Because the log-normal distribution is inherently asymmetric, with a high-mass tail, it does not need to widen as much as the Gaussian distribution did.

Zoom In Zoom Out Reset image size

The confidence limits on the parameters for the parametric models of the underlying mass distribution are displayed in Table 5 (compare to Table 2).

Table 5. Quantiles of the Marginalized Distribution for Each of the Parameters in the Models Discussed in Section 3.3 when the High-mass Samples are Included in the Analysis (Compare to Table 2)

Model	Parameter	5%	15%	50%	85%	95%
Power law (Equation (7))	Mmin	4.87141	5.29031	5.85019	6.26118	6.45674
	Mmax	19.1097	23.4242	31.5726	37.7519	39.3369
	α	−5.04879	−4.30368	−3.23404	−2.31365	−1.77137
Exponential (Equation (11))	Mmin	4.0865	4.60236	5.32683	5.94097	6.22952
	M0	2.82924	3.41139	4.70034	6.52214	7.92979
Gaussian (Equation (13))	μ	7.86599	8.33118	9.20116	10.2493	10.9836
	σ	2.23643	2.58899	3.33545	4.17886	4.67881
Two Gaussian (Equation (15))	μ1	6.741	7.02724	7.48174	8.0139	8.46626
	μ2	13.5534	16.202	20.3839	24.9259	27.9481
	σ1	0.742824	0.913941	1.31244	1.94862	2.50238
	σ2	0.511159	1.5025	4.39824	7.04612	8.25905
	α	0.575692	0.670978	0.798227	0.891522	0.932143
Log normal (Equation (17))	〈M〉	8.00086	8.51192	9.6264	11.1851	12.3986
	σM	2.19262	2.8137	4.16742	6.25101	8.11839

Note. We indicate the 5%, 15%, 50% (median), 85%, and 95% quantiles.

Download table as:  ASCIITypeset image

4.2.6. Histogram Models

The non-parametric (histogram; see Section 3.4) models also show evidence of a long tail from the inclusion of the high-mass samples. Table 6 displays confidence limits on the histogram parameters for the analysis including the high-mass systems; compare to Table 3.

Table 6. The 5%, 15%, 50% (Median), 85%, and 95% Quantiles for the Bin Boundaries in the One-through Five-bin Histogram Models Discussed in Section 3.4 in an Analysis Including the High-mass, Wind-fed Systems

Bins	Boundary	5%	15%	50%	85%	95%
1	w0	2.22294	3.12695	4.2456	5.15132	5.58265
	w1	15.93	16.2535	17.7836	20.5449	22.5836
2	w0	3.87202	4.49983	5.41234	6.08334	6.35933
	w1	7.22163	8.25079	8.93669	9.71551	10.4287
	w2	18.4762	19.9798	24.941	32.5972	36.8615
3	w0	3.39289	4.24509	5.41694	6.15087	6.42822
	w1	6.41849	6.71984	7.47263	8.2942	8.61785
	w2	8.41449	8.64664	9.17056	10.4075	12.2718
	w3	18.5705	21.0481	27.1494	34.7753	38.0652
4	w0	2.42094	3.69875	5.2596	6.25449	6.54316
	w1	5.83725	6.2836	6.84987	7.8033	8.27706
	w2	6.94919	7.43628	8.38531	9.13401	9.91845
	w3	8.50371	8.75188	9.86694	17.1848	22.1086
	w4	18.5823	21.4628	28.367	35.8118	38.5278
5	w0	1.73691	3.19184	4.89769	5.9547	6.35522
	w1	5.46124	5.95881	6.59431	7.26795	7.91821
	w2	6.63468	6.9804	7.93239	8.60918	9.06926
	w3	7.89654	8.35634	8.91766	10.6568	13.9644
	w4	8.74064	9.42672	15.8004	22.7101	27.6399
	w5	20.0202	22.9065	29.6307	36.6606	38.8573

Note. The tails evident in Figure 12 are apparent here as well; compare to Table 3.

Download table as:  ASCIITypeset image

4.2.7. Model Selection for the Combined Sample

Repeating the model selection analysis discussed in Section 4.1.7 for the sample including the high-mass systems, we find that the model probabilities have changed with the inclusion of the extra five systems. As before, we assume for this analysis that the model priors are equal.

Reversible-jump MCMC calculations of the model probabilities are displayed in Figure 18; compare to Figure 11. The relative model probabilities are given in Table 7. The exponential model is the most favored model for the combined sample, with the two-Gaussian model the second-most favored. The ranking of models differs significantly from the low-mass samples. The improvement of the exponential model relative to the low-mass analysis is encouraging for theoretical calculations that attempt to model the entire population of X-ray binaries with this mass model. Note also that the increasing structure of the mass distribution favors histogram models with three bins over those with fewer bins.

Zoom In Zoom Out Reset image size

Table 7. Relative Probabilities of the Various Models from Section 3 Implied by the Combined Sample of Systems

Model	Relative Evidence
Exponential (Section 3.3.2)	0.346944
Two Gaussian (Section 3.3.3)	0.304923
Power law (Section 3.3.1)	0.120313
Log normal (Section 3.3.4)	0.102536
Histogram (3 Bin, Section 3.4)	0.0473464
Histogram (4 Bin, Section 3.4)	0.0282086
Histogram (2 Bin, Section 3.4)	0.0210994
Histogram (5 Bin, Section 3.4)	0.0179703
Gaussian (Section 3.3.3)	0.00901719
Histogram (1 Bin, Section 3.4)	0.00164214

Note. These probabilities have been computed from reversible-jump MCMC samples using the efficient jump proposal algorithm in Appendix B. (See also Figure 18.)

Download table as:  ASCIITypeset image

It is interesting to use our models for the underlying black hole mass distribution in X-ray binaries to place constraints on the minimum black hole mass implied by the present sample. Bailyn et al. (1998) addressed this question in the context of a "mass gap" between the most massive neutron stars and the least massive black holes. The more recent study of Ozel et al. (2010) also looked for a mass gap using a subset of the models and systems presented here. Both works found that the minimum black hole mass is significantly above the maximum neutron star mass (Kalogera & Baym 1996) of ∼3 M_☉ (though Ozel et al. 2010 only state their evidence for a gap in terms of the maximum-posterior parameters and not the full extent of their distributions).

The distributions of the minimum black hole mass from the analysis of the low-mass samples are displayed in Figure 19. The minimum black hole mass is defined as the 1% mass quantile, $M_{1\%}$, of the black hole mass distribution (i.e., the mass lying below 99% of the mass distribution). (A quantile-based definition is necessary in the case of those distributions that do not have a hard cutoff mass; even for those that do, like the power-law model, it can be useful to define a "soft" cutoff in the event that the lower-mass hard cutoff becomes an irrelevant parameter as discussed in Section 3.3.1.) For each mass distribution parameter sample from our MCMC, we can calculate the distribution's minimum black hole mass; the collection of these minimum black hole masses approximates the distribution of minimum black hole masses implied by the data in the context of that distribution. Figure 19 plots histograms of the minimum black hole mass samples.

Zoom In Zoom Out Reset image size

We find that the best-fit model for the low-mass systems (the power law) has $M_{1\%} > 4.3\,M_\odot$ in 90% of the MCMC samples (i.e., at 90% confidence). This is significantly above the maximum theoretically allowed neutron star mass, ∼3 M_☉ (e.g., Kalogera & Baym 1996). Hence, we conclude that the low-mass systems show strong evidence of a mass gap.

The distribution of minimum black hole masses for the analysis of the combined sample (i.e., including the high-mass systems) is shown in Figure 20. For the most favored model, the exponential, we find that $M_{1\%} > 4.5\,M_\odot$ with 90% confidence. We therefore conclude that there is strong evidence for a mass gap in the combined sample as well.

Zoom In Zoom Out Reset image size

Table 8 gives the 10%, 50% (median), and 90% quantiles for the minimum black hole mass implied by the low-mass sample; Table 9 gives the same, but for the combined sample of systems.

We have presented a Bayesian analysis of the mass distribution of stellar-mass black holes in X-ray binary systems. We considered separately a sample of 15 low-mass, Roche lobe filling systems and a sample of 20 systems containing the 15 low-mass systems and 5 high-mass, wind-fed X-ray binaries. We used MCMC methods to sample the posterior distributions of the parameters implied by the data for five parametric models and five non-parametric (histogram) models for the mass distribution. For both sets of samples, we used reversible-jump MCMCs (exploiting a new algorithm for efficient jump proposals in such calculations) to perform model selection on the suite of models. The consideration of a broad range of models and the model selection analysis, along with consideration of the full posterior distribution on the minimum black hole mass, significantly expand earlier statistical analyses of black hole mass measurements (Bailyn et al. 1998; Ozel et al. 2010).

For the low-mass systems, we found the limits on model parameters in Tables 2 and 3. The relative model probabilities from the model selection are given in Table 4. The most favored model for the low-mass systems is a power law. The equivalent limits on the model parameters for the combined systems are given in Tables 5 and 6. Unlike the low-mass systems, the most favored model for the combined sample is the exponential model. This difference indicates that the low-mass subsample is not consistent with being drawn from the distribution of the combined population.

We found strong evidence for a mass gap between the most massive neutron stars and the least massive black holes. For the low-mass systems, the most favored, power-law model gives a black hole mass distribution whose 1% quantile lies above 4.3 M_☉ with 90% confidence. For the combined sample of systems, the most favored, exponential model gives a black hole mass distribution whose 1% quantile lies above 4.5 M_☉ with 90% confidence. Although the study methodology was different, the existence of a mass gap was pointed out first by Bailyn et al. (1998) and most recently by Ozel et al. (2010) (who did not consider a power-law model, and applied both Gaussian and exponential models to the low-mass systems, where the exponential is strongly disfavored compared to our power-law model).

Theoretical expectations for the black hole mass distribution have been examined in Fryer & Kalogera (2001). They considered results of supernova explosion and fallback simulations (Fryer 1999) applied to single star populations; they also included a heuristic treatment of the possible effects of binary evolution on the black hole mass distribution. It is interesting that we find the most favored model for the combined sample to be an exponential, as discussed by Fryer & Kalogera (2001). On the other hand, we find the most favored model for the low-mass sample to be a power law, with the exponential model strongly disfavored for this sample. In agreement with Bailyn et al. (1998) and Ozel et al. (2010), we too conclude that both the low-mass and combined samples require the presence of a gap between 3 and 4–4.5 M_☉.

Fryer (1999) discussed two possible causes of such a gap: (1) a step-like dependence of supernova energy on progenitor mass or (2) selection biases. Current simulations of core collapse in massive stars may shed light on the dependence of supernova energy on progenitor mass. Selection biases can occur because the X-ray binaries with very low mass black holes systems are more likely to be persistently Roche lobe overflowing, preventing dynamical mass measurements. Ozel et al. (2010) conclude that the presence of such biases is not enough to account for the gap, arguing that the number (Equation (26)) of observed persistent X-ray sources not known to be neutron stars is insufficient to populate the 2–5 M_☉ region of any black hole mass distribution that rises toward low masses. Population synthesis models incorporating sophisticated treatment of binary evolution and transient behavior (e.g., Fragos et al. 2008, 2009) could help shed light on this possibility.

W.M.F., N.S., and V.K. are supported by NSF grants CAREER AST-0449558 and AST-0908930. A.C., L.K., and C.D.B. are supported by NSF grant NSF/AST-0707627. I.M. acknowledges support from the NSF AAPF under award AST-0901985. Calculations for this work were performed on the Northwestern Fugu cluster, which was partially funded by NSF MRI grant PHY-0619274. We thank Jonathan Gair for helpful discussions.

MCMC methods produce a Markov chain (or sequence) of parameter samples, $\lbrace \vec{\theta }_i \, | \, i = 1, \ldots \rbrace$, such that a particular parameter set, $\vec{\theta }$, appears in the sequence with a frequency equal to its probability according to a posterior, $p(\vec{\theta }|d)$. A Markov chain has the property that the transition probability from one element to the next, $p(\vec{\theta }_i \rightarrow \vec{\theta }_{i+1})$, depends only on the value of $\vec{\theta }_i$, not on any previous values in the chain.

One way to produce a sequence of MCMC samples is via the following algorithm, first proposed by Metropolis et al. (1953) and used widely in the physical sciences thereafter.

• 1.  
  Begin with the current sample, $\vec{\theta }_i$.
• 2.  
  Propose a new sample, $\vec{\theta }_p$, by drawing randomly from a "jump proposal distribution" with probability $Q(\vec{\theta }_i \rightarrow \vec{\theta }_p)$. Note that $Q(\vec{\theta }_i \rightarrow \vec{\theta }_p)$ can depend on the current parameters, $\vec{\theta }_i$, and any other "constant" data, but cannot examine the history of the chain beyond the most recent point. This is necessary to preserve the Markovian property of the chain.
• 3.  
  Compute the "acceptance" probability,

  $$\begin{equation} p_{\rm{accept}} \equiv \frac{p(\vec{\theta }_p|d)}{p(\vec{\theta }_i|d)} \frac{Q(\vec{\theta }_p \rightarrow \vec{\theta }_i)}{Q(\vec{\theta }_i \rightarrow \vec{\theta }_p)}. \end{equation} \tag{ A1 }$$
• 4.  
  With probability min (1, p_accept) "accept" the proposed $\vec{\theta }_p$, setting $\vec{\theta }_{i+1} = \vec{\theta }_p$; otherwise set $\vec{\theta }_{i+1} = \vec{\theta }_i$.

This algorithm is more likely to accept a proposed jump when it increases the posterior (the first factor in Equation (A1)) and when it is to a location in parameter space from which it is easy to return (the second factor in Equation (A1)); the combination of these influences in Equation (A1) ensures that the equilibrium distribution of the chain is $p(\vec{\theta }|d)$. As i → ∞ the samples $\vec{\theta }_i$ are distributed according to $p(\vec{\theta }|d)$.

In practice the number of samples required before the chain appropriately samples $p(\vec{\theta }|d)$ depends strongly on the jump proposal distribution; proposal distributions that often propose jumps toward or within regions of large $p(\vec{\theta }|d)$ can be very efficient, while poor proposal distributions can require prohibitively large numbers of samples before convergence.

There is no foolproof test for the convergence of a chain. In this work, we test the convergence of our chains in several ways. The most basic is by comparing the statistics calculated from the entire chain to statistics calculated from only the first half of the chain; when the chain has converged, the two calculations agree. This is a necessary, but not sufficient, condition for convergence.

We also have examined the sample traces from our chains, to see that the chains have densely and randomly sampled parameter space. A representative sample trace from our MCMC for the power-law model applied to the low-mass systems appears in Figure 21. Sample traces from MCMCs with other models are similar.

Zoom In Zoom Out Reset image size

Finally, for our most quantitative test of convergence we use the gibbsit code to implement the Raftery–Lewis convergence test for our quantile measures (Raftery & Lewis 1992a, 1992b, 1995). The most extreme quantile is the most difficult to determine accurately because—by design—there will be fewer samples in the tail than in the main body of a distribution obtained from an MCMC. Accordingly, we focus on the 90% confidence limit on the minimum black hole mass. For a quantile, q, the Raftery–Lewis test attempts to estimate how many samples from an MCMC are needed to determine q to within ±r at a confidence s. We use r = 0.0125 and s = 0.95. The Raftery–Lewis test approximates the MCMC chain as a two-state Markov chain, the two states being "within the quantile in question" and "outside the quantile in question." The 2 × 2 transition matrix for this two-state Markov chain and the associated uncertainty can be calculated analytically (Raftery & Lewis 1992a, 1992b, 1995), allowing the algorithm to determine the number of sample points required before the quantile of interest is determined sufficiently accurately.

For our chains, in the worst case (the power law on the lower-mass samples, as shown in Figures 4 and 21), we have twice as many samples as the Raftery–Lewis convergence test estimates we need to determine the 90% confidence level on the minimum mass; for all the other chains, we have about 20 times as many samples as the Raftery–Lewis criterion estimates are required. We suspect that the slow convergence of the power-law model on the lower-mass systems is due to the long tails in the mass parameters and the width of the distribution on the power-law exponent. In any case, the Raftery–Lewis test indicates that all our chains are converged sufficiently to determine the 90% quantile to within ∼1%.

We begin the chain at an arbitrary point in parameter space; this is equivalent to taking a finite section of an infinite chain that begins with the chosen point. Every point in parameter space occurs in an infinite chain, and no section of the chain is better than any other, so a sufficiently long, but finite, section of the infinite chain chosen in this manner can be representative of the statistics of the chain as a whole. However, because consecutive samples in a chain are correlated with each other, the beginning of our finite chain has a "memory" of the starting point; we discard enough points at the beginning of the finite chain that we can be confident it does not retain a memory of the arbitrary starting point. The points discarded in this way are commonly called "burn-in" points.