# A Systematic View of Ten New Black Hole Spins

Paul A. Draghis<EMAIL_ADDRESS>Department of Astronomy, University of
Michigan, 1085 South University Avenue, Ann Arbor, MI 48109, USA  Jon M.
Miller Department of Astronomy, University of Michigan, 1085 South University
Avenue, Ann Arbor, MI 48109, USA  Abderahmen Zoghbi Department of Astronomy,
University of Maryland, College Park, MD 20742, USA  HEASARC, Code 6601,
NASA/GSFC, Greenbelt, MD 20771, USA  CRESST II, NASA Goddard Space Flight
Center, Greenbelt, MD 20771, USA  Mark Reynolds Department of Astronomy,
University of Michigan, 1085 South University Avenue, Ann Arbor, MI 48109, USA
Department of Astronomy, Ohio State University, 140 W 18th Avenue, Columbus,
OH 43210, USA  Elisa Costantini SRON Netherlands Institute for Space Research,
Niels Bohrweg 4, 2333 CA Leiden, The Netherlands  Anton Pannekoek Astronomical
Institute, University of Amsterdam, P.O. Box 94249, 1090 GE Amsterdam, the
Netherlands  Luigi C. Gallo  Department of Astronomy & Physics, Saint Mary’s
University, 923 Robie Street, Halifax, Nova Scotia, B3H 3C3, Canada  John A.
Tomsick Space Sciences Laboratory, 7 Gauss Way, University of California,
Berkeley, CA 94720-7450, USA

###### Abstract

The launch of NuSTAR and the increasing number of binary black hole (BBH)
mergers detected through gravitational wave (GW) observations have
exponentially advanced our understanding of black holes. Despite the
simplicity owed to being fully described by their mass and angular momentum,
black holes have remained mysterious laboratories that probe the most extreme
environments in the Universe. While significant progress has been made in the
recent decade, the distribution of spin in black holes has not yet been
understood. In this work, we provide a systematic analysis of all known black
holes in X-ray binary systems (XB) that have previously been observed by
NuSTAR, but have not yet had a spin measurement using the “relativistic
reflection” method obtained from that data. By looking at all the available
archival NuSTAR data of these sources, we measure ten new black hole spins:
IGR J17454-2919 – $a=0.97^{+0.03}_{-0.17}$; GRS 1758-258 –
$a=0.991^{+0.007}_{-0.019}$; MAXI J1727-203 – $a=0.986^{+0.012}_{-0.159}$;
MAXI J0637-430 – $a=0.97\pm 0.02$; Swift J1753.5-0127 –
$a=0.997^{+0.001}_{-0.003}$; V4641 Sgr – $a=0.86^{+0.04}_{-0.06}$; 4U 1543-47
– $a=0.98^{+0.01}_{-0.02}$; 4U 1957+11 – $a=0.95^{+0.02}_{-0.04}$; H 1743-322
– $a=0.98^{+0.01}_{-0.02}$; MAXI J1820+070 – $a=0.988^{+0.006}_{-0.028}$ (all
uncertainties are at the $1\sigma$ confidence level). We discuss the
implications of our measurements on the entire distribution of stellar mass
black hole spins in XB, and we compare that with the spin distribution in BBH,
finding that the two distributions are clearly in disagreement. Additionally,
we discuss the implications of this work on our understanding of how the
“relativistic reflection” spin measurement technique works, and discuss
possible sources of systematic uncertainty that can bias our measurements.

##  1 Introduction

It is believed that black holes (BH) were first thought of almost 240 years
ago (Michell 1784), and then rediscovered 132 years later, in 1915, when Karl
Schwarzschild found a solution of general relativity that would describe a
black hole (Schwarzschild 1999). In 1972, 57 years later, multiple independent
studies concluded that the system Cygnus X-1 must harbor a black hole (Webster
& Murdin 1972; Bolton 1972; Tananbaum et al. 1972), making this the first BH
detected through electromagnetic radiation. In 2016, 44 years later,
gravitational wave signals from the merging event of two black holes was first
detected (Abbott et al. 2016), and only three years later, the first direct
image of the shadow of a Kerr black hole was obtained (Event Horizon Telescope
Collaboration et al. 2019). Even though our understanding of black holes has
grown exponentially in the past few decades and even though black holes are
intrinsically simple objects, described entirely by their mass and angular
momentum, many questions remain unanswered. This includes, perhaps most
interestingly, a complete understanding of the origin and magnitude of black
hole rotation.

The third Gravitational Wave Transient Catalog (GWTC-3, The LIGO Scientific
Collaboration et al. 2021a) contains 90 merging events detected by Advanced
LIGO (aLIGO) and Advanced Virgo (AdV). Based on this sample, The LIGO
Scientific Collaboration et al. (2021b) inferred that the spin distribution of
the most rapidly spinning BH of the two in the binary black hole (BBH) system
peaks near 0.4111the dimensionless spin parameter is given by $a=cJ/GM^{2}$,
where, in general, $-0.998\leq a\leq 0.998$ (Thorne 1974). However, for the
BHs in BBH, the absolute magnitude of the spin and the orientation of the spin
vector are reported., with the 1st and 99th percentiles at
$0.07^{+0.05}_{-0.03}$ and $0.80\pm 0.08$, while the spin distribution of the
less rapidly spinning BHs in BBH peaks around 0.2, with 99% of the
distribution below $0.54^{+0.09}_{-0.08}$. The spin distribution in these
systems is inconsistent with the spin distribution in X-ray binary systems at
more than 99.9% confidence (Fishbach & Kalogera 2022), hinting at different
formation mechanisms between the two populations of BHs. However, it is
important to mention that the BH spins obtained based on GW measurements are
not directly measured from the observation, rather they are calculated using
various models and theoretical assumptions based on the observable quantities
in the GW events.

Current estimates (Abbott et al. 2020a) predict that with KAGRA joining the
aLIGO/AdV collaboration and with the detectors reaching their designed
sensitivity, $79^{+89}_{-44}$ BBH mergers are expected to be detected during
O4, which is expected to begin in March 2023. In the near future, the
distribution of spin measurements both pre and post merger will become
increasingly statistically significant. In contrast, there have been $\sim 70$
XB BH candidates discovered since the beginning of the era of X-ray astronomy,
of which only $\sim 20$ have been confirmed dynamically (see the updated
results of Corral-Santana et al. 2016). In addition to the $\sim 1-3$ new BH
candidates discovered every year, this number is expected to increase
significantly in the future, in the era of 30 meter optical telescopes and
X-ray probes such as AXIS (Mushotzky et al. 2019). By looking at the
distribution of black hole spins in X-ray binaries (XB) and comparing it to
that of the spins in BBH mergers, we begin to look into the final destiny of
stars, the origin of black holes, and the nature of BBH systems. Fully
describing the two distributions of BH spins is essential to providing a
unified understanding of stellar-mass black hole formation and evolution. This
will help us to understand the entire BH population, whether or not the two
seemingly distinct BH populations have a common origin, and what factors might
lead to different median properties. It is therefore important to measure
spins in the greatest possible number of stellar-mass black holes in XB.

The spins of BHs in XB are most often measured through either continuum
fitting (see e.g., Gou et al. 2009; Steiner et al. 2012; Sreehari et al. 2020;
Zhao et al. 2021) or relativistic disk reflection (see e.g., Brenneman &
Reynolds 2006; Miller et al. 2009; Draghis et al. 2020b; Dong et al. 2022).
For a detailed review of BH spin measurement techniques, see Reynolds (2021).
Continuum fitting models the shape of the emission from the accretion disk
around a black hole, but requires prior knowledge of the BH mass, distance to
the system, accretion rate, and inclination of the inner accretion disk.
Additionally, continuum fitting measurements depend on the assumption of
spectral hardening factor and the nature of the hard emission (i.e. direct and
reflected coronal emission) in the observations.

Relativistic reflection (George & Fabian 1991) models the distortion of
spectral features (mainly the Fe K$\alpha$ spectral line, present at 6.4 keV
for neutral gas, and the Compton hump present above $\sim 20$ keV) in order to
infer the proximity of the emitting matter to the BH. The fluorescent Fe K
line is produced when one of the two K-shell electrons of an Fe atom is
ejected by an ionizing X-ray photon incoming from the compact corona.
Following the photoionization, an L-shell electron can decay to the K shell by
releasing 6.4 keV of energy either as a photon (with 34$\%$ probability) or as
an Auger electron (with 66$\%$ probability). As a Fe atom becomes more
ionized, the reduced number of screening electrons in the upper shells
increases the energy difference between the K and L shells, therefore causing
the fluorescent photons to be emitted with progressively higher energies, up
to 6.97 keV for H-like Fe XXVI, where due to the lack of electrons in the
K-shell, the Auger effect cannot occur and the Fe K line is produced through
recombination with a free electron. Unlike the soft incident X-ray photons
which are preferentially absorbed by the accretion disk, the hard X-ray
photons are preferentially Compton scattered, leading to a broad spectral
shape above $\sim 20~{}\rm keV$, referred to as the “Compton hump”.

By virtue of it being a relative measurement, this method infers distances
from the BH in units of gravitational radii ($r_{\rm g}=GM/c^{2}$), making it
an ideal tool for measuring the spin of BHs for which there are no mass or
distance estimates. However, while relativistic reflection does not require
prior information about the BH in the system, similarly to the continuum
fitting method, it also depends on the assumption that the optically thick,
geometrically thin accretion disk (Shakura & Sunyaev 1973) extends all the way
to the Innermost Stable Circular Orbit (ISCO) and that matter within the ISCO
is on plunging orbits and has no time to contribute significantly to the
observed flux. As the size of the ISCO is determined by gravity (Bardeen et
al. 1972; Novikov & Thorne 1973), relativistic reflection measures the size of
the BH ISCO, directly probing the BH spin. The assumption of an inner disk
edge consistent with the ISCO is motivated by numerical simulations (see e.g.,
Reynolds & Fabian 2008; Shafee et al. 2008; Schnittman et al. 2016) and
observational results (see e.g., Steiner et al. 2010; Salvesen et al. 2013;
García et al. 2015) suggesting that for Eddington fractions between
$10^{-3}\lesssim L_{disk}/L_{Edd}\lesssim 0.3$, there is a sharp inner disk
boundary, consistent with the ISCO. It is important to note that works such as
Tomsick et al. (2009) and Xu et al. (2020a) find evidence of disk truncation
at Eddington fractions of 0.14% and 0.18%, suggesting that while possible in
the low Eddington regime, it is not always a requirement for the accretion
disk to extend to the ISCO.

Additionally, relativistic reflection depends on assumptions regarding the
properties of the compact corona that produces the radiation incident of the
disk, that is then reprocessed (“reflected”). Most current models either
assume a “lamp-post” geometry, or parameterize the emissivity of the corona in
terms of radial emissivity indices, which remove the need of a geometrical
description of the corona. Lastly, the inclination of the inner disk regions
also has an effect on determining spectral features, and is usually treated as
a free parameter in reflection models. While some studies assume that the
inclination of the inner disk is the same as the orbital inclination, that is
not necessarily the case, and the two can be different due to disk tearing and
the Bardeen-Petterson effect (Bardeen & Petterson 1975; Nealon et al. 2015;
Liska et al. 2021).

The launch of NuSTAR (Harrison et al. 2013) in 2012 revolutionized
relativistic reflection spin measurements, with the size of the sample of
measured stellar-mass BH spins more than tripling through the addition of
measurements made using NuSTAR data. This improvement is attributable to the
wide band pass that NuSTAR offers, the high sensitivity and spectral
resolution at high energies, and the fact that the instrument does not suffer
from effects such as pile-up. The top panels of Figure 1 show the product of
the effective area and the live-time fraction of some major X-ray
observatories and the bottom panels show the typical number of counts as a
function of the number of resolution elements in the 5-7 keV (left) and 20-40
keV (right) bands, with the size of the points being proportional to the
signal to noise ratio (SNR). The number of counts was estimated for a source
with a flux of 1 Crab, accounting for the effective area of the instrument,
the live-time fraction, and an estimated typical duration of an observation
obtained through an analysis of archival observations. In the 5-7 keV band,
NuSTAR offers a combination of high SNR, large number of counts, and
sufficient number of resolution elements to accurately observe the shape of
the Fe K line. In the 20-40 keV, when compared to other instruments, NuSTAR
provides unparalleled capabilities to describe the Compton hump feature of
reflection spectra. All these advantages make NuSTAR the most appropriate
instrument for spin measurements through relativistic reflection.

There are (on the order of $\sim 40$) examples of individual black hole spin
measurements in papers adopting different methods. Additionally, most spin
measurements report only statistical uncertainties, while the systematic
uncertainties of the method have not yet been properly understood, due to the
effects mentioned above. This is particularly important for measurements such
as that of MAXI J1803-298, where Feng et al. (2021a) found $a=0.991\pm 0.001$,
or of EXO 1846-031, where Draghis et al. (2020b) measured
$a=0.997^{+0.001}_{-0.002}$, where even though the accuracy is probably very
high, the precision of the measurement is likely overestimated.

We aim to make a comprehensive, uniform treatment of black hole spin with
entirely consistent methods and systematic errors, which is needed in order to
completely understand the events that create black holes, the evolution of
black holes in binary systems, accretion onto black holes, the relativistic
effects onto the matter surrounding them, and the origins of the population of
massive black holes in binary mergers observed through gravitational waves.
Therefore, we analyzed the entire sample of archival NuSTAR observations of
X-ray binary systems harboring BHs that did not previously have a spin
measurement obtained through relativistic reflection from NuSTAR data. The
sources were selected initially based on the BHs identified in the WATCHDOG
(Tetarenko et al. 2016) and BackCAT (Corral-Santana et al. 2016) catalogs,
followed by a survey of journal literature indicating BH candidates. We
excluded sources that presented evidence or thermonuclear bursts (or burst-
like behavior) in order to exclude neutron star candidates. We analyzed 24
sources across 115 observations, and measured 10 new BH spins. We describe our
data analysis methods in Section 2, we present the analysis that led to the 10
new measurements in Section 3, and discuss the sources for which a spin
measurement was not possible in Appendix C. In Section 4, we present the
updated distribution of BH spins in XB. Lastly, we summarize our results and
discuss their implications in Section 5. For clarity, we grouped all figures
and tables in Appendix D.

##  2 Data Analysis

All the NuSTAR spectra were extracted using the nustardas v2.1.1 routines in
Heasoft v6.29c, using the NuSTAR CALDB version 20211103. The spectra from the
two NuSTAR FPM detectors were extracted from circular regions with 120” radius
centered at the location of the source, taken from the header of the event
file. Background spectra were extracted from annuli with an inner radius of
200” and an outer radius of 300”, concentric with the source region, and
subtracted from the source spectra. Our choice of background region does not
significantly change the background spectrum when compared to the spectrum
extracted from a circular background region selected on a different location
on the detector. This is especially the case since the observed sources are
generally bright, with count rates larger than a few counts per second. For
additional explanation regarding this choice of source and background regions,
together with the impact that this might have on the reported results, see
Appendix A. The spectra are then binned using the “optimal” binning scheme
presented by Kaastra & Bleeker (2016), using the “ftgrouppha” ftool, and
truncated at the energy where background emission begins to dominate over the
source spectrum.

The spectral analysis was performed using Xspec v12.12.0g (Arnaud 1996), using
$\chi^{2}$ statistics with standard weighting, and using the leven
minimization method. Throughout the analysis, we used the wilm abundances
(Wilms et al. 2000) and the vern photoelectric absorption cross-sections
(Verner et al. 1996).

When fitting NuSTAR spectra from the two FPM detectors, it is customary to
link all parameters of the model for the spectra from the two detectors and to
allow a multiplicative constant offset between the spectra. However, this
sometimes produces different residuals when fitting the two spectra,
especially at low energies $\leq 5\;\rm keV$ (see e.g., Draghis et al. 2022).
Therefore, instead of adopting a multiplicative constant to account for the
difference between the spectra from the two FPM NuSTAR detectors, we allowed
the normalizations of the model components to be different for the models
describing the two spectra, linking all other parameters. This method not only
accounts for differences in absolute flux sensitivity between the detectors,
but also for possible variations in spectral response. For example, instead of
fitting the two spectra with the model constant*TBabs*(diskbb+powerlaw) and
linking all parameters except for the constant component, we fit the two
spectra using the model TBabs*(diskbb+powerlaw), but allowing the
normalization of the diskbb and powerlaw components to vary between the two
detectors and linking all other parameters, therefore introducing one extra
free parameter. Despite having one extra free parameter, the quality of the
fits is almost always significantly improved by allowing free normalizations
when compared to having a multiplicative constant account for the difference
between spectra, and the parameters of the models are not generally strongly
influenced. By using this method, the fits can better constrain the shape of
the continuum without forcing the reflection features to account for
instrumental differences between the detectors. This allows more clearly
distinguishing reflection features from the underlying continuum.

###  2.1 Choice of Models

In order to test the possible systematic effects introduced through our choice
of models, throughout our analysis we start by fitting spectra with a baseline
array of seven models. First, all models include a TBabs component in order to
describe galactic absorption along the line of sight, using the value
estimated by the web version of the $\rm N_{H}$ Column Density
FTOOL222https://heasarc.gsfc.nasa.gov/cgi-bin/Tools/w3nh/w3nh.pl for each
source as a starting point. Second, in order to describe the emission from the
accretion disk, we include a diskbb component to all models, even those
describing observations that occurred during particularly hard states. If the
diskbb component is not required by the model, the temperature of the disk
and/or the normalization of the component will be driven to values indicating
that the component does not contribute significantly to the model. Lastly, to
describe the coronal emission, we add a powerlaw component. While this does
not probe the effects of relativistic reflection, it allows us to quantify
statistically if reflection is present: we can compare this to the quality of
the fit with models that account for reflection in order to probe if
reflection is indeed present. This gives us the first model variation, in
Xspec language TBabs*(diskbb+powerlaw).

In order to model reflection features, we used the relxill v.1.4.3 (Dauser et
al. 2014; García et al. 2014) family of models. For a complete description of
the models, their variations, and of all model parameters, see the website of
the model333http://www.sternwarte.uni-erlangen.de/~dauser/research/relxill/ ,
Section 3.1 in Draghis et al. 2021b, or Appendix A in Draghis et al. 2022. We
note that since the analysis was performed, an updated version of relxill was
launched (v2.0 - Dauser et al. 2022) that includes the effects of returning
radiation on the observed spectra. However, since the main focus of the paper
is to measure the spins of the compact objects in the systems and since the
addition of returning radiation does not influence spin measurements
significantly (see Figure 12b in Dauser et al. 2022), we chose to continue
using version 1.4.3 in our analysis.

To probe the extent of the theoretical assumptions of the relxill family of
models, we replace the powerlaw component with six variations of the relxill
components. In order to test the effect of the assumption regarding the shape
of the underlying continuum, we test both the simple relxill variant and
relxillCp. To test the effect of the assumed coronal geometry, we also test
the relxilllp variant, which assumes a “lamp-post” geometry, while the
previous two models attempt to measure the coronal emissivity directly,
without assuming a “lamp-post” geometry. Lastly, in order to probe the effect
of the density of the accretion disk, we also tested the relxillD variant.
However, when allowing the disk density to vary, it is almost always
unconstrained by fits. Therefore, we test 3 variations of the relxillD
component, by fixing the disk density at $n=10^{15}$, $10^{17}$, and
$10^{19}\;\rm cm^{-3}$. This gives us the array of six reflection components
that we used to describe relativistic reflection, by replacing the powerlaw
component in our first model variation. We leave all other parameters in the
models free, with the exception of the outer disk radius, which we fix at
$r_{\rm out}=990\;r_{g}$. As the radial emissivity profiles of the coronal
emission are generally steep, the outer regions of the disk contribute
insignificantly to the total observed emission, making the outer disk radius
nearly impossible to constrain for such systems. Therefore, we fixed this
parameter at nearly the maximum value allowed in the model. Additionally,
unless otherwise noted, we fix the inner disk radius at $r_{\rm in}=r_{\rm
ISCO}$ for observations that occur while the source was at an Eddington
fraction between $10^{-3}\lesssim L_{disk}/L_{Edd}\lesssim 0.3$, calculated
using the flux in the $0.1-10\;\rm keV$ band and the measurements in
literature for the BH mass and distance to the systems. Where no mass or
distance measurements were available, we analyze the observations that fall
within this Eddington fraction range when assuming generic values for these
parameters (e.g., $M_{\rm BH}\sim 8\pm 4M_{\odot}$ and $d\sim 10\pm 5\rm
kpc$).

###  2.2 MCMC analysis

After fitting the data with our baseline array of models and finding the best-
fit solution for each relxill variant tested, we ran a Markov Chain Monte
Carlo (MCMC) analysis on the models performing best in terms of the $\chi^{2}$
statistic. In cases where identifying the best performing models was not
trivial, we calculated the Akaike Information Criterion (AIC - Akaike 1974)
and the Bayesian Information Criterion (BIC - Schwarz 1978) using the
following formulas:

$AIC=Nln(\chi^{2}/N)+2N_{var}$ (1) $BIC=Nln(\chi^{2}/N)+ln(N)N_{var}\\\ $ (2)

where $N$ is the number of data bins and $N_{var}$ is the number of variables
in the model.

We used the “best-fit” parameter combination to generate a Gaussian proposal
distribution for initiating the walkers in the MCMC analysis, and used uniform
priors in the range allowed by the parameters in the model. For running the
MCMC analysis, we used the Xspec emcee implementation written by A. Zoghbi444
https://zoghbi-a.github.io/xspec_emcee/. We chose the number of walkers in the
analysis as a large integer, on the order of a few (4-5) times the number of
free parameters in the analysis. In oder to ensure chain convergence, we ran
the chains for $3-4000\times\tau_{f}$, where $\tau_{f}$ represents the
integrated autocorrelation time (Sokal 1996 recommend running the chains with
more steps than $\sim 1000\times\tau_{f}$ in order to ensure convergence).
Therefore, we ran the MCMC chains with 200 walkers for a total of $2.5\times
10^{6}$ steps, and disregarded the first $5\times 10^{5}$ steps of the chains
as a “burn-in” phase.

When running the MCMC analysis on more than one model, we also computed the
Deviance Information Criterion (DIC, Spiegelhalter et al. 2002) and reported
the model that performs best in terms of DIC. We computed the DIC using:

$DIC=\overline{D(\theta)}+p_{D}$ (3)

where $\theta$ represents the unknown parameters in the model, $D(\theta)$
represents the deviance, computed as $D(\theta)=-2\log(p(y|\theta))+C$ (where
$y$ represents the data, $p(y|\theta)$ is the likelihood function, and $C$ is
some constant which cancels out when comparing models), and $p_{D}$ represents
the effective number or parameters of the model, computed as
$p_{D}=\frac{1}{2}\overline{{\rm var}(D(\theta))}$ (Gelman et al. 2004). This
information criterion has the advantage of probing not only the “best-fit”
solution, but also the shape of the entire posterior distribution of the
parameter space around the best solution.

The MCMC analysis was run for each individual observation. For sources with
multiple observations, this produces a number of measurements for each
parameter in the models equal to the number of observations. We report the
modes of the posterior distributions and the $1\sigma$ highest posterior
density credible intervals of our measurements in the tables at the end of the
paper (see Section D). We chose to report the modes of the posterior
distributions instead of the medians in order to highlight the value that
occurs most often in the MCMC analysis for each parameter, suggesting a higher
preference for the specific value.

Additionally, throughout the paper, we chose to report the $1\sigma$ highest
posterior density credible intervals (i.e. the shortest interval in the
posterior distribution of a parameter that contains $68.3\%$ of the posterior
samples) instead of an equal-tailed confidence interval (i.e. having an equal
number of samples below the interval as above it) as this is better suited for
distributions that are clustered around the edge of a parameter space. As an
example, this is often the case for extremely high spin measurements.

###  2.3 Combining Measurements

When there are multiple measurements of the same parameter (i.e. BH spin or
inclination), it can be difficult to report a single value. While one could
choose to report the parameter combination from a single observation, such as
the one with highest total number of counts in the observation, the one
happening during the hardest spectral state, or the one where the reflection
strength is the highest relative to the continuum, this inevitably leads to a
loss of the information from the other observations. An alternative that takes
advantage of all available data is to fit all spectra simultaneously, and link
parameters that are expected to remain constant throughout the observations,
such as the BH spin and the inclination of the inner disk, ensuring that the
values measured are constrained by all the spectra at the same time. While
this method would be the ideal way forward to ensure that all data is fully
used, jointly fitting multiple NuSTAR spectra with complicated models that
include multiple components becomes a very slow, computationally intensive
process. In practice, fitting more than 2-3 observations, each having two
spectra, becomes nearly impossible, with each iteration in the $\chi^{2}$
minimization algorithm in Xspec lasting on the order of tens of seconds to
minutes. With the time required to run each iteration increasing significantly
through the addition of extra spectra, fitting $\sim 10$ observations
simultaneously (for sources such as 4U 1957+11, H 1743-322, or MAXI J1820+070)
is currently impossible due to limited computing power. Additionally, even if
one were to find the parameter combination for which the ”best fit” is
achieved, running the MCMC analysis described in Subsection 2.2 would be an
even more computationally demanding process, by many orders of magnitude.

Therefore, due to these limitations, we ran our analysis by fitting the
spectra produced by each observation individually. We did, however,
simultaneously fit the two spectra produced by the two FPM NuSTAR detectors
during the same observation, as described at the beginning of Section 2. We
implemented a Bayesian method of combining measurements from multiple
observations into a single value, that we report. While most parameters in the
models are not expected to stay constant in time, the BH spin, the inclination
of the inner accretion disk, and the Fe abundance are likely to not change in
time. However, due to possible correlations between the assumed disk density
in the model and the Fe abundance (see e.g., Tomsick et al. 2018), we do not
report a combined Fe abundance measurement, but we do report the individual Fe
abundance obtained from each observation in the tables in Section D for future
exploration of the models. Therefore, throughout the paper, for sources that
have been observed more than once, we report the spin and inclination obtained
by combining the measurements from individual observations, and weight the
contribution of each observation to the final measurement through the ratio of
the reflected to total flux in the 3-79 keV band. For a detailed explanation
of the method used to combine the results from multiple measurements and a
comparison between this method and a joint, simultaneous fit of multiple
observations, see Appendix B.

##  3 New Spins

###  3.1 IGR J17454-2919

At the time of the analysis, one NuSTAR observation of IGR J17454-2919 (ObsID
80001046002) was available in the archive. For this observation, we fit the
entire 3-79 keV spectrum. Fitting with TBabs*powerlaw gives
$\chi^{2}/\nu=544/433$. Reflection features and soft emission consistent with
an accretion disk are present. The spectra and residuals of the fits in terms
of $\sigma$ (i.e. $\frac{\rm data-model}{\rm error}$) are shown in Figure 2.

The best performing model, both in terms of $\chi^{2}$ and DIC, is
TBabs*(diskbb+relxillCp), with $\chi^{2}/\nu=445.48/421$. Based on the MCMC
analysis performed on this model, the measured spin is
$a=0.97^{+0.03}_{-0.17}$ and the inclination $\theta=60^{+8}_{-14}$ degrees.
The full set of parameters determined through the MCMC analysis is presented
in Table 1.

###  3.2 GRS 1758-258

Fitting the entire 3-79 keV NuSTAR spectrum of GRS 1758-258 from ObsID
30401030002 using TBabs*(diskbb+powerlaw) gives $\chi^{2}/\nu=700.72/483$,
with the residuals clearly indicating reflection features (see Figure 3).

The three best performing reflection models in terms of $\chi^{2}$ and DIC are
TBabs*(diskbb+relxillCp) ($\chi^{2}/\nu=442.56/474$, DIC = 467.31),
TBabs*(diskbb+relxillD) with Log(N)=19 ($\chi^{2}/\nu=445.95/475$, DIC =
478.26), and TBabs*(diskbb+relxill) ($\chi^{2}/\nu=449.65/474$, DIC = 475.2).
Of the three models, TBabs*(diskbb+relxillCp) performs best both in terms of
$\chi^{2}$ and DIC. Therefore, we report the measurements of this model, but
note that all measurements are consistent throughout the three models. Through
the MCMC analysis we find $a=0.991^{+0.007}_{-0.019}$ and
$\theta=66^{+8}_{-5}$ degrees. The entire set of parameters produced by the
MCMC analysis is presented in Table 1.

We note that Jana et al. (2022) published a spin measurement for GRS 1758-258
shortly before the submission of this paper. They measure $a>0.92$, and our
measurement is within good agreement. Additionally, we find a similar
ionization parameter, Fe abundance, and reflection fraction. However, the
inclination measurements and emissivity profiles are different, likely owing
to the fact that their analysis did not include a diskbb component, which they
argue was not statistically required by the data. We find that for all three
models presented above, the diskbb component is required, both in terms of
reduced $\chi^{2}$, AIC and BIC. Therefore we continue to include the diskbb
component in our analysis. The consistent spin measurements regardless of
choice of continuum is encouraging when considering the systematic effects
that act on these spin measurements. Since relativistic reflection works by
disentangling the reflection features from the underlying continuum, it is
encouraging that the models are able to isolate and characterize the shape of
reflection even when slightly different continua are chosen. Nevertheless, due
to the temporal proximity of the result of Jana et al. (2022) to the
submission of our manuscript, we chose to not exclude this source from this
sample.

###  3.3 MAXI J1727-203

MAXI J1727-203 was observed by NuSTAR once, in 2018 (ObsID 90401329002). We
fit the entire 3-79 keV spectra from the observation. Fitting with
TBabs*(diskbb+powerlaw) produces $\chi^{2}/\nu=567.22/451$, with the residuals
indicating clear reflection features. The spectra and residuals of the fits
are shown in Figure 4.

The two best performing models are TBabs*(diskbb+relxillCp)
($\chi^{2}/\nu=438.87/442$, DIC=688.31) and TBabs*(diskbb+relxill)
($\chi^{2}/\nu=440.59/442$, DIC=473.09). As it performs much better in terms
of DIC, we report the values measured from the MCMC run of the
TBabs*(diskbb+relxill) model. Through the MCMC analysis we find
$a=0.986^{+0.012}_{-0.159}$ and $\theta=64^{+10}_{-7}$ degrees, and the full
set of parameters is reported in Table 1.

###  3.4 Swift J1753.5-0127

There are 5 NuSTAR observations of Swift J1753.5-0127. Of those, only one
shows clear signs of reflection: ObsID 80001047002. Therefore, we use the
spectra extracted from the observation 80001047002 in the entire 3-79 keV
band. Fitting the spectra with TBabs*(diskbb+powerlaw) returns
$\chi^{2}/\nu=539.63/394$, with clear indications of relativistic reflection.
All reflection models with no coronal geometric assumption perform similarly,
while the lamp-post model appears to converge to a different, statistically
disfavored solution. The best performing reflection model, both in terms of
$\chi^{2}$ and DIC, is TBabs*(diskbb+relxillCp), with
$\chi^{2}/\nu=385.18/385$. The posterior MCMC samples are nearly entirely
concentrated at high spin values, measuring $a=0.997^{+0.001}_{-0.003}$ and
$\theta=74^{+3}_{-3}$ degrees. These measurements are consistent with the MCMC
run of other models without assumptions about the coronal geometry. We show
the spectra and residuals of the fits in Figure 5, and the modes of the
posterior distribution for each parameter in the MCMC analysis together with
their 1$\sigma$ credible intervals are presented in Table 1.

The MCMC run of the lamppost model provides similar inclination constraints,
but a much worse spin constraint $a=0.4^{+0.5}_{-0.8}$. This model produces a
fit worse by $\Delta\chi^{2}=18.7$ for 2 fewer free parameters. The relxillCp
flavor is preferred over the relxilllp both in terms of AIC (15.86 vs. 30.76)
and BIC (79.76 vs. 86.67). An F-test returns a probability of $10^{-4}$,
indicating that the non-lamppost models are indeed a statistically significant
improvement. Additionally, the relxilllp model performs worse than the
relxillCp model by $\Delta\rm DIC=13.57$. Lastly, visually inspecting the
residuals of the fit clearly shows that the Fe K line is not fit as well by
the relxilllp when compared to the relxillCp model. Therefore, we report the
results of our best performing model (TBabs*(diskbb+relxillCp)):
$a=0.997^{+0.001}_{-0.003}$ and $\theta=74^{+3}_{-3}$ degrees.

###  3.5 MAXI J0637-430

At the time of the analysis, there were 8 NuSTAR observations of MAXI
J0637-430, (ObsID 80502324002, 80502324004, 80502324006, 80502324008,
805023240010, 80502324012, 80502324014, 80502324016). The first three
(80502324002, 80502324004, and 80502324006) were taken during the soft state
of the late 2019 outburst of the source and are the only ones that show clear
reflection features. We fit the spectra from ObsID 80502324002 in the 3–50 keV
energy band, from ObsID 80502324004 in the 3–60 keV energy band, and from
ObsID 80502324006 in the 3–40 keV energy band, as the spectra are background
dominated at energies higher than the intervals mentioned. The unfolded
spectra of the three observations are shown in the top panels of Figure 6, and
the residuals of the fit to the TBabs*(diskbb+powerlaw) model are shown in the
middle panels of Figure 6, indicating signs of relativistic reflection.

When fitting the spectra with the models discussed in Subsection 2.1, the
relxill and relxillD with $\log(N)=19$ variants perform best in terms of
$\chi^{2}$ for all three observations. After running the MCMC analysis on the
two mentioned models for all three observations, the relxill variant performs
best in terms of DIC for ObsID 80502324004 and 80502324006, while the relxillD
with $\log(N)=19$ variant produces the better DIC for ObsID 80502324002. The
residuals of the fits using these models are shown in the bottom panels of
Figure 6, together with the fit statistic.

The modes of the posterior distributions for each parameter in the MCMC
analysis, together with their $1\sigma$ credible intervals are shown in Table
1. The histograms of the posterior distributions for spin and inner disk
inclination from the MCMC analysis on the three observations are shown in
Figure 7. Running the combining algorithm described in Subsection 2.3, we
obtain the distributions shown through the black curves in Figure 7, with
medians and credible intervals $a=0.97\pm 0.02$ and $\theta=62^{+3}_{-4}$
degrees.

###  3.6 V4641 Sgr

At the time of the analysis, four NuSTAR observations were public (ObsID
80002012002, 80002012004, 90102011002, and 90601302002). Using a mass of the
black hole of $M=6.4\pm 0.6\;\rm M_{\odot}$ and a distance to the system of
$d=6.2\pm 0.7\;\rm kpc$ (MacDonald et al. 2014), we calculate that only 3 of
the 4 observations happen while the source had a luminosity within the
Eddington range for which we expect the accretion disk to extend to the ISCO,
with ObsID 80002012004 falling outside of that range. Therefore, we continue
our analysis on the remaining 3 observations. We fit ObsID 80002012002 and
90102011002 in the 3–25 keV range, and ObsID 90601302002 in the 3–20 keV
range. The top panels in Figure 8 show the unfolded spectra from the three
observations, and the middle panels show the residuals of the fit to the
spectra when using TBabs*(diskbb+powerlaw). The fits are poor, and the models
clearly require additional components.

Simply replacing the powerlaw component with the six variations of the relxill
model that we test throughout our analysis improves the quality of the fit,
but not significantly. The quality of the fit is drastically improved through
the addition of an apec (Smith et al. 2001) component describing the emission
from a collisionally-ionized, diffuse gas in the vicinity of the source, which
is characterized by the plasma temperature, metal abundance, redshift, and
normalization. Additionally, we allow for the presence of ionized partial
covering of the diskbb component through a zxipcf (Reeves et al. 2008)
multiplicative component, which describes the effects of partially covering a
fraction _f_ of a source by a photoionized absorber. For all three
observations, the addition of the two new components improves the quality of
the fit significantly. The choice for the relxill component does not influence
the quality of the fit, so we continue using the default relxill variant,
making the complete model TBabs*(apec+zxipcf*diskbb+relxill). The residuals of
the fits are shown in the bottom panels of Figure 8, together with the fit
statistic.

We ran the MCMC analysis starting from the best-fit solutions for each of the
three observations. The modes of the posterior distributions for all the
parameters in the MCMC runs, together with their $1\sigma$ credible intervals
are presented in Table 2. Interestingly, in all 3 fits, the apec component
requires a blueshift $z\sim-0.03$, far above the terminal velocity of the
companion wind. This could indicate the presence of a disk wind, which has
been previously hinted at by Shaw et al. (2022). The quality of the fit
becomes significantly worse when attempting to fix the redshift of the
component to zero. While the Fe abundance $A_{\rm Fe}$ in the relxill
component is generally high, the abundance of metals in the apec component is
low, below unity.

Most importantly for this work, the spin of the compact object is generally
high and the inclination is broadly consistent between the three observations.
We combined the spin and inclination measurements from the three observations
as described in subsection 2.3. The results are shown in Figure 9. We measure
$a=0.86^{+0.04}_{-0.06}$ and $\theta=66^{+3}_{-4}$ degrees.

Interestingly, in Figure 9, we can see that ObsID 90601302002 does not
contribute strongly to the spin or inclination measurements, with both the
spin and inclination being determined by the other two observations analyzed.
While for inclination, the two observations produce agreeing results, the spin
measurement averages the two independent measurements. The uncertainty
reported on the parameters is simply statistical, but the difference between
the two posterior distributions that lead to the spin measurement of V4641 Sgr
highlight the importance of understanding the magnitude of the systematic
errors of spin measurements. A uniform treatment of the XB BH sample is
crucial for this goal.

The inclination of the inner accretion disk in V4641 Sgr is well determined by
our measurement, $\theta=66^{+3}_{-4}$ degrees, in good agreement with the
orbital inclination of the system (Orosz et al. 2001). Our measurement is
incompatible with a low inclination, suggesting that V4641 Sgr cannot be a
microblazar. However, Orosz et al. (2001) suggested that the jet angle must be
$\lesssim 10^{\circ}$. Such a discrepancy between the inclination determined
based on the jet angle and the inclination of the inner accretion disk was
found by Draghis et al. (2021b) in XTE J1908+094, where it was suggested that
the complexity of the local environment of the source together with an
inability to correctly “phase” different approaching/receding ejections can
alter our view of inclinations determined based on the radio jet.

###  3.7 4U 1543-47

At the time of the analysis, there were 10 public NuSTAR observations of 4U
1543-47 (ObsID 80702317002, 80702317004, 80702317006, 80702317008,
90702326002, 90702326004, 90702326006, 90702326008, 90702326010 and
90702326012), all taken during the 2021 outburst, with the source flux
decreasing in each new observation. However, assuming a BH mass of $M=9.4\pm
2\;\rm M_{\odot}$ and a distance to the system of $d=7.5\pm 1\;\rm kpc$ (Park
et al. 2004), all 10 observations occur while the source is at an Eddington
fraction larger than 0.3. Still, we analyze the last 4 observations taken
(ObsID 90702326006, 90702326008, 90702326010 and 90702326012), as they happen
while the source was at the lowest Eddington fraction of the 10 observations,
and also had the highest hardness of the 10 available observations. Within the
allowed range of BH mass and distance to the system, the Eddington fraction
during the last four observations varies between $\sim 0.32$ for the most
favorable parameter combination for the last observation (ObsID 90702326012)
and $\sim 0.8$ for the worst parameter combination for the first of the four
observations analyzed (ObsID 90702326006). It is important to note that the
first of the ten observations (ObsID 80702317002) gives an estimate between
$2-5$ for the Eddington fraction during the observation, suggesting that
perhaps our mass and distance estimates are not accurate. As the source was
outside of the Eddington ratio range for which we expect that the accretion
disk extends to the ISCO, we allowed the inner disk radius $r_{\rm in}$ to
vary during our spectral fitting. For all four observations we fit the spectra
over the entire NuSTAR pass band of 3–79 keV.

Figure 10 shows the unfolded spectra of the four observations in the top sub-
panels of each panel. The middle sub-panels show the residuals and statistics
produced when fitting the spectra with the simple TBabs*(diskbb+powerlaw)
model, clearly indicating the presence of relativistic reflection. Similarly
to the case of V4641 Sgr presented in subsection 3.6, simply replacing the
powerlaw component with any of the six variations of the relxill family that
we discussed in subsection 2.1 does not lead to a good fit.

We started fitting the observation with the hardest spectrum (ObsID
90702326008). The initial fit using TBabs*(diskbb+powerlaw) returns
$\chi^{2}/\nu=2074.52/431$, while replacing the powerlaw component with
relxill returns $\chi^{2}/\nu=623.42/422$ showing additional broad residuals.
To fit additional features, models are often expanded to include an additional
xillver component that has the role to account for distant, unblurred
reflection (see e.g., Miller et al. 2018). In this case, the addition of a
xillver component does not improve the quality of the fit, with the solution
simply converging to the same parameter combination as before by fitting the
reflection fraction of the xillver component to $R=0$.

Since the residuals indicate the presence of broad spectral features, we
replaced the xillver component with a second relxill component in the model,
with the role of mimicking a torn accretion disk. To do that, we linked the
black hole spin between the two components, the Fe abundance, the power law
indices $\Gamma$, the high-energy cutoff of the underlying power law
continuum, and the normalization of the two components. The reflection
fractions of the two components were allowed to vary, with the “inner” relxill
component taking positive values for $R$, while the “outer” component taking
only negative reflection fractions. Positive values of reflection fraction in
relxill components ensure that the model includes the direct coronal emission,
while negative reflection fraction values force the components to only include
the reflected emission. This way, the underlying continuum is included by the
“inner” relxill component. We allowed the inner disk radius of the “inner”
relxill component to vary, the outer disk radius of the “outer” relxill
component to $r_{\rm out}=1000\;\rm r_{g}$, and linked the outer radius of the
first relxill component to the inner radius of the second relxill component.
The disk inclination in the two relxill components was allowed to vary. This
model improves the quality of the fit, producing $\chi^{2}/\nu=540.68/416$.
The residuals of the fit clearly indicate a narrow absorption feature around
7.5 keV, and the measured inclination of the two relxill components is
different, but not well constrained: the inner component measures $\theta\sim
3^{\circ}$ while the outer component measures $\theta\sim 85^{\circ}$. When
accounting for the absorption feature using a gaussian additive component, and
also linking the inclination in the two relxill components, the fit improves
much further, producing $\chi^{2}/\nu=501.13/416$. With the inclination of the
two reflection components linked, the main difference between them is the
ionization, with the inner component taking large values $\log(\xi)\sim 4.6$,
while the outer one measures $\log(\xi)\sim 0$.

An alternative treatment that can be used to test for variable ionization
throughout the accretion disk is the relxilllpion model. We replaced the two
reflection components in our model with this new model, making the model
TBabs*(diskbb+relxilllpion). This also produces an improved fit when compared
to the the simple relxill or relxilllp variants, returning
$\chi^{2}/\nu=523.85/424$. Lastly, similarly to the case of V4641 Sgr, we
tested the addition of ionized partial covering (through zxipcf) and re-
emission (through apec). While the apec component did not improve the quality
of the fits in this case, the addition of the zxipcf component to the default
model improves the quality of the fit significantly
($\chi^{2}/\nu=508.48/418$) for the model TBabs*zxipcf*(diskbb+relxill). The
absorber requires moderate ionization ($\log\xi\sim 2$) and covering fraction
$f\sim 0.4$, consistent with a disk wind. The inner disk radius was allowed to
vary, however the fit converges to a value of $r_{\rm in}$ consistent with the
size of the ISCO.

We tested the three mentioned models on all four observations. In all four
cases, the last model (including zxipcf) was preferred in terms of AIC and
BIC, likely due to it producing a better statistic than the variable
ionization model and having a smaller number of free parameters than the model
with two reflection components. We mention that modifying the relxill variant
in these models does not strongly affect the quality of the fit or the
measured spin. Therefore we continue our analysis using the model
TBabs*zxipcf*(diskbb+relxill). The residuals produced by this model are shown
in the bottom panels of Figure 10. We ran the MCMC analysis on these four
observations with the given model. The modes and $1\sigma$ credible intervals
for all parameters in the model are presented in Table 2 for each observation.
The 1D posterior distributions of the spin and inclination parameters are
shown in Figure 11. By combining the posterior distributions as explained in
Section 2.3, we measure $a=0.98^{+0.01}_{-0.02}$ and $\theta=68^{+3}_{-4}$
degrees. Interestingly, this spin measurement disagrees with previous
measurements, namely $a=0.67^{+0.15}_{-0.08}$, measured through relativistic
reflection on RXTE data (Dong et al. 2020) and $a=0.80\pm 0.05$, measured
through continuum fitting (Shafee et al. 2006).

###  3.8 4U 1957+11

Using NuSTAR data, Sharma et al. (2021) estimated the black hole mass and
distance to 4U 1957+11 to be $M=5\pm 1\;\rm M_{\odot}$ and $d\sim 7\;\rm kpc$
and measured a spin of $a\sim 0.85$. Using Suzaku data, Nowak et al. (2012)
estimated a BH spin $a>0.90$ for $M=3\;\rm M_{\odot}$ and $d\sim 10\;\rm kpc$,
and Maitra et al. (2014) measured a spin $a>0.98$. All measurements were
performed using the continuum fitting method. 4U 1957+11 is a particularly
interesting source, due to its relatively low mass estimates and large
distance to the system, which makes this source consistently soft while at a
relatively low Eddington fraction. Regardless of the choice of black hole mass
and distance to the system of the two combinations, all 10 existing NuSTAR
observations are taken while the source was within the Eddington ratio limits
for which we expect the accretion disk to extend to the ISCO. We fit all
existing observations in the following energy bands: 30001015002 (3–20 keV),
30402011002 (3–20 keV), 30402011004 (3–50 keV), 30402011006 (3–50 keV),
30502007002 (3–45 keV), 30502007004 (3–50 keV), 30502007006 (3–60 keV),
30502007008 (3–60 keV), 30502007010 (3–25 keV), 30502007012 (3–20 keV). The
residuals when fitting the 10 spectra with TBabs*(diskbb+powerlaw) are shown
in the top panels in Figure 12, along with the statistic. The spectra of 4U
1957+11 are generally dominated by disk emission and the reflection features
are not always immediately obvious, but fitting with models that account for
relativistic reflection always improves the quality of the fit.

Replacing the powerlaw component in our model with the six variants of relxill
discussed in Subsection 2.1 produces better fits. We ran the MCMC analysis on
the two best performing models that account for reflection and selected the
one that produced the best DIC. The residuals of the best performing models in
terms of DIC are shown in the bottom panels of Figure 12, and the modes and
$\pm 1\sigma$ credible intervals of the posterior distributions for each
parameter are presented in Table 3. The 1D posterior distributions for the
spin and inclination are shown in Figure 13, with the transparency of the
lines corresponding to each observation being proportional to the weights used
when combining the measurements. We measure a spin of $a=0.95^{+0.02}_{-0.04}$
and an inclination of $\theta=52^{+4}_{-5}$ degrees. For this source in
particular, it is extremely important to weight the posterior distributions
from each individual measurement when combining them into a single
measurement, and the combined spin and inclination values are dominated by the
distributions inferred from the observations where reflection was strongest.
This spin measurement broadly agrees with the measurements made using
continuum fitting, and it could be used to place better constraints on the
mass of the black hole and the distance to the system.

###  3.9 H 1743-322

There were ten NuSTAR archival observations of H 1743-322 at the time of the
analysis. However, the spectra from ObsID 80002040004 are well fit by
TBabs*(diskbb+powerlaw) and the source is not detected in ObsID 80002040006.
Using a black hole mass of $M=11.21^{+1.65}_{-1.96}\;\rm M_{\odot}$ (Molla et
al. 2017) and a distance to the system of $d=8.5\pm 0.8\;\rm kpc$ (Steiner et
al. 2012), all remaining eight observations (80001044002, 80001044004,
80001044006, 80002040002, 80202012002, 80202012004, 80202012006, 90401335002)
fall within the $10^{-3}\lesssim L_{disk}/L_{Edd}\lesssim 0.3$ Eddington
fraction range. Therefore, we continue analyzing the eight NuSTAR observations
of H 1743-322, in the entire $3-79\;\rm keV$ NuSTAR band.

Using the TBabs*(diskbb+powerlaw) model to fit the spectra from the eight
NuSTAR observations indicates the presence of relativistic reflection. The
residuals of the fits are shown in the top panels of Figure 14, together with
the fit statistic. We accounted for reflection through the six models
discussed in Subsection 2.1, ran the MCMC analysis on the best performing two
models for each observation, and selected the best performing model of the two
(in terms of DIC) for each observation. The residuals of the fits using the
best reflection model for each observation are shown in the bottom panels in
Figure 14, together with the $\chi^{2}$ statistic produced by the models. The
modes and $1\sigma$ credible intervals of the posterior distributions for each
parameter produced by the reflection models returning the best DIC for each
observation are presented in Table 4, and the posterior distributions for the
black hole spin and inclination of the inner accretion disk are shown in
Figure 15. Combining the spin and inclination measurements produced by each
observation as described in Subsection 2.3 returns a spin of
$a=0.98^{+0.01}_{-0.02}$ and an inclination of $\theta=56\pm 4$ degrees.

Both the spin and inclination measurements of our analysis disagree with those
determined by Steiner et al. 2012 through continuum fitting ($a=0.2\pm 0.3$)
and from radio and X-ray jets ($\theta=75\pm 3$). While it is expected that
the jet axis coincides with the black hole spin axis and while that often is
the case in observations (see e.g., Cyg X-1 - Krawczynski et al. 2022, or MAXI
J1820+070 in Subsection 3.10), disparities between the inclination inferred
through jet morphology and through relativistic reflection have previously
been found (see e.g., XTE J1908+094 - Draghis et al. 2021b). Interestingly,
Tursunov & Kološ 2018 estimated that the black hole in H 1743-322 has a spin
of $a=0.4\pm 0.2$ from quasiperiodic oscillations (QPO). However, also from
QPOs, Mondal 2010 find that the black hole in H 1743-322 must have a spin
$a>0.68$. It is important to note that there is no commonly accepted model to
predict BH masses and spins based on QPOs, and models often disagree (Stefanov
& Tasheva 2017). Furthermore, as QPO frequencies often shift, it is likely
that they do not properly trace the BH ISCO.

###  3.10 MAXI J1820+070

NuSTAR observed MAXI J1820+070 27 times in 2018 and 2019. Using a black hole
mass of $M=8.48\pm 0.79\;\rm M_{\odot}$ (Torres et al. 2020) and a distance of
$d=2.96\pm 0.33\;\rm kpc$ (Atri et al. 2020), we find that 19 of the 27
observations happened while the source was at an Eddington fraction between
$10^{-3}\lesssim L_{disk}/L_{Edd}\lesssim 0.3$: ObsID 90401309002,
90401309004, 90401309006, 90401309008, 90401309010, 90401309012, 90401309013,
90401309014, 90401309016, 90401309018, 90401309019, 90401309021, 90401309023,
90401309025, 90401309026, 90401309027, 90401309033, 90401309035, and
90401324002. We ran our analysis on these 19 observations, fitting in the
entire $3-79\;\rm keV$ NuSTAR band pass in each case. The residuals produced
by fitting the spectra from each observation with the model
TBabs*(diskbb+powerlaw), and the fit statistic are shown in the top panels in
Figures 16 and 17.

The residuals generally indicate the presence of relativistic reflection, so
we fit the spectra from the 19 observations with the six models discussed in
Subsection 2.1. Due to the large number of observations, we only ran the MCMC
analysis on the best performing models for each observation in terms of
$\chi^{2}_{r}=\chi^{2}/\nu$, where $\nu$ represents the number of degrees of
freedom in the model. The bottom panels in Figure 16 and continued in Figure
17 show the residuals of the best performing reflection models. We ran the
MCMC analysis on the fits to the spectra from the 19 observations. Table 5
shows the modes of the posterior distributions for each parameter in the MCMC
analysis for the 19 observations, along with the best performing relxill
variant, that we used in our analysis. The posterior distributions for spin
and inclination are shown in Figure 18. When combining the posterior
distributions while weighting them by the ratio of reflected to total flux in
the $3-79$ keV band, we obtain a spin of $a=0.988^{+0.006}_{-0.028}$ and an
inclination of $\theta=64^{+3}_{-4}$ degrees.

Our inclination measurement is in good agreement with that determined by Atri
et al. (2020) from the radio jet ($\theta=63\pm 3^{\circ}$). However, our spin
measurement disagrees with the values found using continuum fitting:
$a=0.14\pm 0.09$ (Zhao et al. 2021) and $a=0.2^{+0.2}_{-0.3}$ (Guan et al.
2021). However, Bhargava et al. (2021) used a Relativistic Precession Model
for QPO frequencies on NICER data to measure a spin of
$a=0.799^{+0.016}_{-0.015}$, and Prabhakar et al. (2022) analyzed multiple
Swift/XRT (Burrows et al. 2005), NICER (Gendreau et al. 2012), NuSTAR, and
AstroSat (Agrawal 2006) observations and found that when fitting the
reflection spectra, the spin was always high, so they fixed it at the maximum
allowed value throughout their analysis. While Guan et al. (2021) argue that
given their low spin measurement, the strong jet seen in MAXI J1820+070 must
be powered by the accretion disk, our measurement allows for the jet to be
powered by the high back hole spin through the Blandford–Znajek mechanism
(Blandford & Znajek 1977). This would be expected, especially given the good
agreement between the inclination of the black hole spin axis measured from
relativistic reflection in this work and the inclination of the observed radio
jet.

##  4 Spin Distribution

The new era of gravitational wave (GW) signals from binary black hole (BBH)
mergers holds tremendous promise. Currently, two challenges keep the full
potential of black hole mergers from being realized: (1) a degeneracy between
the mass ratio in the binary and a combination of the spins of the black
holes, and (2) a second degeneracy between the two spins themselves, which
makes it difficult to measure individual spins (Pürrer et al. 2016). It was
estimated that the degeneracies can be broken for signal to noise ratios (SNR)
of 100, but only a few such observations are likely to occur every year
(Pürrer et al. 2016). Moreover, even in this situation, the lowest mass black
hole in the system can only be classified as having ”high” or ”low” spin, as
opposed to an actual measurement. Given current event rates, the most
pragmatic path forward is to obtain informative priors, based on a robust
derivation of the spins of stellar-mass black holes in X-ray binary systems:
the posterior distributions for GW spin measurements are correlated to the
assumed prior distribution (Vitale et al. 2017; Abbott et al. 2019, 2020b),
making it crucial to have educated predictions for the prior distributions.

We compiled all XB BH spin measurements in literature obtained through
continuum fitting and relativistic reflection in Table 6 and plotted them on
the same scale in Figure 19, with the orange points showing measurements made
through continuum fitting, yellow points representing relativistic reflection
measurements made on data from instruments other than NuSTAR, blue points
represent reflection measurements obtained using NuSTAR data, and the green
points showing the 10 measurements of this paper.

A few interesting points can be made regarding the existing sample of BH spins
in XB. Firstly, it is worth nothing that there are no precise negative BH spin
measurements, only upper limit measurements suggesting $a\leq 0$. Of the seven
measurements consistent with negative spins, five are made using continuum
fitting and two are made using relativistic reflection, but not using NuSTAR
data. Also, only four of the seven measurements that allow negative spins
exclusively require negative spins. A few of those constraints are
particularly interesting. For GX 339-4, the continuum fitting method finds
$a<0.9$ (Kolehmainen & Done 2010), while relativistic reflection on NuSTAR
data finds $a=0.95_{-0.08}^{+0.02}$ (Parker et al. 2016). While the continuum
fitting measurement nominally allows low and negative spins, simultaneously
considering the two measurements made using the two methods indicates that the
BH spin in GX339-4 is indeed high. Another particularly interesting source is
IGR J17091-3624, for which continuum fitting measures $a<0.2$ (Rao & Vadawale
2012) while relativistic reflection on NuSTAR data finds $a=0.07\pm 0.20$
(Wang et al. 2018). Again, while considering only the continuum fitting
measurement negative spins are allowed, treating the two measurements together
indicates the BH is likely slowly rotating, with spin consistent with zero.

Secondly, there appears to be a trend of decreasing measurement uncertainty
with increasing BH spin. This trend applies both to spins obtained through the
reflection method and through continuum fitting. As the size of the ISCO is
set by the BH spin, a truncated disk around a rapidly spinning BH would
produce similar spectral features as a disk that extends all the way to the
ISCO around a slowly rotating BH (e.g., a maximally spinning BH with a disk
truncated at $6\;r_{g}$ would produce similar features as a non-spinning BH
with a disk that extends to the ISCO). However, since XB spin measurements
assume inner disk radii consistent with the ISCO and since possible disk
truncation would have the effect of biasing spin measurements to lower values,
the high measured spins cannot be systematically mistaken and are likely
accurate. The increased precision of high measured spins could be owed to
higher both direct and reflected emissivity with increased proximity to the
BH, making the spectral features easier to observe and characterize, leading
to smaller uncertainties on spin measurements. Works such as Bonson & Gallo
(2016) and Kammoun et al. (2018) test the accuracy of spin measurements in
AGN, finding that for values above 0.8, the spin is better constrained both in
terms of accuracy and precision.

Thirdly, it seems like the overall distribution of spin measurements is
concentrated around high values. Using the same method that we used to combine
multiple measurements into a single distribution, we combined the multiple
independent spin measurements in Table 6 and Figure 19 in an attempt to
understand the entire observed BH spin distribution in XB. It is crucial to
acknowledge that this distribution is likely influenced by observational and
selection biases, which can be better understood through a uniform treatment
of the entire spin sample. Additionally, it is important to note that the
measurements in literature vary widely in terms of the energy range and
resolution of the instrument used and of theoretical assumptions and numerical
resolution of the models used, which can likely lead to biases in the
measurements.

We plot the combined distributions in Figure 20. The orange curve shows the
distribution of XB BH spin measurements obtained through continuum fitting.
The green curve shows the distribution inferred based on the 10 measurements
presented in this paper, which is a subset of the blue curve, showing the
distribution inferred based on all measurements made using relativistic
reflection on NuSTAR data. The blue curve is a subset of the yellow curve,
which shows the distribution of all XB spin measurements obtained through
relativistic reflection. The black curve shows the combined distribution of
all XB spin measurements, through both relativistic reflection (yellow) and
continuum fitting (orange). For reference, we compare this with the spin
distribution inferred by The LIGO Scientific Collaboration et al. (2021b)
based on all the GW signals from BBH mergers presented in GWTC-3, shown in
purple.

The spin distributions in BBH and in XB are clearly in disagreement. Fishbach
& Kalogera (2022) compared the distribution of spins in the previous version
of the GWTC (GWTC-2, Abbott et al. 2021) to the distribution of BH spins in
high mass X-ray binaries (HMXB), inferred from 3 sources (LMC X-1, Cyg X-1,
and M33 X-7). Despite the very limited sample size and despite only using the
values measured through continuum fitting, Fishbach & Kalogera (2022) conclude
that the HMXB and BBH spin distributions disagree at the $\geq 99.9\%$ level.
However, they conclude that the spin distribution of BHs in low Mass X-ray
binaries (LMXB) agrees with that of HMXB. There is value in attempting to
compare the BBH sample to HMXB only, since HMXB are candidates to evolving to
produce BBH systems. However, it is worth mentioning that studies have shown
that neither Cyg X-1 (Belczynski et al. 2011) nor M33 X-7 (Ramachandran et al.
2022) are likely to evolve to produce a BBH system. Additionally, works such
as Gallegos-Garcia et al. (2022) find that at most 20% of BBH originate from
HMXB systems and at most 11% of HMXB evolve into BBH systems. Here, we
highlight the comparison of all BH spins in XB obtained through continuum
fitting and relativistic reflection to the spins in BBH. We performed a two-
sample Kolmogorov-Smirnov test (K-S test - Smirnov 1939) in order to quantify
the probability that the two observed distributions of spins in BBH and XB are
drawn from the same unknown underlying distribution. Using the distributions
shown in Figure 20, we generated $10^{5}$ samples of 50 “measurements”
randomly drawn from the XB population and 200 “measurements” drawn from the
BBH population, simulating the existing sample sizes. Then, we performed the
two-sample K-S test and analyzed the distribution of probabilities produced.
When comparing the BBH distribution with all spin measurements in XB, we find
that the highest probability across the $10^{5}$ tests is $p\sim 2\times
10^{-5}$, while when comparing it with the distribution of spins obtained
using relativistic reflection, the highest probability was $p\sim 5\times
10^{-13}$. The median probability across the $10^{5}$ samples in the two cases
agrees, $p\sim 9\times 10^{-16}$. These results clearly state that the
distribution of spins in BBH is different from the observed distribution of BH
spins in XB. Since the formation channels that lead to the different spin
distributions are likely different, it is important to consider the need to
explain both distributions when attempting to provide a unified formation
mechanism for all black holes.

##  5 Discussion

In this work, we analyzed the archival NuSTAR observations of a sample of 24
X-ray binary sources containing BHs that do not have previous spin constraints
using relativistic reflection on NuSTAR data. We measured 10 new BH spins. In
14 cases, we were unable to obtain a BH spin measurement. The 10 new
measurements represent an increase in the sample size of BH spin measurements
obtained through the relativistic reflection method using NuSTAR data of
nearly 40%, and an increase of nearly 25% in the sample size of all
relativistic reflection measurements. Four of the ten sources wherein we
measure the spin using reflection have previous spin constraints from
continuum fitting. Interestingly, in three of the four sources, our
relativistic reflection measurements disagree with the spins obtained through
continuum fitting (MAXI J1820+070, H 1743-322, and 4U 1543-47), while for 4U
1957+11 our measurement is consistent within the uncertainties with that
obtained through continuum fitting.

With the exception of V4641 Sgr, the new spin measurements have very high
values, $a>0.95$, with all values being larger than 0.8 within $1\sigma$. This
is consistent with the spin distribution in X-ray binaries (Fishbach &
Kalogera 2022). It is unlikely that the relativistic reflection method is
biased to measure high black hole spins, as low spins produce spectral shapes
that are less blurred, and therefore easier to distinguish from the underlying
continuum. Additionally, low and moderate spins are still detected by works
such as Wang et al. (2018); Draghis et al. (2021b); Jana et al. (2021a) or Jia
et al. (2022). However, one cannot exclude the possibility that there is an
observational bias, making highly spinning BHs easier to detect. Vasudevan et
al. 2016 explain that due to selection effects of flux-limited observations,
the BH spin distribution measured in Active Galactic Nuclei (AGN) is biased
toward high values, but they also argue that this effect should not be
directly applicable to stellar mass galactic BHs. Hirai & Mandel (2021)
explain that high mass X-ray binaries (HMXB) should only become observable if
the companion star has a Roche lobe filling factor above 0.8-0.9. A BH
rotating in retrograde direction or slowly in the prograde direction will have
a large ISCO radius; when paired with a high Roche lobe filling factor of the
companion, this might not allow enough room for the formation of an accretion
disk around the BH, with matter plunging directly towards the compact object.
If an accretion disk cannot be formed, the X-ray binary will not be observed,
biasing our observations toward high spin values.

At the same time, one cannot exclude the possibility that the intrinsic spin
distribution of BHs in X-ray binaries strongly favors high spins. Due to the
short lifetimes of HMXBs (van den Heuvel 1976; Iben et al. 1995), we know that
the observed high spins (see e.g., Parker et al. 2015; Jana et al. 2021a) must
be natal. Works such as Fragos & McClintock (2015) argue that the high spins
in galactic low mass X-ray binaries (LMXB) can be explained through long
episodes of highly efficient accretion. However, it is important to
acknowledge that a BH must accrete 80$\%$ of its initial mass to increase its
spin from $a=0$ to $a=0.9$ and 122$\%$ of its initial mass to increase its
spin from $a=0$ to $a=0.99$ (Bardeen 1970; Belczynski et al. 2008), so a
limited reservoir of matter might prevent significant spin changes even in
LMXB. Additionally, works such as Draghis et al. (2022) find that even the BHs
in intermediate or low mass X-ray binaries can be born with high spins, not
requiring significant accretion.

We measured ten high BH spins, with nine of them approaching the maximal
possible value. However, only three of the ten sources show resolved radio
jets or ”cannonball” ejections: MAXI J1820+070 (Homan et al. 2020), H 1743-322
(Hao & Zhang 2009), and GRS 1758-258 (Martí et al. 2002). This suggests that a
high BH spin is not sufficient for production of resolvable radio jets, and
that the observability of the ejecta is likely more strongly influenced by
system properties such as the mass accretion rate in the system, the orbital
period of the binary, the surrounding interstellar medium (ISM), and even the
distance to the system.

The reported uncertainties of our measurements are purely statistical, and the
systematic uncertainties are likely to be smaller, yet not negligible. As the
BH spin and inclination of the inner disk influence spectral features
similarly, we show all the individual spin and inclination measurements for
all observations treated in this work in Figure 21, aiming to highlight any
possible correlation between the two parameters in the models. While there
appears to be a positive trend, that is largely owed to the measurements of
MAXI J1820+070.

In Figure 22, we show the spin and inclination measurements obtained from the
19 observations of MAXI J1820+070, with the transparency of the points being
proportional to the ratio of reflected to total flux in the observation, which
was used to weight the measurements when computing the combined measurement.
The combined measurements for spin and inclination are shown through solid
lines, while their uncertainties are presented through the shaded regions.
Previously, any of those measurements could have been reported in literature,
but by looking at all the observations we obtain a value that is not biased by
observation-specific peculiarities and more clearly reflects the true value of
the parameters. MAXI J1820+070 was one of the brightest X-ray binary systems
ever observed, and with it being extensively observed in a comprehensive
manner across the electromagnetic spectrum, the existing data set makes this
source an anchor in our understanding of electromagnetic emission from stellar
mass BHs.

Figure 21 also shows that few spin measurements take low or negative values,
most originating from observations of 4U 1957+11. In Figure 23, we show the
inclination vs spin parameter space for the 10 observations of 4U 1957+11,
with the transparency of the points being proportional to the ratio of
reflected to total flux. Similar to Figure 22, the solid lines and shaded
regions represent the values and uncertainties obtained while combining the
individual measurements. While there are some low and negative spin
measurements, those come from observations in which reflection is not as
strong, with the combined measurements being dominated by the observations in
which reflection is more strongly present. Figures 22 and 23 highlight the
importance of analyzing all available observations with the aim of reducing
any possible systematic uncertainties of our measurements.

Having such a large data set produced through consistent models and
assumptions and with uniform systematic uncertainties also allows us to look
into the systematic effects of parameters in the models on the spin
measurements. In Figures 24 and 25, we show all the individual spin
measurements presented in this paper, in relation with the Fe abundance
($A_{\rm Fe}$) and reflection fraction ($R$) predicted by the model,
respectively. While it appears that the models prefer high Fe abundance and
low reflection fraction, there is no apparent direct correlation between these
values and the spin measurement. Additionally, we can begin to look into the
general behavior of the models by highlighting correlations between other
parameters. For example, in Figure 26 we show the Fe abundance in all
measurements in comparison with the reflection fraction. In this case, no
direct correlation is present. However, such parameter correlations can also
be treated at the level of each individual fit through the posterior
distributions obtained through the MCMC runs. While inspecting the corner
plots showing trends in the posterior distributions of parameters, no obvious
correlations affecting the spin were found. Still, on occasion, degeneracies
between parameters were present. Such examples include the correlation between
the outer emissivity index $q_{2}$ and the breaking radius $R_{\rm br}$
parametrizing the coronal emissivity, which often lead to unphysically low
values of $q_{2}$ and high values of $R_{\rm br}$ due to the inability of the
data to constrain the parameters. Due to the large number of parameters, the
scale of a complete treatment of the parameter space is substantial and beyond
the scope of this work and will be treated in a future paper, in which we also
plan to expand the sample size of spin measurements with existing measurements
reanalyzed using the same assumptions of this paper.

While based on the entire GWTC-3 sample, The LIGO Scientific Collaboration et
al. (2021b) concluded that most BHs in BBH systems are slowly rotating.
However, a few rapidly rotating BHs were identified in this sample, such as
GW190412 (either the spin of the more massive component is
$a_{1}=0.44^{+0.16}_{-0.22}$ \- Abbott et al. 2020c, or the spin of the less
massive component is $a_{2}=0.88^{+0.11}_{-0.24}$ \- Mandel & Fragos 2020),
GW190517_055101 ($a_{1}=0.86^{+0.13}_{-0.35}$ \- Abbott et al. 2021), GW190403
051519 ($a_{1}=0.92^{+0.07}_{-0.22}$ \- The LIGO Scientific Collaboration et
al. 2021c), or GW191109_010717 ($a_{1}=0.83^{+0.15}_{-0.58}$ \- The LIGO
Scientific Collaboration et al. 2021a). The formation of BBH systems is a
topic of debate, and the current preferred view is that the entire observed
distribution cannot be explained through a single formation channel (see e.g.,
Zevin et al. 2021 for a detailed view on BBH formation channels). Fishbach et
al. (2022) find that based on the existing sample of GW measurements, at most
26% of the underlying BBH population originates from ”hierarchical mergers”
(repeated mergers of smaller BHs). While works such as Olejak & Belczynski
(2021) argue that the observed distribution of BBH is consistent with
formation from isolated binaries, Gallegos-Garcia et al. (2022) find that at
most 20% of BBH originate from HMXB systems and at most 11% of HMXB systems
evolve into BBH, suggesting that the spin distributions in BBH and XB are not
only observed to be different (Fishbach & Kalogera 2022), but intrinsically
different. In the future, expanding and comparing the two spin samples will
allow placing better constraints on BH formation and evolution mechanisms.

AZ is supported by NASA under award number 80GSFC21M0002. JAT acknowledges
partial support from NASA through NuSTAR Guest Observer grants 80NSSC20K0644
and 80NSSC22K0059.
Software: Astropy (Astropy Collaboration et al. 2013, 2018), emcee (Foreman-
Mackey et al. 2013), numpy (Harris et al. 2020), matplotlib (Hunter 2007),
scipy (Virtanen et al. 2020), pandas (Reback et al. 2022; Wes McKinney 2010),
corner (Foreman-Mackey 2016), iPython (Pérez & Granger 2007), Xspec (Arnaud
1996), relxill (Dauser et al. 2014; García et al. 2014).

## References

  * Abbott et al. (2020a) Abbott, B. P., Kagra Collaboration, Ligo Scientific Collaboration, & VIRGO Collaboration. 2020a, Living Reviews in Relativity, 23, 3, doi: 10.1007/s41114-020-00026-9
  * Abbott et al. (2019) Abbott, B. P., LIGO Scientific Collaboration, & Virgo Collaboration. 2019, Physical Review X, 9, 031040, doi: 10.1103/PhysRevX.9.031040
  * Abbott et al. (2016) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2016, Phys. Rev. Lett., 116, 061102, doi: 10.1103/PhysRevLett.116.061102
  * Abbott et al. (2020b) Abbott, R., LIGO Scientific Collaboration, & Virgo Collaboration. 2020b, arXiv e-prints, arXiv:2010.14527.  https://arxiv.org/abs/2010.14527
  * Abbott et al. (2020c) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020c, Phys. Rev. D, 102, 043015, doi: 10.1103/PhysRevD.102.043015
  * Abbott et al. (2021) —. 2021, Physical Review X, 11, 021053, doi: 10.1103/PhysRevX.11.021053
  * Agrawal (2006) Agrawal, P. C. 2006, Advances in Space Research, 38, 2989, doi: 10.1016/j.asr.2006.03.038
  * Akaike (1974) Akaike, H. 1974, IEEE Transactions on Automatic Control, 19, 716
  * Arnaud (1996) Arnaud, K. A. 1996, in Astronomical Society of the Pacific Conference Series, Vol. 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes, 17
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Atri et al. (2020) Atri, P., Miller-Jones, J. C. A., Bahramian, A., et al. 2020, MNRAS, 493, L81, doi: 10.1093/mnrasl/slaa010
  * Bardeen (1970) Bardeen, J. M. 1970, Nature, 226, 64, doi: 10.1038/226064a0
  * Bardeen & Petterson (1975) Bardeen, J. M., & Petterson, J. A. 1975, ApJL, 195, L65, doi: 10.1086/181711
  * Bardeen et al. (1972) Bardeen, J. M., Press, W. H., & Teukolsky, S. A. 1972, ApJ, 178, 347, doi: 10.1086/151796
  * Barret et al. (2018) Barret, D., Lam Trong, T., den Herder, J.-W., & et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10699, Space Telescopes and Instrumentation 2018: Ultraviolet to Gamma Ray, ed. J.-W. A. den Herder, S. Nikzad, & K. Nakazawa, 106991G, doi: 10.1117/12.2312409
  * Belczynski et al. (2011) Belczynski, K., Bulik, T., & Bailyn, C. 2011, ApJ, 742, L2, doi: 10.1088/2041-8205/742/1/L2
  * Belczynski et al. (2008) Belczynski, K., Taam, R. E., Rantsiou, E., & van der Sluys, M. 2008, ApJ, 682, 474, doi: 10.1086/589609
  * Bhargava et al. (2021) Bhargava, Y., Belloni, T., Bhattacharya, D., Motta, S., & Ponti., G. 2021, MNRAS, 508, 3104, doi: 10.1093/mnras/stab2848
  * Blandford & Znajek (1977) Blandford, R. D., & Znajek, R. L. 1977, MNRAS, 179, 433, doi: 10.1093/mnras/179.3.433
  * Bolton (1972) Bolton, C. T. 1972, Nature Physical Science, 240, 124, doi: 10.1038/physci240124a0
  * Bonson & Gallo (2016) Bonson, K., & Gallo, L. C. 2016, MNRAS, 458, 1927, doi: 10.1093/mnras/stw466
  * Brenneman & Reynolds (2006) Brenneman, L. W., & Reynolds, C. S. 2006, ApJ, 652, 1028, doi: 10.1086/508146
  * Burrows et al. (2005) Burrows, D. N., Hill, J. E., Nousek, J. A., et al. 2005, Space Sci. Rev., 120, 165, doi: 10.1007/s11214-005-5097-2
  * Chen et al. (2016) Chen, Z., Gou, L., McClintock, J. E., et al. 2016, ApJ, 825, 45, doi: 10.3847/0004-637X/825/1/45
  * Corral-Santana et al. (2016) Corral-Santana, J. M., Casares, J., Muñoz-Darias, T., et al. 2016, A&A, 587, A61, doi: 10.1051/0004-6361/201527130
  * Dauser et al. (2014) Dauser, T., Garcia, J., Parker, M. L., Fabian, A. C., & Wilms, J. 2014, MNRAS, 444, L100, doi: 10.1093/mnrasl/slu125
  * Dauser et al. (2022) Dauser, T., García, J. A., Joyce, A., et al. 2022, MNRAS, 514, 3965, doi: 10.1093/mnras/stac1593
  * Dong et al. (2020) Dong, Y., García, J. A., Steiner, J. F., & Gou, L. 2020, MNRAS, 493, 4409, doi: 10.1093/mnras/staa606
  * Dong et al. (2022) Dong, Y., Liu, Z., Tuo, Y., et al. 2022, MNRAS, 514, 1422, doi: 10.1093/mnras/stac1466
  * Draghis et al. (2020a) Draghis, P., Miller, J. M., Balakrishnan, M., et al. 2020a, The Astronomer’s Telegram, 13665, 1
  * Draghis et al. (2021a) Draghis, P. A., Miller, J. M., Balakrishnan, M., et al. 2021a, The Astronomer’s Telegram, 14512, 1
  * Draghis et al. (2020b) Draghis, P. A., Miller, J. M., Cackett, E. M., et al. 2020b, ApJ, 900, 78, doi: 10.3847/1538-4357/aba2ec
  * Draghis et al. (2021b) Draghis, P. A., Miller, J. M., Zoghbi, A., et al. 2021b, ApJ, 920, 88, doi: 10.3847/1538-4357/ac1270
  * El-Batal et al. (2016) El-Batal, A. M., Miller, J. M., Reynolds, M. T., et al. 2016, ApJ, 826, L12, doi: 10.3847/2041-8205/826/1/L12
  * Event Horizon Telescope Collaboration et al. (2019) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2019, ApJ, 875, L1, doi: 10.3847/2041-8213/ab0ec7
  * Fabrika et al. (2015) Fabrika, S., Ueda, Y., Vinokurov, A., Sholukhova, O., & Shidatsu, M. 2015, Nature Physics, 11, 551, doi: 10.1038/nphys3348
  * Feng et al. (2021a) Feng, Y., Zhao, X., Li, Y., et al. 2021a, arXiv e-prints, arXiv:2112.02794.  https://arxiv.org/abs/2112.02794
  * Feng et al. (2021b) Feng, Y., Zhao, X., Gou, L., et al. 2021b, arXiv e-prints, arXiv:2109.07357.  https://arxiv.org/abs/2109.07357
  * Feng et al. (2021c) —. 2021c, arXiv e-prints, arXiv:2112.03479.  https://arxiv.org/abs/2112.03479
  * Fishbach & Kalogera (2022) Fishbach, M., & Kalogera, V. 2022, ApJ, 929, L26, doi: 10.3847/2041-8213/ac64a5
  * Fishbach et al. (2022) Fishbach, M., Kimball, C., & Kalogera, V. 2022, ApJ, 935, L26, doi: 10.3847/2041-8213/ac86c4
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 1, 24, doi: 10.21105/joss.00024
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Fragos & McClintock (2015) Fragos, T., & McClintock, J. E. 2015, ApJ, 800, 17, doi: 10.1088/0004-637X/800/1/17
  * Gallegos-Garcia et al. (2022) Gallegos-Garcia, M., Fishbach, M., Kalogera, V., Berry, C. P. L., & Doctor, Z. 2022, arXiv e-prints, arXiv:2207.14290.  https://arxiv.org/abs/2207.14290
  * García et al. (2014) García, J., Dauser, T., Lohfink, A., et al. 2014, ApJ, 782, 76, doi: 10.1088/0004-637X/782/2/76
  * García et al. (2015) García, J. A., Steiner, J. F., McClintock, J. E., et al. 2015, ApJ, 813, 84, doi: 10.1088/0004-637X/813/2/84
  * García et al. (2018) García, J. A., Steiner, J. F., Grinberg, V., et al. 2018, ApJ, 864, 25, doi: 10.3847/1538-4357/aad231
  * Gelman et al. (2004) Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B. 2004, Bayesian data analysis, 2nd edn., Texts in Statistical Science Series (Chapman & Hall/CRC, Boca Raton, FL), xxvi+668
  * Gendreau et al. (2012) Gendreau, K. C., Arzoumanian, Z., & Okajima, T. 2012, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 8443, Space Telescopes and Instrumentation 2012: Ultraviolet to Gamma Ray, ed. T. Takahashi, S. S. Murray, & J.-W. A. den Herder, 844313, doi: 10.1117/12.926396
  * George & Fabian (1991) George, I. M., & Fabian, A. C. 1991, MNRAS, 249, 352, doi: 10.1093/mnras/249.2.352
  * Gou et al. (2010) Gou, L., McClintock, J. E., Steiner, J. F., et al. 2010, ApJ, 718, L122, doi: 10.1088/2041-8205/718/2/L122
  * Gou et al. (2009) Gou, L., McClintock, J. E., Liu, J., et al. 2009, ApJ, 701, 1076, doi: 10.1088/0004-637X/701/2/1076
  * Gou et al. (2014) Gou, L., McClintock, J. E., Remillard, R. A., et al. 2014, ApJ, 790, 29, doi: 10.1088/0004-637X/790/1/29
  * Grefenstette et al. (2021) Grefenstette, B. W., Ludlam, R. M., Thompson, E. T., et al. 2021, ApJ, 909, 30, doi: 10.3847/1538-4357/abe045
  * Guan et al. (2021) Guan, J., Tao, L., Qu, J. L., et al. 2021, MNRAS, 504, 2168, doi: 10.1093/mnras/stab945
  * Hao & Zhang (2009) Hao, J. F., & Zhang, S. N. 2009, ApJ, 702, 1648, doi: 10.1088/0004-637X/702/2/1648
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Harrison et al. (2013) Harrison, F. A., Craig, W. W., Christensen, F. E., & et al. 2013, ApJ, 770, 103, doi: 10.1088/0004-637X/770/2/103
  * Hiemstra et al. (2011) Hiemstra, B., Méndez, M., Done, C., et al. 2011, MNRAS, 411, 137, doi: 10.1111/j.1365-2966.2010.17661.x
  * Hirai & Mandel (2021) Hirai, R., & Mandel, I. 2021, PASA, 38, e056, doi: 10.1017/pasa.2021.53
  * Homan et al. (2020) Homan, J., Bright, J., Motta, S. E., et al. 2020, ApJ, 891, L29, doi: 10.3847/2041-8213/ab7932
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Iben et al. (1995) Iben, Icko, J., Tutukov, A. V., & Yungelson, L. R. 1995, ApJS, 100, 217, doi: 10.1086/192217
  * Ishida et al. (2004) Ishida, M., Morio, K., & Ueda, Y. 2004, ApJ, 601, 1088, doi: 10.1086/380780
  * Jana et al. (2022) Jana, A., Chang, H.-K., Chatterjee, A., Naik, S., & Safi-Harb, S. 2022, arXiv e-prints, arXiv:2208.01399.  https://arxiv.org/abs/2208.01399
  * Jana et al. (2021a) Jana, A., Naik, S., Chatterjee, D., & Jaisawal, G. K. 2021a, MNRAS, 507, 4779, doi: 10.1093/mnras/stab2448
  * Jana et al. (2021b) Jana, A., Jaisawal, G. K., Naik, S., et al. 2021b, Research in Astronomy and Astrophysics, 21, 125, doi: 10.1088/1674-4527/21/5/125
  * Jia et al. (2022) Jia, N., Zhao, X., Gou, L., et al. 2022, MNRAS, 511, 3125, doi: 10.1093/mnras/stac121
  * Kaastra & Bleeker (2016) Kaastra, J. S., & Bleeker, J. A. M. 2016, A&A, 587, A151, doi: 10.1051/0004-6361/201527395
  * Kammoun et al. (2018) Kammoun, E. S., Nardini, E., & Risaliti, G. 2018, A&A, 614, A44, doi: 10.1051/0004-6361/201732377
  * King et al. (2014) King, A. L., Walton, D. J., Miller, J. M., et al. 2014, ApJ, 784, L2, doi: 10.1088/2041-8205/784/1/L2
  * Kolehmainen & Done (2010) Kolehmainen, M., & Done, C. 2010, MNRAS, 406, 2206, doi: 10.1111/j.1365-2966.2010.16835.x
  * Koljonen & Tomsick (2020) Koljonen, K. I. I., & Tomsick, J. A. 2020, A&A, 639, A13, doi: 10.1051/0004-6361/202037882
  * Krawczynski et al. (2022) Krawczynski, H., Muleri, F., Dovčiak, M., et al. 2022, arXiv e-prints, arXiv:2206.09972.  https://arxiv.org/abs/2206.09972
  * Liska et al. (2021) Liska, M., Hesp, C., Tchekhovskoy, A., et al. 2021, MNRAS, 507, 983, doi: 10.1093/mnras/staa099
  * Liu et al. (2008) Liu, J., McClintock, J. E., Narayan, R., Davis, S. W., & Orosz, J. A. 2008, ApJ, 679, L37, doi: 10.1086/588840
  * Lo & Lin (2011) Lo, K.-W., & Lin, L.-M. 2011, ApJ, 728, 12, doi: 10.1088/0004-637X/728/1/12
  * MacDonald et al. (2014) MacDonald, R. K. D., Bailyn, C. D., Buxton, M., et al. 2014, ApJ, 784, 2, doi: 10.1088/0004-637X/784/1/2
  * Madsen et al. (2017) Madsen, K. K., Christensen, F. E., Craig, W. W., et al. 2017, arXiv e-prints, arXiv:1711.02719.  https://arxiv.org/abs/1711.02719
  * Madsen et al. (2018) Madsen, K. K., Harrison, F., Broadway, D., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10699, Space Telescopes and Instrumentation 2018: Ultraviolet to Gamma Ray, ed. J.-W. A. den Herder, S. Nikzad, & K. Nakazawa, 106996M, doi: 10.1117/12.2314117
  * Maitra et al. (2014) Maitra, D., Miller, J. M., Reynolds, M. T., Reis, R., & Nowak, M. 2014, ApJ, 794, 85, doi: 10.1088/0004-637X/794/1/85
  * Mandel & Fragos (2020) Mandel, I., & Fragos, T. 2020, ApJ, 895, L28, doi: 10.3847/2041-8213/ab8e41
  * Margon & Anderson (1989) Margon, B., & Anderson, S. F. 1989, ApJ, 347, 448, doi: 10.1086/168132
  * Marshall et al. (2013) Marshall, H. L., Canizares, C. R., Hillwig, T., et al. 2013, ApJ, 775, 75, doi: 10.1088/0004-637X/775/1/75
  * Martí et al. (2002) Martí, J., Mirabel, I. F., Rodríguez, L. F., & Smith, I. A. 2002, A&A, 386, 571, doi: 10.1051/0004-6361:20020273
  * McCollough et al. (2016) McCollough, M. L., Corrales, L., & Dunham, M. M. 2016, ApJ, 830, L36, doi: 10.3847/2041-8205/830/2/L36
  * Michell (1784) Michell, J. 1784, Philosophical Transactions of the Royal Society of London, 74, 35.  http://www.jstor.org/stable/106576
  * Middleton et al. (2014) Middleton, M. J., Miller-Jones, J. C. A., & Fender, R. P. 2014, MNRAS, 439, 1740, doi: 10.1093/mnras/stu056
  * Middleton et al. (2012) Middleton, M. J., Sutton, A. D., Roberts, T. P., Jackson, F. E., & Done, C. 2012, MNRAS, 420, 2969, doi: 10.1111/j.1365-2966.2011.20145.x
  * Miller et al. (2009) Miller, J. M., Reynolds, C. S., Fabian, A. C., Miniutti, G., & Gallo, L. C. 2009, ApJ, 697, 900, doi: 10.1088/0004-637X/697/1/900
  * Miller et al. (2013) Miller, J. M., Parker, M. L., Fuerst, F., et al. 2013, ApJ, 775, L45, doi: 10.1088/2041-8205/775/2/L45
  * Miller et al. (2015) Miller, J. M., Tomsick, J. A., Bachetti, M., et al. 2015, ApJ, 799, L6, doi: 10.1088/2041-8205/799/1/L6
  * Miller et al. (2018) Miller, J. M., Gendreau, K., Ludlam, R. M., et al. 2018, ApJ, 860, L28, doi: 10.3847/2041-8213/aacc61
  * Molla et al. (2017) Molla, A. A., Chakrabarti, S. K., Debnath, D., & Mondal, S. 2017, ApJ, 834, 88, doi: 10.3847/1538-4357/834/1/88
  * Mondal (2010) Mondal, S. 2010, ApJ, 708, 1442, doi: 10.1088/0004-637X/708/2/1442
  * Mori et al. (2019) Mori, K., Hailey, C. J., Mandel, S., et al. 2019, ApJ, 885, 142, doi: 10.3847/1538-4357/ab4b47
  * Morningstar et al. (2014) Morningstar, W. R., Miller, J. M., Reis, R. C., & Ebisawa, K. 2014, ApJ, 784, L18, doi: 10.1088/2041-8205/784/2/L18
  * Mushotzky et al. (2019) Mushotzky, R., Aird, J., Barger, A. J., et al. 2019, in Bulletin of the American Astronomical Society, Vol. 51, 107.  https://arxiv.org/abs/1903.04083
  * Nealon et al. (2015) Nealon, R., Price, D. J., & Nixon, C. J. 2015, MNRAS, 448, 1526, doi: 10.1093/mnras/stv014
  * Novikov & Thorne (1973) Novikov, I. D., & Thorne, K. S. 1973, in Black Holes (Les Astres Occlus), ed. C. DeWitt & B. DeWitt, 343–450
  * Nowak et al. (2012) Nowak, M. A., Wilms, J., Pottschmidt, K., et al. 2012, ApJ, 744, 107, doi: 10.1088/0004-637X/744/2/107
  * Olejak & Belczynski (2021) Olejak, A., & Belczynski, K. 2021, ApJ, 921, L2, doi: 10.3847/2041-8213/ac2f48
  * Orosz et al. (2001) Orosz, J. A., Kuulkers, E., van der Klis, M., et al. 2001, ApJ, 555, 489, doi: 10.1086/321442
  * Pahari et al. (2018) Pahari, M., Bhattacharyya, S., Rao, A. R., et al. 2018, ApJ, 867, 86, doi: 10.3847/1538-4357/aae53b
  * Park et al. (2004) Park, S. Q., Miller, J. M., McClintock, J. E., et al. 2004, ApJ, 610, 378, doi: 10.1086/421511
  * Parker et al. (2015) Parker, M. L., Tomsick, J. A., Miller, J. M., et al. 2015, ApJ, 808, 9, doi: 10.1088/0004-637X/808/1/9
  * Parker et al. (2016) Parker, M. L., Tomsick, J. A., Kennea, J. A., et al. 2016, ApJ, 821, L6, doi: 10.3847/2041-8205/821/1/L6
  * Pérez & Granger (2007) Pérez, F., & Granger, B. E. 2007, Computing in Science and Engineering, 9, 21, doi: 10.1109/MCSE.2007.53
  * Picchi et al. (2020) Picchi, P., Shore, S. N., Harvey, E. J., & Berdyugin, A. 2020, A&A, 640, A96, doi: 10.1051/0004-6361/202037960
  * Pike et al. (2022) Pike, S. N., Negoro, H., Tomsick, J. A., et al. 2022, ApJ, 927, 190, doi: 10.3847/1538-4357/ac5258
  * Prabhakar et al. (2022) Prabhakar, G., Mandal, S., Athulya, M. P., & Nandi, A. 2022, MNRAS, 514, 6102, doi: 10.1093/mnras/stac1176
  * Pürrer et al. (2016) Pürrer, M., Hannam, M., & Ohme, F. 2016, Physical Review D, 93, 084042, doi: 10.1103/PhysRevD.93.084042
  * Ramachandran et al. (2022) Ramachandran, V., Oskinova, L. M., Hamann, W. R., et al. 2022, arXiv e-prints, arXiv:2208.07773.  https://arxiv.org/abs/2208.07773
  * Rao & Vadawale (2012) Rao, A., & Vadawale, S. V. 2012, ApJ, 757, L12, doi: 10.1088/2041-8205/757/1/L12
  * Reback et al. (2022) Reback, J., jbrockmendel, McKinney, W., et al. 2022, pandas-dev/pandas: Pandas 1.4.3, v1.4.3, Zenodo, doi: 10.5281/zenodo.6702671
  * Reeves et al. (2008) Reeves, J., Done, C., Pounds, K., et al. 2008, MNRAS, 385, L108, doi: 10.1111/j.1745-3933.2008.00443.x
  * Reis et al. (2012) Reis, R. C., Miller, J. M., Reynolds, M. T., Fabian, A. C., & Walton, D. J. 2012, ApJ, 751, 34, doi: 10.1088/0004-637X/751/1/34
  * Reis et al. (2013) Reis, R. C., Reynolds, M. T., Miller, J. M., et al. 2013, ApJ, 778, 155, doi: 10.1088/0004-637X/778/2/155
  * Reynolds (2021) Reynolds, C. S. 2021, ARA&A, 59, 117, doi: 10.1146/annurev-astro-112420-035022
  * Reynolds & Fabian (2008) Reynolds, C. S., & Fabian, A. C. 2008, ApJ, 675, 1048, doi: 10.1086/527344
  * Rout et al. (2020) Rout, S. K., Vadawale, S., & Méndez, M. 2020, ApJ, 888, L30, doi: 10.3847/2041-8213/ab629e
  * Salvesen et al. (2013) Salvesen, G., Miller, J. M., Reis, R. C., & Begelman, M. C. 2013, MNRAS, 431, 3510, doi: 10.1093/mnras/stt436
  * Schnittman et al. (2016) Schnittman, J. D., Krolik, J. H., & Noble, S. C. 2016, ApJ, 819, 48, doi: 10.3847/0004-637X/819/1/48
  * Schwarz (1978) Schwarz, G. 1978, Annals of Statistics, 6, 461
  * Schwarzschild (1999) Schwarzschild, K. 1999, arXiv e-prints, physics/9905030.  https://arxiv.org/abs/physics/9905030
  * Shafee et al. (2006) Shafee, R., McClintock, J. E., Narayan, R., et al. 2006, ApJ, 636, L113, doi: 10.1086/498938
  * Shafee et al. (2008) Shafee, R., Narayan, R., & McClintock, J. E. 2008, ApJ, 676, 549, doi: 10.1086/527346
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33
  * Sharma et al. (2021) Sharma, P., Sharma, R., Jain, C., Dewangan, G. C., & Dutta, A. 2021, Research in Astronomy and Astrophysics, 21, 214, doi: 10.1088/1674-4527/21/9/214
  * Shaw et al. (2022) Shaw, A. W., Miller, J. M., Grinberg, V., et al. 2022, MNRAS, doi: 10.1093/mnras/stac2213
  * Smirnov (1939) Smirnov, N. V. 1939, Bull. Math. Univ. Moscou, 2, 3
  * Smith et al. (2001) Smith, R. K., Brickhouse, N. S., Liedahl, D. A., & Raymond, J. C. 2001, ApJ, 556, L91, doi: 10.1086/322992
  * Sokal (1996) Sokal, A. 1996, in Monte Carlo Methods in Statistical Mechanics: Foundations and New Algorithms Note to the Reader
  * Spiegelhalter et al. (2002) Spiegelhalter, D. J., Best, N. G., Carlin, B. P., & Van Der Linde, A. 2002, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 64, 583, doi: https://doi.org/10.1111/1467-9868.00353
  * Sreehari et al. (2020) Sreehari, H., Nandi, A., Das, S., et al. 2020, MNRAS, 499, 5891, doi: 10.1093/mnras/staa3135
  * Stecchini et al. (2020) Stecchini, P. E., D’Amico, F., Jablonski, F., Castro, M., & Braga, J. 2020, MNRAS, 493, 2694, doi: 10.1093/mnras/staa417
  * Stefanov & Tasheva (2017) Stefanov, I. Z., & Tasheva, R. P. 2017, arXiv e-prints, arXiv:1711.11398.  https://arxiv.org/abs/1711.11398
  * Steiner et al. (2014) Steiner, J. F., McClintock, J. E., Orosz, J. A., et al. 2014, ApJ, 793, L29, doi: 10.1088/2041-8205/793/2/L29
  * Steiner et al. (2012) Steiner, J. F., McClintock, J. E., & Reid, M. J. 2012, ApJ, 745, L7, doi: 10.1088/2041-8205/745/1/L7
  * Steiner et al. (2010) Steiner, J. F., McClintock, J. E., Remillard, R. A., et al. 2010, ApJ, 718, L117, doi: 10.1088/2041-8205/718/2/L117
  * Steiner et al. (2016) Steiner, J. F., Walton, D. J., García, J. A., et al. 2016, ApJ, 817, 154, doi: 10.3847/0004-637X/817/2/154
  * Steiner et al. (2011) Steiner, J. F., Reis, R. C., McClintock, J. E., et al. 2011, MNRAS, 416, 941, doi: 10.1111/j.1365-2966.2011.19089.x
  * Tananbaum et al. (1972) Tananbaum, H., Gursky, H., Kellogg, E., Giacconi, R., & Jones, C. 1972, ApJ, 177, L5, doi: 10.1086/181042
  * Tao et al. (2019) Tao, L., Tomsick, J. A., Qu, J., et al. 2019, ApJ, 887, 184, doi: 10.3847/1538-4357/ab5282
  * Tashiro et al. (2018) Tashiro, M., Maejima, H., Toda, K., & et al. 2018, in Space Telescopes and Instrumentation 2018: Ultraviolet to Gamma Ray, ed. J.-W. A. den Herder, S. Nikzad, & K. Nakazawa, Vol. 10699, International Society for Optics and Photonics (SPIE), 520 – 531, doi: 10.1117/12.2309455
  * Tetarenko et al. (2016) Tetarenko, B. E., Sivakoff, G. R., Heinke, C. O., & Gladstone, J. C. 2016, ApJS, 222, 15, doi: 10.3847/0067-0049/222/2/15
  * The LIGO Scientific Collaboration et al. (2021a) The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration, et al. 2021a, arXiv e-prints, arXiv:2111.03606.  https://arxiv.org/abs/2111.03606
  * The LIGO Scientific Collaboration et al. (2021b) —. 2021b, arXiv e-prints, arXiv:2111.03634.  https://arxiv.org/abs/2111.03634
  * The LIGO Scientific Collaboration et al. (2021c) The LIGO Scientific Collaboration, the Virgo Collaboration, Abbott, R., et al. 2021c, arXiv e-prints, arXiv:2108.01045.  https://arxiv.org/abs/2108.01045
  * Thorne (1974) Thorne, K. S. 1974, ApJ, 191, 507, doi: 10.1086/152991
  * Tomsick et al. (2009) Tomsick, J. A., Yamaoka, K., Corbel, S., et al. 2009, ApJ, 707, L87, doi: 10.1088/0004-637X/707/1/L87
  * Tomsick et al. (2018) Tomsick, J. A., Parker, M. L., García, J. A., et al. 2018, ApJ, 855, 3, doi: 10.3847/1538-4357/aaaab1
  * Torres et al. (2020) Torres, M. A. P., Casares, J., Jiménez-Ibarra, F., et al. 2020, ApJ, 893, L37, doi: 10.3847/2041-8213/ab863a
  * Tursunov & Kološ (2018) Tursunov, A. A., & Kološ, M. 2018, Physics of Atomic Nuclei, 81, 279, doi: 10.1134/S1063778818020187
  * van den Heuvel (1976) van den Heuvel, E. P. J. 1976, in Structure and Evolution of Close Binary Systems, ed. P. Eggleton, S. Mitton, & J. Whelan, Vol. 73, 35
  * Vasudevan et al. (2016) Vasudevan, R. V., Fabian, A. C., Reynolds, C. S., et al. 2016, MNRAS, 458, 2012, doi: 10.1093/mnras/stw363
  * Verner et al. (1996) Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996, ApJ, 465, 487, doi: 10.1086/177435
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Vitale et al. (2017) Vitale, S., Gerosa, D., Haster, C.-J., Chatziioannou, K., & Zimmerman, A. 2017, Phys. Rev. Lett., 119, 251103, doi: 10.1103/PhysRevLett.119.251103
  * Walton et al. (2017) Walton, D. J., Mooley, K., King, A. L., et al. 2017, ApJ, 839, 110, doi: 10.3847/1538-4357/aa67e8
  * Wang et al. (2018) Wang, Y., Méndez, M., Altamirano, D., et al. 2018, MNRAS, 478, 4837, doi: 10.1093/mnras/sty1372
  * Webster & Murdin (1972) Webster, B. L., & Murdin, P. 1972, Nature, 235, 37, doi: 10.1038/235037a0
  * Wes McKinney (2010) Wes McKinney. 2010, in Proceedings of the 9th Python in Science Conference, ed. Stéfan van der Walt & Jarrod Millman, 56 – 61, doi: 10.25080/Majora-92bf1922-00a
  * Wilms et al. (2000) Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914, doi: 10.1086/317016
  * Xu et al. (2020a) Xu, Y., Harrison, F. A., Tomsick, J. A., et al. 2020a, ApJ, 893, 42, doi: 10.3847/1538-4357/ab7cdb
  * Xu et al. (2020b) —. 2020b, ApJ, 893, 30, doi: 10.3847/1538-4357/ab7dc0
  * Xu et al. (2018) Xu, Y., Harrison, F. A., Kennea, J. A., et al. 2018, ApJ, 865, 18, doi: 10.3847/1538-4357/aada03
  * Zdziarski et al. (2013) Zdziarski, A. A., Mikolajewska, J., & Belczynski, K. 2013, MNRAS, 429, L104, doi: 10.1093/mnrasl/sls035
  * Zevin et al. (2021) Zevin, M., Bavera, S. S., Berry, C. P. L., et al. 2021, ApJ, 910, 152, doi: 10.3847/1538-4357/abe40e
  * Zhao et al. (2021) Zhao, X., Gou, L., Dong, Y., et al. 2021, ApJ, 916, 108, doi: 10.3847/1538-4357/ac07a9

##  Appendix A Choice of source and background regions

The choice of background region adopted in this work is somewhat uncommon when
handling NuSTAR data. The common practice in the community is to adopt a
circular source region of a certain radius ($\sim 100^{\prime\prime}$) and a
circular background region of the same size, placed in an area of the detector
that is not contaminated by the source flux or any stray light that is being
detected (for a discussion on stray light in NuSTAR data, see Madsen et al.
2017; Grefenstette et al. 2021). This process requires the user to manually
identify the position of the source and background region in relation to the
image coordinates on the detector.

In order to simplify, automate, and increase the reproducibility of the
spectral extraction process, we adopted circular regions of radius equal to
120” for extraction of the source spectrum, centered at the position of the
source in the header of events file of the observation, for each of the
detectors. This choice for the position of the source region was made in order
to avoid the possibility that the coordinates available in literature for some
sources are uncertain and therefore bias the spectral extraction process. For
the background region, we used annular regions centered at the position of the
source, with an inner radius of 200” and an outer radius of 300”. Firstly,
this reduces the arbitrariness of the choice of position of the background
region. Secondly, this method samples a much larger detector area when
constructing the background spectrum. However, for bright sources, the 200”
inner radius of the annulus might be too small, leading to source photons
being included in the background. This is potentially important for sources
that have a low flux at high energies, as an increase in flux in the
background spectrum would lead to the source spectra becoming background
dominated at lower energies.

We tested this effect by extracting source and background spectra from the
observation 80502324004 of MAXI J0637-430, using five combinations of source
and background regions. The spectra obtained are shown in Figure 27, together
with the residuals of the fit to the best performing reflection model for this
observation (see Subsection 3.5). The error bars on the background points were
omitted for visual clarity. Additionally, we binned the spectra using three
different methods: the “optimal” binning scheme (Kaastra & Bleeker 2016), or
by requiring minimum a signal to noise ratio of 5 and of 10 in each spectral
bin. The five combinations of source and background regions are as follows:
Case 1 - 120” radius circular source region and annular background region with
200” inner radius and 300” outer radius, centered at the position of the
source; Case 2 - 120” radius circular source region and 120” radius circular
background region, placed as far as possible on the detector; Case 3 - 90”
radius circular source region and 90” radius circular background region,
placed as far as possible on the detector; Case 4 - 150” radius circular
source region and 150” radius circular background region, placed as far as
possible on the detector; Case 5 - 120” radius circular region and annular
background region with 300” inner radius and 600” outer radius, centered at
the position of the source.

The source spectrum is practically unchanged regardless of the size of the
source region, with no significant differences even at high energies. The
background spectra show changes between the five different cases, mainly at
low energies. The background spectra extracted using the region in Case 1
(used throughout the analysis in this paper) produces a background rate at 3
keV an order of magnitude larger than e.g., Case 2 and 3 (i.e. a 120” or 90”
circular background region placed far on the detector from the source region).
However, when comparing to Cases 4 and 5, the background rate is only $\sim 5$
times higher in Case 1. Given that the sources treated in this paper have high
fluxes at low energies, the difference in background spectra should not
influence the quality of the data and the obtained results. Given the goal of
automatizing of the spectral extraction process, Case 5 would produce the
ideal combination of source and background spectra. However, in Case 5, with
the significantly increased size of the annular extraction region, the risk of
the region overlapping with regions on the detector dominated by stray light
during observations is increased, likely influencing the background spectra
more than the difference between Case 1 and 5 would.

At high energies, the background spectra are consistent regardless of choice
of extraction region. However, different binning schemes can lead to the
appearance of the source spectra becoming background dominated at lower
energies, the clearest being Case 3 when requiring all bins to have a minimum
signal to noise ratio of 10. For the same case, different binning techniques
clearly indicate that the background spectra dominate over the source spectra
above $\sim 60$ keV. Interestingly, the binning techniques requiring a minimum
signal to noise ratio often have difficulties in treating the background
spectra, creating apparent streaks. However, the error bars on the fluxes in
the background bins are larger than the variation trend seen in Cases 1-4.
Nevertheless, it appears that optimal binning better handles background
spectra.

With the different cases of source and background region and different binning
schemes producing visually no difference at high energies and changes of
little importance at low energies, it is important to quantify the difference
produced when fitting the data with a model. We fit all spectra shown in
Figure 27 with the best performing relxill model for this particular
observation in our analysis. The best-fit values for the temperature of the
accretion disk, BH spin, inner accretion disk inclination, power-law index,
and disk ionization are shown in Figure 28. The different colors of the points
represent the five cases presented in Figure 27, and the different marker
represent the binning scheme used, with circles, squares, and diamonds
representing the optimal binning scheme, the signal to noise ratio of 5 and of
10, respectively. The filled markers indicate the results of fitting the
spectra in the 3-60 keV band, while the empty markers indicate the results of
fitting the spectra in the entire 3-79 keV NuSTAR band. The left panels show
only the value of the parameters for which the best-fit solution was achieved,
while the right panels also show the uncertainty of the parameters, computed
as the square root of the diagonal element in the covariance matrix of the
fit. This uncertainty is used to produce the proposal distribution for the
walkers in the MCMC analysis.

It is important to mention that these measurements are made only by fitting
one of the two NuSTAR spectra obtained during one observation (FPMA in this
case), and fitting both spectra jointly is likely to reduce the variation
between the parameters and to decrease the size of the uncertainties shown in
the right panels. Truncating the spectra at 60 keV as opposed to 79 keV has no
obvious impact on the best-fit parameters, with the only noticeable difference
being that fitting the spectra up to 79 keV produces slightly smaller
uncertainties on some parameters. The difference between the best-fit
parameters found when fitting spectra extracted using the five different cases
for the region size, shape, and location is comparable with the difference
produced when using different binning techniques. However, the variability in
the best-fit parameters is negligible when compared to the size of the errors.
As these fit parameters are only used for the proposal distribution in the
MCMC analysis, the results will not be significantly influenced by this
variability induced by the different choices of source and background region
and of binning scheme used.

Therefore, throughout the analysis, we continue using Case 1 for extraction of
source and background spectra, as it increases the ability to automate our
spectra extraction pipeline when compared to Cases 2-4, and it reduces the
risk of producing biased background spectra due to stray light when compared
to Case 5, at the cost of producing increased background spectra at low
energies, which appear to have no significant impact on the result. We
continue using the optimal binning scheme, as it leads to similar results with
the other binning techniques, while producing spectra with fewer spectral bins
compared to the other techniques, which helps reduce the time required to fit
the spectra. We mention that the optimal binning scheme still produces bins
with enough counts that using $\chi^{2}$ statistics is still appropriate.

##  Appendix B Algorithm to combine measurements

The script used to combine measurements takes the posterior distributions of
the individual parameters obtained through the MCMC runs and maps them to the
$[0,1]$ interval. Then, it then generates a beta distribution with the
probability density function described by:

$PDF(x)=\frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha,\beta)}$ (B1)

where
$B(\alpha,\beta)=\frac{\Gamma(\alpha)\Gamma(\beta)}{\Gamma(\alpha+\beta)}$ and
$\Gamma$ is the Gamma function, $\Gamma(a)=\int_{0}^{\infty}t^{a-1}e^{-a}dt$.
This distribution has two parameters, $\alpha$ and $\beta$, that we determine
through a Bayesian analysis in order to describe the combined measurement. We
assume uniform priors for $\alpha$ and $\beta$, but we require that
$\alpha>1$, $\beta>1$, and $(\alpha+\beta)<171.5$. The first two requirements
come from the definition of the beta distribution, and the third is dictated
by numerical limitations of the used programming language: $\Gamma(171.5)\sim
9.5\times 10^{307}$. This number is so large that when computing the value of
$B(\alpha,\beta)$, any combination of $\alpha,\beta$ for which
$\alpha+\beta\geq 171.5$ returns $B(\alpha,\beta)=0$ due to division by a very
large number. This then leads to division by zero when computing the value of
the PDF at a given location.

The total likelihood of a parameter combination is evaluated by summing the
(logarithm of the) likelihood of the parameter combination describing each
individual measurement. As the posterior distributions for the parameter are
normalized in the same way as the Beta distribution, we quantify how well the
Beta distribution together with a parameter combination describes the spin
measurement by taking the ratio of the are under the minimum of the two curves
to the area under the maximum of the two curves, evaluated over the $[0,1]$
interval. When computing the total likelihood by summing the (logarithm of
the) individual likelihoods for each measurement, we weight each measurement
in the sum by the ratio of the reflected flux to the total flux in the 3-79
keV band during the observation. This way, we weight observations with
stronger reflection more than those with weaker reflection, as they are more
likely to allow us to separate the reflection component from the underlying
continuum and better constrain the parameters.

As an example of how our likelihood function behaves, Figure 29 illustrates
our attempt to fit a Gaussian function defined on the $[0,1]$ interval, with a
mean of $\mu=0.8$ and a standard deviation $\sigma=0.1$. The left panels show
the Gaussian function (in black) and three attempted models shown in green: a
beta distribution (top panel), Delta distribution (middle panel), and a
uniform distribution (bottom panel). All curves are normalized in the same
way. On the right, we show the minimum of the two curves in blue and the
maximum of the two in red. In order to maximize likelihood, our function
attempts to maximize the minimum of the two curves and to minimize the maximum
of the two. The method becomes more intuitive upon inspection of the top right
panel in Figure 29: our code aims to increase the blue curve (minimum of the
two) and to decrease the red one (maximum of the two). This way, the minimum
and maximum are most similar, translating to a better fit of our beta
distribution to the initial assumed distribution (in this example) and to the
data (in real usage).

To test the ability of our script to combine measurements and recover the
underlying distribution, we adopted a Gaussian distribution of some given mean
and standard deviation and generated a number of Gaussian measurements with
position probabilities given by the initial assumed probability density
function. We then used our script to combine all “measurements” in order to
infer the underlying distribution, with the main goal of recovering its mode
and credible interval.

Four examples are shown in the four panels of Figure 30. The black curves
represent the assumed underlying Gaussian distributions, and the dashed green
curves show the Beta distributions representing the best fit of our algorithm
to the distribution underlying all “measurements”. The top left panel
simulates an attempt to combine ten measurements around a mean of $\mu=0.7$,
with a standard deviation $\sigma=0.1$. Both the mode and the credible
interval of our inferred distribution matches the true values very closely. We
note that due to numerical artifacts caused by the limited number of
computation iterations, the truth mode and credible intervals of the assumed
distribution (black curve) differ very slightly from the values assumed when
generating the distributions. This example is well applicable to combining
multiple inclination measurements. While the example shows distributions
defined on the $[0,1]$ interval, one could combine inclination measurements
through a change of variable.

The top right panel shows an attempt to combine five measurements generated
around a mean of $\mu=0.97$ and standard deviation $\sigma=0.05$. This example
illustrates the usability of combining spin measurements clustered around the
maximal value. In this case, the inferred distribution (green dashed line)
again closely recovers the assumed underlying distribution, with the mode and
credible intervals of the assumed and derived distributions being in good
agreement. In this case, the derived distribution appears to be slightly
shifted toward lower spin values, likely due to the requirement of the Beta
distribution to return to a value of zero at the edge of the interval.

The bottom panels of Figure 30 show an attempt to combine three measurements
drew from a Gaussian of mean $\mu=0.5$ and standard deviation $\sigma=0.05$.
The left panel shows the combined distribution when all three measurements are
weighted equally, well recovering the assumed distribution. The right panel
weights the second measurement (shown through the solid green line) eight
times more than the other two measurements. The inferred distribution much
better recovers the second “measurement”, which was more heavily weighted,
while also slightly accounting for the other two measurements.

Often, when multiple observations of the same source are available, an option
is to fit the spectra jointly, linking parameters that are expected to stay
constant between observations. In this case, the BH spin and the inclination
of the inner accretion disk could be linked between observations. The
constraints on the parameters would simultaneously account for all the data
that is being fit, providing more credible measurements. However, fitting
multiple NuSTAR spectra with complicated models quickly becomes a task too
computationally intensive, which requires unfeasible times to run. In order to
test the robustness of our result combining method against the possibility of
fitting the spectra jointly with parameters linked, we fit two of the three
observations of MAXI J0637-430 used in the analysis presented in Subsection
3.5 jointly. We fit the spectra of the two observations in the 3-50 keV band,
linking the BH spin, the inclination of the inner disk, and the column density
along the line of sight, and allowed all other parameters to vary freely, as
described in Subsection 2.1.

Fitting the spectra jointly using TBabs*(diskbb+powerlaw) returns
$\chi^{2}/\nu=1651.15/634$, with clear signs of reflection. Similarly to the
individual fits, the best performing models in terms of $\chi^{2}$ that
accounts for refection are TBabs*(diskbb+relxill), returning
$\chi^{2}/\nu=707.67/619$ and TBabs*(diskbb+relxillD) with Log(N)=19,
returning $\chi^{2}/\nu=710.5/622$. The latter performs significantly better
in terms of DIC, with $\Delta$DIC=200 better than the former. Therefore, we
report the measurements based on the TBabs*(diskbb+relxillD) with Log(N)=19
MCMC run. We find $a=0.980^{+0.005}_{-0.009}$ and $\theta=59^{+3}_{-2}$
degrees. In comparison, when fitting the spectra individually and combining
the results using our algorithm, we find $a=0.97\pm 0.02$ and
$\theta=62^{+3}_{-4}$ degrees. The results formally agree well within their
uncertainties, suggesting that our method produces results that are robust and
precise. However, it is worth noting that the combined result also takes
advantage of the third available observation. The quoted statistical
uncertainties on both the spin and inclination measurements are larger when
combining the measurements when compared to the measurements based on a joint
fit. This is linked to the fact that the errors quoted are simply statistical,
and the combined measurement begins to include systematical errors, which are
likely to contribute significantly.

##  Appendix C Sources wherein no measurement was possible

###  C.1 XTE J1859+226

NuSTAR observed XTE J1859+226 once (ObsID 90701305002) for 41 ks. This
observation was described in Draghis et al. (2021a). The spectra are
background dominated above $\sim 20\;\rm keV$. Fitting with the power law
model returns $\chi^{2}/\nu=187/208$ when the data is binned to have a minimum
signal/noise ratio SNR=5, and $\chi^{2}/\nu=304/340$ when using optimal
binning. No reflection features are apparent and a spin measurement is not
possible.

###  C.2 4U 1755-338

4U 1755-338 was observed once by NuSTAR (ObsID 90601313002), and the
observation was first described by Draghis et al. (2020a). During the
observation, the source was in a “high/soft” state, with the emission being
dominated by the accretion disk. The source is detected up to an energy $\sim
25\;\rm keV$. Fitting the spectra using diskbb+powerlaw produces a good fit,
returning $\chi^{2}/\nu=416/374$. No reflection features are apparent.

###  C.3 M31 ULX-1 (CXOM31 J004253.1+411422)

The position of M31 ULX-1 was covered by 11 NuSTAR observations: ObsID
50026001002, 50026001004, 50111001002, 50110001002, 50110001004, 50101001002,
50110001006, 50201001002, 50302001002, 50302001004, 50302001006. The source is
not detected in any of the 11 observations.

###  C.4 M31 ULX-2 (XMMU J004243.6+412519)

There are 4 NuSTAR observations of the region containing M31 ULX-2: ObsIDs
50302001002, 50302001004, 50302001006, and 30365002002. The source is not
detected in any of the four observations.

###  C.5 XTE J1748-288

NuSTAR observation 40012019001 was pointed at $\sim 8^{\prime}$ from XTE
J1748-288, placing the known location of the source close to the edge of the
detectors. However, the source was quiescent during the NuSTAR observation and
not detected.

###  C.6 SS 433

SS 433 is an X-ray binary system that is best known for launching extended
relativistic jets (v = 0.26c, Margon & Anderson 1989) that are loaded with
baryons (see, e.g., Marshall et al. 2013). The X-ray spectrum of the source is
dominated by emission from the variable and precessing jets, complicating
efforts to detect emission from the inner disk. Moreover, the jets may be
launched because SS 433 harbors a compact object that is accreting far above
its Eddington limit (e.g. Fabrika et al. 2015). Especially in this case, the
central engine may be enshrouded within an optically-thick spherical flow, or
at least partially obscured by a vertically extended disk that further serves
to obscure signatures from the innermost accretion flow. Therefore, despite
emerging evidence that SS 433 may harbor a black hole rather than a neutron
star (e.g., Picchi et al. 2020), SS 433 fails several of our criteria for
examination, and it is excluded from our analysis.

###  C.7 XTE J0421+560

XTE J0421+560 is the X-ray counterpart of CI Camelopardalis, which is believed
to be a B[e] star. The nature of the compact object in the system is not
constrained, and Ishida et al. (2004) argued that the emission from the system
during quiescence does not allow placing constraints on the nature of the
compact object, but that the ASCA observations taken during the 1998 outburst
of the source broadly point toward the compact object being a white dwarf.
Constraints on the nature of a compact object can be placed through spin
measurements: a spin higher than $a=0.7$ (Lo & Lin 2011) indicates that it
must be a black hole (see e.g., Draghis et al. 2022 in prep for an example
analysis). Therefore, we chose to analyze this source despite a lack of clear
evidence of it harboring a black hole.

The spectra from ObsID 90201040002 of XTE J0421+560 do not require a diskbb
component. We fit the spectra in the 3-20 keV band. Fitting with powerlaw
gives $\chi^{2}/\nu=103.2/95$. The two best performing models are
TBabs*relxillD with Log(N)=19 with $\chi^{2}/\nu=92.37/87$ and TBabs*relxill
with $\chi^{2}/\nu=93.41/87$.

Running MCMC on both models, the chains produce similar DIC, 109.07 and 109.42
respectively. In both cases, the spin is low, $a<0.26$. However, running
steppar between -0.998 and 0.998 in spin shows differences of
$\Delta\chi^{2}=0.1$ across the entire parameter space. Even though reflection
is evident in this observation, the spin cannot be constrained given the
quality of the data.

###  C.8 IC 10 X-1

There is one NuSTAR observation of IC 10 X-1 (ObsID 30001014002). Using the
continuum fitting method, Steiner et al. (2016) measured the spin of the BH in
IC 10 X-1 to be $a=0.85^{+0.04}_{-0.07}$ (90% confidence).

Fitting the 3-20 keV spectrum with TBabs*(diskbb+powerlaw) produces a good
fit, with $\chi^{2}/\nu=111.91/96$. Replacing the power-law component with
models from the relxill family significantly improves the fit, with the best
performing model TBabs*(diskbb+relxilCp) returning $\chi^{2}/\nu=81.10/86$. We
find that the fits are unable to constrain the BH spin. 68% of the posterior
samples in the MCMC run take spin values $a\leq 0.17$. This measurement could
be biased towards low values by constraining the accretion disk to extend all
the way to the ISCO. Still, assuming the mass and distance estimates from
Steiner et al. (2016) ($M\sim 15\;\rm M_{\odot}$ and $d\sim 800~{}\rm kpc$ and
given the 3-20 keV flux in the NuSTAR observation, this would result in an
Eddington fraction of $\sim 3\%$, well within the expected theoretical range
for which the inner disk is expected to extend to the ISCO. Additionally, the
environment around the black hole is likely to be significantly more
complicated than the simplistic assumption of this model. As the signal to
noise of the observation is low, we do not pursue more complicated models.

###  C.9 M33 X-7

There are 2 NuSTAR observations that cover the position of the source (ObsID
50310001002 and 50310001004). M33 X-7 is not detected during the two
observations, taken in March and July 2017. No spin measurement was possible.

###  C.10 M33 X-6

In the same two NuSTAR observations as for M33 X-7 (ObsID 50310001002 and
50310001004), the source M33 X-6 is clearly detected. For ObsID 50310001002 we
fit the spectra in the 3-15 keV band, while for ObsID 50310001004 we fit in
the 3-10 keV band. For the former, fitting the 3-15 keV spectra with
TBabs*powerlaw returns $\chi^{2}/\nu=95.91/80$, while the best performing
reflection model TBabs*relxilllp returns $\chi^{2}/\nu=90.75/73$. Similarly,
for ObsID 50310001004, fitting with TBabs*powerlaw produces a good fit,
returning $\chi^{2}/\nu=54.72/50$, while the best reflection model
TBabs*relxill returns $\chi^{2}/\nu=50.48/41$, with the improvement in
statistic not being justified by the increase in the number of free
parameters. In both cases, the addition of a reflection component is not
justified statistically.

###  C.11 M33 X-8

There are 2 NuSTAR observations of M33 X-8 (ObsID 50310002001 and 50310002003)
in which the source is clearly detected. Fitting the 50310002001 spectra in
the 3-15 keV band using TBabs*(diskbb+powerlaw) returns
$\chi^{2}/\nu=118.97/89$, while the best performing reflection model
TBabs*(diskbb+relxilCp) returns $\chi^{2}/\nu=116.20/80$. For the 50310002003
observation, also fitting in the 3-15 keV band, the TBabs*(diskbb+powerlaw)
model returns a good fit, with $\chi^{2}/\nu=60.72/86$, while the best
performing reflection model TBabs*(diskbb+relxilD) with Log(n)=19 fixed
returns $\chi^{2}/\nu=55.89/78$. The statistical improvement owed to the
addition of reflection is not justified by the increase in the number of free