Systematically Revisiting All NuSTAR Spins of Black Holes in X-Ray Binaries

We extend our recent work on black hole spin in X-ray binary systems to include an analysis of 189 archival NuSTAR observations from 24 sources. Using self-consistent data reduction pipelines, spectral models, and statistical techniques, we report an unprecedented and uniform sample of 36 stellar-mass black hole spin measurements based on relativistic reflection. This treatment suggests that prior reports of low spins in a small number of sources were generally erroneous: our comprehensive treatment finds that those sources tend to harbor black holes with high spin values. Overall, within 1σ uncertainty, ∼86% of the sample are consistent with a ≥ 0.95, ∼94% of the sample are consistent with a ≥ 0.9, and 100% are consistent with a ≥ 0.7 (the theoretical maximum for neutron stars; a = cJ/GM 2). We also find that the high-mass X-ray binaries (those with A-, B-, or O-type companions) are consistent with a ≥ 0.9 within the 1σ errors; this is in agreement with the low-mass X-ray binary population and may be especially important for comparisons to black holes discovered in gravitational wave events. In some cases, different spectra from the same source yield similar spin measurements but conflicting values for the inclination of the inner disk; we suggest that this is due to variable disk winds obscuring the blue wing of the relativistic Fe K emission line. We discuss the implications of our measurements, the unique view of systematic uncertainties enabled by our treatment, and future efforts to characterize black hole spins with new missions.


INTRODUCTION
The "relativistic reflection" spin measurement technique (Fabian et al. 2000) has been extensively used to measure black hole (BH) spin1 across the entire mass range, from stellar-mass BHs in X-ray binary (XB) systems (see, e.g.Brenneman & Reynolds 2006;Miller et al. 2009;Draghis et al. 2023c) to the BHs in active galactic nuclei (AGN; see, e.g., Gallo et al. 2015;Keck et al. 2015).For a review of BH spin measurement tech-niques, see Reynolds (2021).Furthermore, NuSTAR (Harrison et al. 2013) observations taken over the past decade have enhanced the findings of the method, by providing spectra ideally suited for such studies, which cover the main features of relativistic reflection (the Fe K complex and the Compton "hump") with high sensitivity.Relativistic reflection studies on NuSTAR data have led to more than 30 BH spin measurements in XB systems, with more sources being discovered, observed, and analyzed every year.However, as the studies occurred over a decade and were performed by independent authors, as understanding of the physical mechanisms present in the analyzed systems grew and our models matured, and as our ability to quantify instrumental effects evolved, a comparison between indepen-dent measurements or a uniform analysis of the entire BH distribution would be unproductive.
In contrast, gravitational wave (GW) observations have so far detected ∼ 90 binary black hole (BBH) merger candidates, translating to roughly twice as many BHs observed and analyzed in the third edition of the Gravitational Wave Transient Catalog (GWTC-3;Abbott et al. 2023a).The efforts of the GW community have produced a uniform description of the entire observed BH sample (Abbott et al. 2023b), leading to an inference on the population of merging BHs.This study suggests that the BHs observed to merge in BBH systems have preferentially low spins, in contrast with the BH spin distribution observed for the BHs in XBs, hinting at the fact that different populations of BHs are observed through GWs and electromagnetic radiation.In this regard, Gallegos-Garcia et al. ( 2022) concluded that at most 20% of BBHs evolve from high-mass X-ray binary (HMXB) systems formed through Case-A mass transfer and fewer than 11% of Case-A HMXB systems evolve into BBHs, further hinting at the different origins of the BH populations.In order to understand BH formation and evolution as a whole, any theory needs to be able to reproduce both the high spins observed in XBs and the low spins observed in BBHs.However, before taking such steps, a uniform analysis of the entire population of BH spins in XB systems is necessary, similarly to the case of BBHs observed through GWs.
A significant topic of discussion within the field of Xray measurements of BH spins is the magnitude of the systematic uncertainties of such measurements, with existing studies reporting only the statistical uncertainties of the measurements.A few of the most significant possible sources of systematic uncertainty come from (1) disagreement between existing families of models, or between different flavors from the same family of models; (2) the validity of the assumptions regarding the physical processes occurring in the observed systems, and the accuracy with which the models capture the extent of the source behavior; (3) small-scale source variability, which is unaccounted for by the models used to fit the data; and (4) an incomplete understanding of the models used, and of possible correlations between sets of parameters which produce broad regions of the parameter space that are not properly explored.Furthermore, when attempting to describe the entire population of BH spins in XBs, it is important to acknowledge that so far, only a small subset of such systems in our Galaxy have been observed, and to consider the existence of possible selection effects and observational biases.Works such as Draghis et al. (2020) and Draghis et al. (2021) attempt to address the first point highlighted above, by exploring a large array of models when fitting the data, composed of multiple combinations of different model components.The second point in the list will be addressed in time, as more sources will be observed with instruments such as NuSTAR or NICER (Gendreau et al. 2016), and in the more distant future with observations from proposed missions such as AXIS (Mushotzky 2018), HEX-P (Madsen et al. 2018), or STROBE-X (Ray et al. 2019).In the near future, pairing NuSTAR and NICER observations with high-resolution spectra obtained using the Resolve instrument on board XRISM (Tashiro et al. 2020) will further help in evolving our understanding of physical mechanisms that might be unaccounted for during a relativistic reflection study, and their importance on the precision of BH spin measurements.
In Draghis et al. (2023b), we began addressing the second possible source of systematic uncertainties listed above, and hinted at the possibility of taking steps in the direction to address the third point listed above.In this paper, we continue the effort of trying to incorporate the information available throughout all existing NuS-TAR observations in BH spin measurements, and lay the groundwork for a large-scale analysis of the entire parameter space across a large number of observations.In that regard, we analyze 189 archival NuSTAR spectra of 24 sources that have previous spin measurements using the relativistic reflection method using NuSTAR data.For consistency with previous work, we performed this analysis in a manner that follows the pipeline highlighted in Draghis et al. (2023b).The best-fit parameters from the array of spectral fits will be analyzed in a future paper to systematically investigate the behavior of the models used under a variety of observational cases.The work presented in this paper, together with previous results, creates a sample of 36 sources with spins measured using all the available data, by using fully consistent methods and making the same sets of assumptions.Through this large data-set of BH spins measured using a uniform treatment, we set the groundwork for properly deriving the BH spin distribution in XB systems.
This paper is structured as follows.In Section 2, we briefly review the methods used to extract spectra from the NuSTAR observations, we explain the criteria used to select the observations analyzed, we discuss the assumptions made by our models, and briefly describe the analysis pipeline.In Section 3, we present an analysis and results of the BH spin measurements for each of the 24 sources treated in this paper.Lastly, in Section 4 we discuss the findings of this large-scale project, their implications, and future directions for improving our understanding of BH spin and of BH formation as a whole.

OBSERVATIONS, MODELS, AND ANALYSIS
For most observations, we extracted spectra using the standard nustardas v2.1.1 routines in Heasoft v6.29c, using the NuSTAR CALDB version 20211103.Unless otherwise stated, we extracted the source spectra from circular regions of 120" radius, centered at the position of the source, and background spectra from an annulus centered on the source, with inner radius of 200" and outer radius of 300".We rebinned the spectra using the "optimal" binning scheme (Kaastra & Bleeker 2016) through the ftgrouppha ftool, and only fit the spectra in the energy band where the source spectra are significantly above the background spectra, also ensuring that usage of the χ 2 statistic is appropriate.Background spectra are subtracted from the source spectra during spectral fitting.We performed the spectral analysis using the Xspec v12.12.0g (Arnaud 1996) software, ensuring wilm abundances (Wilms et al. 2000) and vern photoelectric absorption cross sections (Verner et al. 1996).Instead of taking the customary approach of allowing a constant offset between the spectra from the two NuS-TAR focal plane module (FPM) detectors, we allowed the normalizations of the components to vary freely.
We only fit observations for which the Eddington fraction of the source during the observation was between 10 −3 ≤ L/L Edd ≤ 0.3, under the assumption that for sources within this Eddington fraction range, the accretion disk is not truncated and extends to the innermost stable circular orbit (ISCO).If measurements of BH mass and distance to the source were available in the literature, we used them to estimate the Eddington fraction of the source during the observation.For the remaining nine sources where estimates of BH mass and distance were unavailable, we assumed generic values: M BH ∼ 8 ± 4M ⊙ and d ∼ 10 ± 5kpc.We note that while the assumption of an accretion disk extending to the ISCO is backed by theoretical, numerical, and observational results, the lower limit of the Eddington fraction for which this behavior is expected is a topic of debate, with many works suggesting significant disk truncation, especially in hard spectral states (see reviews such as Done et al. 2007;Bambi et al. 2021).If disk truncation is present in the faint hard state, the choice of Eddington fraction for which this assumption holds is unlikely to significantly influence the final reported results for the measured spin, as low flux observations generally produce poor spin constraints, which contribute less to the overall spin probability distribution obtained by combining information from all available observations.Further-more, it is likely that the lower limit of the Eddington fraction for which the disk extends to the ISCO is not a universal number, applicable for all sources.Therefore, for consistency with numerous amounts of existing work, we continue to employ this assumption, despite acknowledging its limitations and potential effects.
We first fit all observations with a "baseline" model, which includes Galactic absorption along the line of sight, simplistic emission from a thermal accretion disk, and power-law emission from a compact corona: TBabs*(diskbb+powerlaw).Then, we fit all spectra with an array of six models that replace the powerlaw component in our baseline model with different flavors of the relxill v.1.4.3 family of models (Dauser et al. 2014;García et al. 2014).We test the quality of the fits using the relxill, relxillCp, relxilllp, and the relxillD flavor, with the accretion disk density fixed to log(n) = 10 15 , 10 17 , and 10 19 cm −3 .It is important to note that the relxill family of models only considers scattering by thermal electrons, and current-generation data are unable to constrain the need for scattering by a population of nonthermal electrons.Furthermore, fits to our current-generation data appear to be unable to constrain the need for a component describing the Comptonization of disk photons, which could impact the ability to disentangle the reflected radiation from the underlying continuum, potentially biasing the spin measurements.In the future, pairing NuSTAR data with high-resolution XRISM spectra will enable probing the effects that these effects have on our ability to constrain BH spin.
For spectra that show the presence of additional narrow absorption and emission features, we introduce the multiplicative component zxipcf to the models, which describes the shape of the spectrum after being partially covered by ionized material.Furthermore, when informed by the data, we introduce the gaussian model component with negative (or positive, when describing narrow emission features) normalization to describe narrow absorption features around 7 keV, likely related to absorption by an ionized outflow (or possible residual distant refection when positive -for Cygnus X-1 and MAXI J1848-015).In a few instances, our models include both the zxipcf and gaussian components.When the diskbb component is not required by the data, the normalization of the component and/or the temperature of the disk take low values, which ensure that no significant contribution from the component is present in the total model.However, in a few situations, in order to reduce the complexity of the model, we remove the diskbb component from the models when it is not required by the data.This happens for a few of the spectra of sources in a hard state, especially when the models also require the zxipcf and gaussian components.
In this analysis, we fit each observation individually.Upon obtaining the best-fit combination of parameters for the models for each spectrum, we run a Markov Chain Monte Carlo (MCMC) analysis as described in Section 2.2 in Draghis et al. (2023b), and combine the individual 1D posterior probability distributions using a method similar to the one described in Draghis et al. (2023b), but as expanded upon in Draghis et al. (2023a).When running the combining algorithm, we weigh the independent measurements by the strength of reflection during the observation.We report the outputs of this Bayesian combining algorithm as our results, for the BH spin and the viewing inclination angle.As all the sources treated in this paper have been analyzed in previous works and as the goal of this paper is to quantify the information regarding the BH spin from all the available observations together, we do not report the full list of measured parameters for each observation.However, this entire data-set will be analyzed in a future paper and will be published at that time, and we note that it can be made available by request to the authors.

AT 2019wey
At the time of the analysis, three NuSTAR observations of AT 2019wey were public, all taken in 2022 (ObsIDs 90601315002, 90601315004, and 90601315006).Feng et al. (2022b) analyzed the third observation using the relxill family of models to find a BH spin of a ∼ 0.96 − 0.97 and an inclination θ < 30 • .Given that the distance estimate to the source is uncertain, 1 ≲ d ≲ 10 kpc (Yao et al. 2021), and that the BH mass is unknown (we assumed M = 8 ± 4 M ⊙ ), all three observations fall within the range of Eddington fractions for which we expect the accretion disk to extend to the ISCO.Therefore, we analyzed all three existing observations of the source.The top panels in Figure 7 show the unfolded FPMA and FPMB NuSTAR spectra obtained during the three observations, while the middle and bottom panels show the residuals produced when fitting the spectra with a power-law model and with the best-fitting reflection models, respectively, together with the statistic produced.
Upon running the MCMC analysis on the bestperforming reflection models, we obtained the posterior probability distributions shown by the colored lines in the left (for spin) and right (for inclination) panels of Figure 8.By combining the individual measurements using our prescription, which accounts for system-atic variations between observations, we obtain a spin measurement of a = 0.906 +0.084 −0.202 and an inclination of θ = 14 +12 −10 [ • ].Our measurements formally agree with those of Feng et al. (2022b), indicating that our spin constraint is accurate, but less precise.The increase in the uncertainty of the measurement comes from incorporating information from all existing observations in a way that better accounts for systematic variations.We note that if only the posterior distribution for the third observation was considered (red curve in the left panel of Figure 8), the spin measurement would be much more narrowly constrained and at nearly maximal values, similar to the result of Feng et al. (2022b).

LMC X-3
At the time of the analysis, there were eight NuSTAR public observations of LMC X-3.As the source is relatively faint and persistently in a disk-dominated state, we only report the analysis of spectra which have shown a statistical improvement in fit statistic when including reflection features in the models: ObsIDs 10601308002, 10601308004, 30402035002, 30402035004, 30402035008, and 30902041002.We note that using a mass of M = 6.98 ± 0.56 M ⊙ and a distance of d = 48.1 ± 2.2 kpc (Orosz et al. 2014), all observations fall within the acceptable Eddington limit regime.Figure 9 shows the spectra from the six observations analyzed, together with the residuals produced when fitting with our baseline model and with the best-performing reflection models.We note that while our baseline model produces acceptable fits and the residuals do not indicate clear reflection features, for the six observations presented the fit statistic improves significantly when including relativistic reflection to the models.
The 1D posterior distributions obtained by running the MCMC analysis on the fits to the spectra from the six observations (shown in Figure 10) show a general agreement for a high measured BH spin.However, the quality of the spectra prohibits a consentaneous inclination measurement.Using our combining algorithm, we obtain a spin of a = 0.928 +0.058 −0.146 and a poorly constrained inclination of θ = 38 +14 −13 [ • ].Interestingly, the measured spin is in strong disagreement with results obtained using continuum fitting, which indicate a spin a ∼ 0.25 (Steiner et al. 2014;Jana et al. 2021a;Svoboda et al. 2024;Majumder et al. 2024).However, it is important to note that Zdziarski et al. (2024) showed that spin measurements using continuum fitting for LMC X-3 and Cygnus X-1 can be heavily model dependent, and rely on assumptions regarding the accretion disk.Furthermore, Svoboda et al. (2024) placed an upper limit on the spin of the source to be a ≤ 0.7 using an X-ray polarime-try study.This disagreement between the methods is of great concern, and should be used as a starting point for identifying potential issues with the methods used by these BH spin measurement techniques.Ultimately, a spectropolarimetric analysis of the source, in which continuum fitting is combined with a relativistic reflection analysis, should lead to consistent results of the BH spin.

LMC X-1
The BH in HMXB LMC X-1 has been well studied, and previous results strongly suggest that it is rapidly rotating.Steiner et al. (2012) used relativistic reflection on Suzaku and RXTE data to measure the spin of the BH in LMC X-1, finding a = 0.97 +0.02 −0.25 .Jana et al. (2021a) used three NuSTAR observations to measure a BH spin of a = 0.935±0.015using relativistic reflection, result confirmed by Bhuvana et al. (2022), who measured a = 0.93 ± 0.01.Continuum fitting measurements find agreeing results, with Gou et al. (2009) reporting a = 0.92 +0.05 −0.07 based on RXTE data, result confirmed by Mudambi et al. (2020) who measure a = 0.93 +0.04 −0.03 .LMC X-1 was observed four times by NuSTAR, under ObsIDs 30001039002, 30001143002, 30201029002, and 90801324002.Using the mass and distance measurements of M = 10.91 ± 1.41 M ⊙ and d = 48.1 ± 2.2 kpc (Orosz et al. 2009), we find that all four observations occurred while the source had an Eddington fraction 0.1% ≤ L/L Edd ≤ 30%.We fit the spectra from the four observations in the 3-25, 3-30, 3-30, and 3-20 keV bands, respectively.The top panels in Figure 11 show the unfolded spectra obtained from the four NuSTAR observations.The middle panels in Figure 11 show the residuals produced when fitting the NuSTAR spectra of LMC X-3 using TBabs*(diskbb+powerlaw), indicating clear signs of relativistic reflection.
The bottom panels in Figure 11 show the residuals produced when fitting the spectra with the bestperforming reflection models for each observation.For three of the four observations, the preferred relxill variant is the one that allows probing higher accretion disk densities (relxillD), with the spectra from ObsIDs 30001143002 and 90801324002 preferring log(N ) = 19 and the spectra from ObsID 30201029002 preferring the relxillD variant with log(N ) = 17.For remaining observation (ObsID 30001039002), the preferred reflection model statistically favors the addition of a narrow Gaussian absorption component around 6.67±0.04keV, likely consistent with a slightly blueshifted neutral Fe Kα absorption caused by the stellar wind launched by the massive companion star in the system.The individual posterior distributions for the BH spin and inclina-tion of the inner disk are shown in Figure 12, with the thickness of the lines being proportional to the ratios of the reflected to the total flux in the 3-79 keV band, which were used as weights when combining the individual measurements.The combined distributions are represented in Figure 12 through the solid black curves, while the vertical solid and dashed lines indicate the mode and 1σ credible intervals of the distributions, respectively.We find a spin of a = 0.897 +0.077 −0.176 , formally consistent with previous results.The inclination was allowed to vary freely in our analysis, and we measured a value of 50 +10 −13 [ • ].This value is slightly higher but formally consistent with the binary inclination of the system, θ = 36.38• ± 1.92 • (Orosz et al. 2009), suggesting that the inner accretion disk is likely close to being aligned with the binary orbit.

MAXI J1348-630
There are nine archival NuSTAR observations of MAXI J1348-630 (ObsIDs 80402315002, 80402315004, 80402315006, 80402315008, 80402315010, 80402315012, 80502304002, 80502304004, and 80502304006).All nine observations were previously analyzed in Jia et al. (2022), who obtained a spin of a = 0.78 ± 0.04 and an inclination of θ = 29.2+0.3 −0.5 [ • ].However, in their analysis, Jia et al. (2022) fixed the inner and outer emissivity indices q 1 = q 2 = 3, which suggests that this measurement should be interpreted as a lower limit on the spin (Fabian et al. 2014).Additionally, Jia et al. (2022) did not test for the possibility of the presence of narrow absorption features in the spectra.In our analysis, we find that narrow absorption gaussian features are statistically preferred in seven out of the nine observations.Using a mass of M = 11 ± 2 M ⊙ and a distance of d = 3.4 ± 0.4kpc (Lamer et al. 2021), we find that all observations are within the assumed constraints of Eddington range for an observation to be used in our analysis.
Panels (a) in Figure 13 show the residuals produced when fitting the NuSTAR spectra from the nine observations with a power-law model, together with the statistic produced.panels (b) show the residuals produced when fitting the spectra with our best-performing models.The individual 1D posterior probability distributions for the BH spin and the viewing inclination are shown in the left and right panels of Figure 14, respectively, using the colored curves.The combined distributions are shown through the solid black curves in Figure 14.Through our analysis, we measure a = 0.977 +0.017 −0.055 and θ = 52 +8 −11 [ • ].The BH spin takes values higher than those measured by Jia et al. (2022), as expected when allowing the emissivity parameters to vary freely.
Through the introduction of the gaussian absorption features, the inclination generally takes higher values, leading to an inclination measurement higher than that reported by Jia et al. (2022).
It is worth focusing on the results from a few observations.For ObsID 80402315012, the fits do not require the addition of an absorption feature, with the statistical improvement being insignificant.This observation has a high reflection strength, leading to the results produced by the fits to the spectra from this observation being weighted more strongly when combining independent results.The fits to the spectra from ObsID 80402315012 produce spin and inclination constraints lower than the combined distributions, peaking around a ∼ 0.84 and θ ∼ 24 • , respectively.However, when including a gaussian component around 7 keV, both the spin and the inclination distributions peak at higher values, more consistent with the other measurements: a ∼ 0.94 and θ ∼ 39 • .Nevertheless, as the improvement in terms of χ 2 was minimal and the deviance information criterion (DIC) computed based on the MCMC runs was worse due to the increased complexity of the model when including the gaussian component, we chose report the results of the analysis which did not include it, even if the inclusion would lead to results in better agreement with the ones from the other observations.In contrast, the fits to the spectra from ObsID 80502304006 produce a low inclination despite the inclusion of the absorption feature in the model.At the same time, the fits to the spectra from ObsID 80402315004 produce a high inclination measurement despite not including the absorption feature around 7 keV.At the same time, the MCMC analysis of the spectra from ObsID 80502304002 produces a posterior distribution for spin that has two peaks, one around a = 0.77 and one for nearly maximal values of a.This is likely due to the inability of the data to distinguish between similarly good solutions.Despite the best fit in terms of χ 2 being achieved for a nearly maximal spin, the solution with lower spin produces a similar statistic.As the parameter space is wider around the lower spin solution, the walkers converge to that solution and have trouble returning to the best-fit solution (with high spin) due to the parameter space being narrower.A similar phenomenon was described in Draghis et al. (2023a).All these systematic variations are encapsulated in the uncertainties reported on our measurements.

GS 1354-645
Two of the existing three NuSTAR observations of GS 1354-645 were analyzed by El-Batal et al. (2016), where a BH spin of a ≥ 0.98 and a high inclination of θ = 75 • ± 2 • were reported.We analyzed all three existing observations (ObsIDs 90101006002, 90101006004, and 90101006006), and fit the spectra in the entire 3-79 keV NuSTAR band.We note that the BH mass and distance to the system are highly uncertain.We used M = 7.6 ± 0.7 M ⊙ and d = 43 ± 18 kpc (Casares et al. 2009).
The top panels in Figure 15 show the unfolded spectra of the three observations, and the middle panels show the residuals produced when fitting the spectra with models that do not account for relativistic reflection.These residuals show clear features of relativistic reflection.The bottom panels of Figure 15 show the residuals produced by the best-fitting reflection models, together with the fit statistic.
The posterior distributions produced by the MCMC analysis for the spin and inclination parameters are shown in the left and right panels of Figure 16, with the width of the lines being proportional to the weighting used when combining the measurements.In the two panels of Figure 16, the black dotted curves show the combined distribution obtained using the method highlighted in Draghis et al. (2023b), and the black solid curves show the probability distribution of the combined measurements obtained using the method presented in Draghis et al. (2023a).We measured a spin of a = 0.849 +0.103  −0.221 and an inclination of θ = 47 +11 −10 [ • ].It is interesting to note that the inclination measurement is lower and inconsistent with that of El-Batal et al. (2016), and that the spin distribution is driven to lower values by the measurements of two of the three available observations.However, the third observation (ObsID 90101006006) produces a posterior distribution that is strongly peaked at very high values and it does produce a measurement consistent with that previous results.When combining the independent measurements from the three observations in a way that more strongly accounts for systematic uncertainties between observations, we obtain a more loose constraint on the spin, which peaks at a lower value than previously reported.

MAXI J1535-571
We fit 16 of the 17 available NuSTAR observations of MAXI J1535-571, excluding the observation during which the source flux was placing the observation outside the Eddington ratio range for which we expect the accretion disk to extend to the ISCO.The fluxes were calculated assuming a BH mass of M = 8.9 ± 1 M ⊙ (Shang et al. 2019) and a distance of d = 4.1 ± 0.6kpc (Chauhan et al. 2019).The residuals produced when fitting the spectra with our baseline models are shown in the top panels of Figure 17, and the residuals produced when fitting the spectra with the best-performing re-flection models are show in the bottom panels in Figure 17.
Similar to the case of a few other sources, the 1D posterior distributions obtained from the MCMC runs on the best-performing reflection models indicate that the inclination of the inner disk is not consistently measured (shown in the right panel of Figure 18).The combining algorithm produces a poorly constrained inclination measurement of θ = 44 +17 −19 [ • ].The large uncertainty on the measurement is caused by some independent measurements preferring very low inclination values, while others preferring nearly maximal inclinations.The fits for a handful of spectra statistically prefer the addition of gaussian absorption component around 7 keV, consistent with absorption caused by ionized winds, which would suggest a high viewing inclination of the source.However, there is no apparent trend that would suggest that the addition of the absorption component leads to preferentially low or high inclination measurements.This result suggests the need to better quantify the ability of the reflection models to constrain the viewing inclination.The spin measurements are all consistent with a high value, with the individual 1D posterior distributions being shown in the left panel of Figure 18.However, a few measurements suggest a BH spin around a ∼ 0.9, while most prefer nearly maximal values.Interestingly, the same observations that seem to prefer a somewhat lower spin also prefer a low inclination, such as the ones produced by fits to the spectra from ObsIDs 80302309006 and 80402302008.Using the individual 1D posterior distributions, our combining algorithm produces a spin measurement of a = 0.979 +0.015 −0.049 , in agreement with the measurement of Dong et al. (2022), who report a spin of a = 0.985 +0.002 −0.004 and an inclination in the range 70 • ≤ θ ≤ 74 • .Interestingly, our inclination measurement is in reasonable agreement with the inclination determined using ballistic motion in the observed radio jet of θ ≤ 45 • (Russell et al. 2019).

MAXI J1631-479
We analyzed the four existing NuSTAR observations of MAXI J1631-479 (ObsIDs 80401316002, 80401316004, 90401372002, and 90501301001).Due to its proximity to the Sun during ObsID 90501301001, the source tracking was unstable at times, causing the position of the source on the detectors to drift.During one such episode of imperfect tracking, the source landed on detector (DET) 1 in FPMA (nominally the source was placed on DET 0), and on the gap between DET 0 and DET 1 in FPMB.The source being located on the gap between the detectors in FPMB for part of the observation led to improper source reconstruction while running the default nupipeline procedures, causing a significant difference between the spectra measured by the FPMA and FPMB NuSTAR detectors.Therefore, we split the Observing Mode 01 (SCIENCE) files produced during the observation based on the camera head unit (CHU) configuration used to reconstruct the science image, similarly to the procedure outlined in the suggestions regarding Mode 06 observations.Using the event files generated while using different CHU combinations, we identified the interval during the observation during which the source drifted into the gap between DET 0 and DET 1 on FPMB, and created user good-time intervals (GTIs) to exclude this segment of the observation, which were used to reextract Mode 01 science products from both the FPMA and FPMB detectors, leading to spectra in significantly better agreement.The spectra from the four observations, together with the residuals produced when fitting them with our baseline model and with the best-fitting reflection models, are shown in Figure 19.
Upon running the MCMC analysis of the best-fitting spectral models for the four individual observations, we obtained the 1D posterior probability distributions shown in Figure 20.We combined the independent measurements using our combining algorithm, and measured a spin of a = 0.951 +0.039 −0.077 and an inclination of It is important to note that the independent probability distributions for inclination produced by the four observations are in disagreement, with the inclination measured using the spectra from ObsID 90501301001 taking a value of θ ∼ 75 • .The discrepancy in inclination measurements likely comes from an inability of the data to constrain the presence and importance of absorption features around ∼ 7 keV, which land on the blue wing of the Fe K line complex, strongly influencing the ability to measure inclination.However, our measurements for both the BH spin and inclination are in good agreement with those of Xu et al. (2020), who measured a spin a ≥ 0.94 and an inclination of θ = 29 • ± 1 • .

4U 1630-472
At the time of the analysis, six NuSTAR observations of 4U 1630-472 were publicly available (ObsIDs 40014009001, 80802313002, 80802313004, 80802313006, 90002004002, and 90002004004).Using the spectra from ObsID 40014009001, King et al. (2014) measured a BH spin of a = 0.985 +0.005 −0.014 , and highlighted the presence of absorption features linked to ionized outflows in the source.We detected similar absorption features in all six observations analyzed.As the goal of the analysis is to describe the shape of the relativistically broad-ened reflection features and to measure the BH spin, we adopted the simplest models which properly account for the absorption features in the spectrum.Therefore, in all models, we added a gaussian absorption feature and included the zxipcf multiplicative component.The unfolded spectra from the six NuSTAR observations are shown in Figure 21, together with the residuals produced using our baseline model, and the best-performing reflection models.
Upon running the MCMC analysis of the bestperforming models, we obtained the 1D posterior probability distributions shown in Figure 22 for spin (left) and inclination (right).We combined the measurements using our combining algorithm and obtained a spin of a = 0.857 +0.095 −0.211 and an inclination of θ = 55 +8 −11 [ • ].The inclination and the spin measured are formally consistent with those measured by King et al. (2014), however taking lower values, and with larger uncertainties.This is a product of analyzing a larger number of observations, which independently lead to varying results.It is important to note that while all observations show a preference for high spin values, the measurement is strongly determined by the results of the analysis on the spectra from ObsIDs 90002004004 and 80802313006.
3.9.Swift J1658.2-4242Xu et al. (2018) used NuSTAR observation 90401307002 to place a limit on the spin of Swift J1658.2-4242 to be a > 0.94, and measured the inclination of the inner disk to be θ = 64 +2 −3 [ • ].In our analysis, we used the eight available NuSTAR observations: ObsIDs 80301301002, 80302302002, 80302302004, 80302302006, 80302302008, 80302302010, 90401307002, and 90401317002.In all observations, a clear absorption feature was detected around ∼ 7.1 keV, so we included a gaussian absorption feature in all the models that we tested.We note that this source shows long-and short-scale dips in the light curves.For consistency with the rest of the analysis, we extracted and fit the time-averaged spectra from all NuSTAR observations.The top panels in Figure 23 show the residuals produced when fitting the spectra from the eight observations with the TBabs*(diskbb+powerlaw) model.The bottom panels show the residuals in terms of σ produced when fitting the spectra with the best-performing reflection models, which include the gaussian absorption feature around ∼ 7.1 keV.No strong residuals are apparent in the lower panels, suggesting that good fits were achieved for all spectra.
The individual 1D posterior distributions for spin and inclination are shown in the left and right panels of Figure 24, respectively.Using our combining algorithm, we obtain a = 0.951 +0.031  −0.069 and θ = 50 +9 −10 [ • ].We note that for one single observation, the inclination measurement is systematically different from the rest.If we fix the inclination in the fits to the spectra from Ob-sID 90401317002 to θ = 50 • , we obtain a fit worse by ∆χ 2 = 36, but the spin measurement is unaffected.This again highlights the importance of analyzing all available data in order to encapsulate and characterize systematic uncertainties due to variations between observations.

GX 339-4
Of the 38 available archival NuSTAR observations, two occurred while the source was at an Eddington fraction below 10 −3 .We computed this using the values estimated by Parker et al. (2016): M = 9 ± 1.5 M ⊙ and d = 8.4 ± 0.9 kpc.We fit the spectra from the other 36 observations, with ObsIDs listed in Appendix A. We show the residuals obtained through fits using our baseline model and the best-performing reflection models in Figure 25.
The 1D posterior distributions for BH spin and inclination are shown in the left and right panels of Figure 26, respectively.We combined the measurements using our combining method, and obtained a spin of a = 0.970 +0.026 −0.076 , in good agreement with the result of Parker et al. (2016) who find a spin of a = 0.95 +0.02 −0.08 .However, the independent 1D posteriors for inclination range throughout the entire parameter space, with our combining algorithm measuring a poorly constrained inclination of θ = 49 ± 14 [ • ].It is unclear what aspects of the spectra or the models lead to the inability to constrain the inclination of the source reliably between observations, and further investigation is required to understand this discrepancy.

IGR J17091-3624
We analyzed the 14 public NuSTAR observations available at the time.We assumed a mass of M = 12.15 ± 3.45 M ⊙ and a distance of d = 12.6 ± 2.0 kpc (Iyer et al. 2015).A few observations produced fluxes in slight excess of the upper limit of Eddington fraction required by our analysis.However, as the flux estimates are model dependent and as the uncertainties of the BH properties are relatively large, we included all observations in the analysis.The residuals produced when fitting the spectra from the 14 observations with our baseline model are shown in the top panels in Figure 27, and the residuals produced by the best-performing reflection models together with the statistic produced are shown in the bottom panels in Figure 27.
Similarly to the case of LMC X-3 (described in Subsection 3.2), the individual spin measurements agree on a high value, and the inclination measurements take values across the entire parameter space.The 1D posterior distributions obtained through the MCMC analysis of the best-performing reflection models are shown in Figure 28.Using our measurement combining algorithm, we obtain a = 0.963 +0.027 −0.085 and θ = 47 +10 −11 [ • ].The high spin measurement for this source is in disagreement with that measured by Wang et al. (2018b), who report a spin in the range −0.13 ≤ a ≤ 0.27.It is worth noting that in their analysis, Wang et al. (2018b) reported that the source favors a low inner emissivity index, q 1 ∼ 3.7, which could potentially lead to an underestimation of the BH spin (see Draghis et al. 2023a and Section 3.13).This source is particularly interesting, as it exhibits 'heartbeat' variability (Altamirano et al. 2011) similar to that observed in GRS 1915+105 (see, e.g.Belloni et al. 1997;Neilsen et al. 2011), for which studies have shown that the BH spin is high.

GRS 1716-249
NuSTAR observed GRS 1716-249 six times during its 2017 outburst.However, ObsID 80201034004 had an effective exposure of only 408 s, and the spectra are well fit using a power-law model.Therefore, we only analyzed the other five observations, with ObsIDs 80201034006, 80201034007, 90202055002, 90202055004, and 90301007002.All observations fall within the required Eddington range regime when assuming a mass of M = 6.5 ± 1.5 M ⊙ (Tao et al. 2019) and a distance of d = 2.4 ± 0.4 kpc (della Valle et al. 1994).The residuals produced when fitting the NuSTAR FPMA and FPMB spectra using a power-law model are shown in the panels (b) in Figure 29, clearly indicating signs of relativistic reflection.Tao et al. (2019) analyzed the last three observations mentioned above, jointly with simultaneous Swift observations.By only fitting the NuSTAR data in the 4.5-79 keV band, accounting for the thermal emission from the accretion disk using the kerrbb model, and by using only the default relxill model flavor while fixing the inner and outer emissivity indices q 1 = q 2 = 3, they constrain the spin to a > 0.92 and the inclination θ to 40 − 50 • .
We used our pipeline on all the existing NuSTAR data in the entire 3-79 keV band pass.The best-fit models that account for reflection features produce the residuals shown in panels (c) in Figure 29, while panels (a) show the unfolded NuSTAR spectra during the observations.The individual posterior probability distributions for BH spin and inner disk inclination produced by the MCMC analysis on the best-performing models are shown in the left and right panels of Figure 30, respectively.Combining the individual measurements produces a spin measurement of a = 0.97 +0.02 −0.06 and an inclination measurement of θ = 59 +7 −12 [ • ].

GRS 1739-278
At the time of the analysis, there were five public NuSTAR observations of GRS 1739-278.ObsID 80002018002 was taken in 2014, and is the only one of the five during which the source had a flux high enough to fall within the Eddington ratio limits for which we expect the accretion disk to extend to the ISCO.We computed this using estimated values for the BH mass, M = 6.75 ± 2.75 M ⊙ , and distance, d = 7.25 ± 1.25 kpc, in agreement with the findings of Wang et al. (2018a).During the other four archival observations, taken in 2015 and 2016, the source had a significantly lower flux.Additionally, we obtained a NuSTAR DDT observation of GRS 1739-278 during its 2023 outburst (Ob-sID 90901323002), during which the source was in a soft spectral state, and again within the preferred Eddington ratio regime.
For ObsID 80002018002, we extracted spectra from circular source regions with a radius of 90", in order to exclude the contribution from the nearby source 2CXO J174231.5-274349.The FPMA observation was contaminated by stay light, so we manually placed a circular background region of the same size in an uncontaminated region, and located the background extraction region for FPMB in a similar location to that in FPMA.For ObsID 90901323002, 2CXO J174231.5-274349does not show any activity, so we used our default sized 120" circular source region for the source.Similarly to the previous observation, the FMPB observation is contaminated by stray light, so we placed background regions of the same size as the source regions in appropriate locations to account for the effect.
The spectra extracted are shown in the top panels of Figure 31.During the 2014 observation, GRS 1739-278 was in a hard state, while the 2023 observation captured the source in a soft state.The central panels show the residuals produced when fitting the spectra with the TBabs*(diskbb+powerlaw) model, indicating clear signs of relativistic reflection.The bottom panels of Figure 31 show the residuals produced by our bestperforming reflection models.
Two aspects of the fits are important to note.The spectra obtained from ObsID 80002018002 do not statistically require the addition of a diskbb component, so our best fit is achieved using the TBabs*relxill model.In the fits to the spectra from ObsID 90901323002, the Galactic column density N H goes to zero, indicating either that during the first observation, the measured column density is mostly due to intrinsic absorption in the system, or that the correlation between the TBabs and the diskbb components in the model cannot be broken given the low energy bound of the NuSTAR spectra of 3 keV.All the fits to the spectra from the two observations prefer high spins and inclinations.Using the spectra from ObsID 80002018002, Miller et al. (2015) measured a high spin, a = 0.8 ± 0.2 and an intermediate inclination The posterior distributions for the spin and inclination obtained by running the MCMC analysis on the bestperforming models are shown through the blue and red curves in the left and right panels of Figure 32.Given the quality of the data from ObsID 90901323002, despite the best fit being achieved for a high spin, a large number of the posterior samples prefer a solution with similar but worse statistic in which the spin is negative and the inner emissivity index q 1 is low.This behavior is frequent when the quality of the data is not good enough for the fit to strongly distinguish between two similar families of solutions.This behavior is described in Draghis et al. (2023a), and the measured low spin is to be interpreted as a lower limit due to the measured low inner emissivity index.Nevertheless, when combining the posterior distributions using our combining algorithm, we measure a spin of a = 0.968 +0.022 −0.074 and an inclination The spin is formally in agreement with the value measured by Miller et al. (2015), but the measured inclination takes a higher value, consistent with the values measured by the fits to the two individual sets of NuSTAR spectra.Visually inspecting the left panel of Figure 32 indicates that the combining algorithm accounts for the contribution of the second independent measurement, which produces a low spin by bringing the lower bound of our uncertainty interval to lower values, suggesting that the method well includes the systematic uncertainties of the variation of the measurements caused by differences between observations.

1E 1740.7-2942
At the time of the analysis, there were four NuSTAR observations of 1E 1740.7-2942available (ObsIDs 10002012001, 10002021001, 10002021003, 90701317002).Stecchini et al. (2020) analyzed the FPMA spectrum from ObsID 10002012001 in addition to XMM-Newton and INTEGRAL observations of the source, and using the relativistic reflection method they find a nearly maximal but poorly constrained BH spin, quoting a ≥ 0.5.Additionally, Stecchini et al. (2020) estimated the BH mass to be around M ∼ 5 M ⊙ and the distance to the system to be comparable to that to the Galactic center.Therefore, for estimating the Eddington fractions, we assumed a BH mass of M = 5±0.5 M ⊙ and a distance of d = 8.5±2 kpc.All four available NuS-TAR observations of the source happened while, given the above assumptions, the source was at an Eddington fraction within the range required for this analysis.
We extracted the spectra from the four NuSTAR observations and fit them with our baseline family of models.The reflection models perform significantly better than the simple TBabs*(diskbb+powerlaw) model, but similarly to models that describe the continuum as a thermally Comptonized continuum: TBabs*(diskbb+nthcomp).This suggests that reflection is not statistically significant in the four NuS-TAR observations.However, for the same observation as the one treated by Stecchini et al. (2020) (Ob-sID 10002012001), the TBabs*(diskbb+relxill) model does reach a solution that produces a statistically significant improvement through reflection, but some of the parameters are unphysical.The required accretion disk temperature is high, T in ∼ 2.5 keV, with the normalization for the FPMB spectrum going to zero.In other words, the FPMA spectrum drives the diskbb component to a high temperature to account for a difference between the two NuSTAR FPM spectra.To test for a possible difference between the detectors artificially induced by our spectral extraction, we reextracted the spectra from this observation, this time manually placing circular background extraction regions instead of the annular background regions.The resulted spectra produce the same fits, suggesting that there is indeed a difference between the two NuSTAR spectra.When fitting the two spectra simultaneously and allowing the normalization of the relxill component and the power-law indices to differ for the two spectra, the quality of the fit is significantly improved and the diskbb component becomes no longer required.
Inspecting the light curves of the observations shows that toward the end of the exposure, the count rate in the FPMB detector drops significantly, while the rate in the FPMA detector remained unchanged.We identified the source of this abnormal behavior to be a shift of the position of the source to the gap between two of the four detectors that make up the NuSTAR FPMs.We generated GTI files that excluded the last part of the observation, during which the decrease in FPMB count rate occurred.We reextracted the spectra using the GTI files, and fit them again with the same models.
If we add reflection to the models, we obtain χ 2 /ν = 365.21/335when fitting with TBabs*(diskbb+relxill).Further improvement is achieved if the power-law index is allowed to vary between the FPMA and FPMB spectra in this model, obtaining χ 2 /ν = 356.22/337,but this model no longer requires the presence of a diskbb component, making the model TBabs*relxill.Similarly, when fitting the spectra with TBabs*relxillCp with free power-law indices, the fit returns χ 2 /ν = 356.25/337.Note that given the high column density along the line of sight, the low energies are heavily obscured in the observed spectra.Therefore, it is unsurprising that the disk component is weakly required by the NuSTAR spectra.However, when the power-law indices are linked between the spectra from the two NuSTAR detectors, removing the diskbb component produces a significantly worse fit, with χ 2 /ν = 378.07/338.Therefore, we ran the MCMC analysis for the mentioned best-performing reflection models.
The top panel in Figure 33 shows the unfolded spectra from ObsID 10002012001.
The middle and bottom panels in Figure 33 show the residuals in terms of σ produced when fitting the spectra with the TBabs*(diskbb+powerlaw) and the TBabs*(diskbb+relxill) models, together with the statistic produced.In this figure, we show the residuals produced by the reflection model in which the power-law indices Γ are linked between the FPMA and FPMB spectra.Despite the models with free Γ being slightly statistically favored in terms of DIC following the MCMC analysis, they both measure low inner disk inclinations, which are unlikely given the clearly observed radio jets in this source (Luque-Escamilla et al. 2015).Furthermore, the broadband fits performed by Stecchini et al. (2020), which include XMM-Newton and Integral spectra, show the presence of a thermal disk component in the system.For these reasons, and in order to maintain consistency with the rest of the analyzed sample, we report the results of the TBabs*(diskbb+relxill) model with Γ linked between the two FPM spectra, despite it being slightly statistically disfavored.However, as seen in Figure 34, the spin measurements are generally consistent between the fits.
The two panels in Figure 34 show the posterior distributions for the spin (left) and inclination (right) for the TBabs*(diskbb+relxill) model with Γ linked (blue line) and for the TBabs*relxill model with Γ free (red line).The dotted and solid black curves show the results produced by our combining algorithm when ran only on the blue posterior distributions.The values measured are a = 0.856 +0.134  −0.443 and θ = 31 +29 −18 [ • ].Note that the inclination of the system is unconstrained by the fit when linking the power-law indices between the spectra, and constrained to low values when allowing them to vary independently.This latter constraint is likely an artifact of an improper characterization of the underlying continuum caused by the lack of inclusion of the diskbb component, which given the quality of the data, is not statistically significant.The modes of the independent distributions for the BH spin are high and much more narrowly constrained, with a = 0.97 +0.3  −0.32 when linking Γ, and a = 0.994 +0.004 −0.059 when allowing Γ to vary independently.The independent measurements likely strongly underestimate the size of the systematic uncertainty, but by applying our combining algorithm even on one fit alone, we obtain a much better estimate of the true size of the uncertainties of the measurement.
Visually inspecting the residuals in Figure 33, the presence of relativistic reflection is not immediately obvious.However, the statistics of the fits indicate that reflection is statistically significant.This raises the valuable point that it is possible that many observations with low signal-to-noise ratios that could still place constraints on the BH spin were simply not analyzed because the residuals did not indicate visually clear reflection features, but the quality of the data might still be good enough to obtain a statistically significant constraint.Lastly, given the shape of the residuals when fitting with TBabs*(diskbb+powerlaw) and TBabs*(diskbb+relxill), another possible interpretation can be offered to the spin constraint: if the spin of the BH was small, the spectra would have shown narrower features, which would have been easier to distinguish from the underlying continuum.Ultimately, we encourage the readers to visit the more in-depth analysis of Stecchini et al. (2020).

Swift J174540.7-290015
Swift J174540.7-290015(hereafter T15) is an XB located near the line of sight to the Galactic center.It was observed to be in an outburst during the NuSTAR observation 90101022002.During the same observation, the source AX J1745.6-2901 was also above quiescence.The left panel in Figure 1 shows the exposure map of the NuSTAR observation.The green circle shows the 50" region used for extraction of the source spectrum.The white dashed circle indicates a 50" region at the position of Swift J174540.2-290037, the other source present in this field, also analyzed in this paper (see Section 3.16), which was in quiescence during this observation.The yellow circle shows a 50" region around the source AX J1745.6-2901.The magenta cross represents the position of the Galactic center.The magenta, red, and cyan circles (labeled Bkg 1, Bkg 2, and Bkg 3, respectively) show three different background regions tested for extraction of background spectra.Regions 1 and 2 were placed around AX J1745.6-2901 at a distance equal to the distance between AX J1745.6-2901 and T15 (our target), with the aim of quantifying not only the contribution from the detector background, but also from AX J1745.6-2901 to our source region.Region 3 was placed away from the sources, and quantifies the detector background.The difference between the background spectra from regions 1 and 2 were minimal, with slightly increased rates in region 1 due to contribution from T15.In order to adopt the most conservative background choice, we continued our analysis with the background spectra from region 1.However, we note that our results are consistent when using the background spectrum obtained from region 2.
T15 was analyzed by Mori et al. (2019), who measured a BH spin of a = 0.94 +0.03 −0.01 and an inclination of θ = 64.2+0.9 Of the models tested on the NuSTAR spectra of T15, the fit using TBabs*(diskbb+relxillD 19+gaussian) perform best in terms of DIC.This model includes a gaussian absorption feature around 7.16 keV.The unfolded spectra of the observation are presented in the top panel of

Swift J174540.2-290037
Similarly to T15 (see Section 3.15), Swift J174540.2-290037(hereafter T37) is located near the Galactic center.The source was detected during the NuSTAR observation 90201026002.The right panel of Figure 1 shows the NuSTAR exposure map during the observation.The cyan and green circles show 30" and 120" regions centered at the position of T37, respectively.The magenta cross shows the position of the Galactic center, and the white and yellow dashed circles show the positions of T15 and AX J1745.6-2901, both in quiescence during the observation.The blue annulus shows the typical background region used throughout our analysis, with an inner radius of 200" and an outer radius of 300".The red and cyan circles respectively show the 120" and 30" circles that were used for extracting background regions.We extracted source and background spectra from the 30" regions, from the 120" source and background regions, and from the 120" source region and the background annulus region.We chose to continue our analysis with the source spectra extracted from the 120" source region and background spectra extracted from the annulus, as this is the most conservative pair of source and background regions, and also ensures consistency with the method used throughout the entire paper.We note that using different pairs of source and background regions produce consistent results.
T37 was analyzed by Mori et al. (2019), who measured a BH spin of a = 0.92 +0.07 −0.05 and an inclination of θ = 21 +2 −3 [ • ].The top panel in Figure 37 shows the unfolded spectra of the NuSTAR observation of T37, while the middle and bottom panels show the residuals in terms of σ produced when fitting the spectra with a model that does not account for reflection, and with our best-performing reflection model, TBabs*(diskbb+relxillCp), respectively.The 1D posterior probability distributions for the spin and inclination are shown through the red curves in the left and right panels of Figure 38, respectively.The dotted and solid black curves show the results of applying the two methods of our combining algorithm on the single posterior distributions.We measure a = 0.774 +0.082 −0.106 and θ = 31 +8 −9 [ • ].

MAXI J1803-298
NuSTAR observed the 2021 outburst of MAXI J1803-298 four times (ObsIDs 80701332002, 90702316002, 90702318002, and 90702318003).We fit the spectra extracted from all four observations (shown in the top panels of Figure 39) in the entire 3-79 keV NuSTAR band pass, with the exception of the spectra from Ob-sID 90702318003, which we only fit in the 3-20 keV band, as they became background dominated at higher energies.Feng et al. (2022a)  The middle and bottom panels show the residuals produced when fitting the spectra with the model that ignores reflection and with our best-performing reflection models, which include narrow absorption features around 7 keV in all cases.Running the MCMC analysis on our best-performing models produces the 1D posterior probability distributions shown in Figure 40

MAXI J1813-095
NuSTAR observed MAXI J1813-095 three times in 2018 (ObsIDs 80402303002, 80402303004, and 80402303006).Jana et al. (2021b) fit the three observations with a model that accounts for relativistic reflection to find a lower limit for the spin parameter of a ≥ 0.76 and an inclination in between 28 • < θ < 45 • .Later, Jiang et al. (2022) used an averaged spectrum obtained from the existing three NuSTAR observations to study the reflection component, and by freezing the spin parameter in their fits to the maximum possible value of a = 0.998, they find an inner disk radius consistent with the ISCO and an inner disk inclination of θ = 20 +7 −10 [ • ].Using the mass (M = 7.4±1.5 M ⊙ ) and distance (d = 6 ± 1 kpc) estimates for the BH in the system obtained by Jana et al. (2021b), we find that all three existing NuSTAR observations occurred while the source was at an Eddington fraction for which we expect the inner disk radius to extend to the ISCO.Therefore, we continue analyzing all three existing observations of MAXI J1813-095 in the entire NuSTAR band pass of 3-79 keV band.The spectra from the three observations together with the residuals produced when fitting the spectra with a model not accounting for reflection are shown in the top and middle panels in Figure 41, indicating clear signs of relativistic reflection.
The residuals produced by fitting the spectra with our best-performing reflection models are shown in the bottom panels of Figure 41, together with the statistic produced.We note that for the third observation (80402303006), the diskbb component was not required by the data, but adding a narrow Gaussian absorption feature around ∼ 6.7 keV produced a statistically significant improvement in the fit.In order to best assess the statistical significance of the improvement in statistic produced by the addition of the narrow absorption component, we removed the diskbb component for simplicity.We tested whether the addition of similar Gaussian features improves the quality of the fit for the other two observations, but the decrease in χ 2 was not statistically significant.
The posterior distributions for the spin and inclination parameters obtained through the MCMC analysis are shown in the left and right panels of Figure 42, respectively.The combined distributions for the two parameters are shown through the solid black curves in Figure 42.We measure a spin of a = 0.88 +0.10 −0.27 and an inclination of θ = 42 +11 −13 [ • ].The three individual spin measurements agree on a high spin measurements, with the third observation producing a significantly narrower spin constraint, which drives the 1σ credible region of the combined distribution to produce a relatively precise spin constraint.The three individual inclination measurements are not in good agreement, hinting perhaps at the fact that the simplicity of the models does not fully capture the complexity of the data.However, this source highlights the benefit of using an averaging function when combining the individual measurements into a single posterior distribution which encapsulates systematic uncertainties not accounted for during the analysis.Our results are in agreement with those of Jana et al. (2021b).

MAXI J1848-015
The two existing NuSTAR observations of MAXI J1848-015 (ObsIDs 90601340002 and 90601341002) have been analyzed by Pike et al. (2022)  For the spectra from ObsID 90601340002, the model that is statistically favored requires the addition of a narrow absorption line around 6.8 keV, while the spectra from ObsID 90601341002 were taken after the source had transitioned to a hard state and do not statistically require the presence of a diskbb component, but they do require the addition of a narrow emission line around 6.3 keV.Additionally, the fits measure a high hydrogen column density in the first observation (N H ∼ 7 cm −2 ) and much lower in the second observation (N H ∼ 2.6 cm −2 ).The change in measured column density and the presence of the narrow emission features were also pointed out and explained in-depth by Pike et al. (2022), so we refrain from physically interpreting these effects.
The 1D posterior distributions for the BH spin and viewing inclination obtained from running the MCMC analysis are shown in Figure 44.By combining the independent posterior distributions using our combining algorithm, we measure a spin of a = 0.753 +0.191  −0.652 and an inclination of θ = 29 +13 −10 [ • ].While the measured inclination agrees with that obtained by Pike et al. (2022), the larger uncertainty better encapsulates the systematic uncertainties of the methods.Formally, our measured spin disagrees with that obtained by Pike et al. (2022).However, visually inspecting the individual posterior probability distributions for the BH spin shown in the left panel of Figure 44, we see that the distribution produced by the MCMC run on the observation taken during the hard state strongly prefers high BH spins, while the soft state observation places much weaker constraints on the spin.Combining the two results while weighting them by the strength of the reflection during the observation produces a relatively loose constraint on the BH spin.This apparent discrepancy should not be interpreted as a clear disagreement, but rather as proof that when encapsulating the variation between the constraining capabilities of multiple independent observations, the uncertainty of this particular measurement is significantly underestimated when not considering systematic effects.

EXO 1846-031
NuSTAR observed EXO 1846-031 six times during its 2019 outburst (ObsIDs 80502303002, 80502303004, 80502303006, 80502303008, 90501334002, and 90501350002).The first observation was analyzed in Draghis et al. (2020), where a spin of a = 0.997 +0.001 −0.002 and an inclination of θ = 73 +2 −1 [ • ] were reported.Draghis et al. (2020) also reported the presence of an absorption feature around 7 keV, which was modeled through the addition of a gaussian component with negative amplitude.Therefore, for the analysis of the six observations of EXO 1846-031, in addition to the models usually tested, we also verified whether the addition of gaussian components was statistically significantly.We find that for three of the six observations, the gaussian component improves the quality of the fit significantly.
The top panels of Figure 45 show the unfolded spectra of the six NuSTAR observations, while the middle supbanels show the residuals produced when fitting the spectra with a power-law model.The bottom panels of Figure 45 show the residuals obtained when fitting the spectra with best-performing models in terms of statistic, together with the χ 2 values produced.The individual posterior distributions for the spin and inclination obtained from the MCMC analysis of the six bestperforming models are shown in the left and right panels of Figure 46 through the colored curves, respectively.The thickness of the lines is proportional to the reflection strength during the observations, which was used as weighting when combining the individual posterior distributions.The combined distributions for the spin and the inclination are shown in Figure 46 through the black solid curves, while the dotted curves show the distribution combined using the method described in Draghis et al. (2023b).The vertical black solid and dashed lines show the mode of the combined distribution, together with the 1σ credible interval.Using the six NuSTAR observations, we measure a spin of a = 0.959 +0.031 −0.077 and an inclination of θ = 62 +9 −10 [ • ].These values are consistent with those reported in Draghis et al. (2020), but in ad-dition to statistical uncertainties, they include the systematic variation originating from the independent measurements on the six different NuSTAR observations.

XTE J1908+094
The two existing NuSTAR observations of XTE J1908+094 were analyzed by Draghis et al. (2021) (Ob-sIDs 90501317002 and 80001014002), who reported a spin of a = 0.55 +0.29  −0.45 and an inclination of θ = 27 +2 −3 [ • ].Draghis et al. (2021) divided the spectra based on the hardness ratios during the observations and fit the four resulting pairs of spectra jointly.We extracted timeaveraged spectra and fit them independently in the 3-50 keV band.The top panels of Figure 47 show the unfolded spectra of the observations.The center panels show the residuals of the models that do not account for reflection, while the lower panels show the residuals produced when accounting for reflection.
We ran the MCMC analysis for the best-performing reflection models, and the posterior distributions for the BH spin and inclination are shown in the left and right panels of Figure 48, respectively.By combining the posterior distributions, we measure a = 0.466 +0.368  −0.517 and θ = 28 ± 11 [ • ].These results are consistent with those found by Draghis et al. (2021), but with larger uncertainties owing to our joint treatment of individual observations.

GRS 1915+105
At the time of the analysis, there were 29 public archival NuSTAR observations of GRS 1915+105.Of the 29 observations, the source was at an Eddingon fraction 10 −3 ≤ L/L Edd ≤ 0.3 only in 18 of them.This was calculated using the mass and distance estimates from Reid et al. (2014): M = 12.4±2 M ⊙ and d = 8.6±2 kpc.Furthermore, due to the complexity of the obscuration of the central engine in this source, we were unable to obtain a good fit using our combination of simplistic models for the spectra from ObsID 30202033004.Therefore, we only fit the remaining 17 observations, listed in Appendix A. Figure 49 shows the residuals produced by fits using our baseline model (top) and the best-performing reflection models (bottom).
Using our combining algorithm on the individual 1D posterior probability distributions inferred using the MCMC analysis on the best-fit models (shown through the colored lines in the left and right panels of Figure 50 for spin and inclination, respectively), we measure a BH spin of a = 0.976 +0.018 −0.056 and an inclination of θ = 60 ± 8 [ • ].Our spin measurement is in good agreement with that of Miller et al. (2013), who performed a reflection study on the NuSTAR spectrum from Ob-sID 10002004001 to measure a = 0.98 ± 0.01.However, our inclination measurement is lower than the value of θ = 72 ± 1 [ • ] determined by Miller et al. (2013).Interestingly, looking at the right panel in Figure 50, we can see that the inclination that we measured from the same observation does indeed take larger values, consistent with that measured by Miller et al. (2013).This result suggests the importance of analyzing all the available observations in order to ensure minimal bias in our measurements originating from physical effects for which our models do not account.

Cygnus X-1
We analyzed all the 34 public NuSTAR observations of Cygnus X-1, as they all occur while the source was at an Eddington fraction for which it is expected that the accretion disk should extend to the ISCO.For this calculation, we used M = 21.2±2.2M⊙ and d = 2.22±0.18kpc (Miller-Jones et al. 2021).We show the residuals produced when fitting the spectra from the individual observations in Figure 51, with the top panels indicating the residuals produced by our baseline model, and the bottom panels indicating the residuals produced by the best-performing reflection models.The source shows signs of variable narrow absorption and emission features, which we account for using the zxipcf multiplicative component, and through the addition of gaussian components with positive or negative normalizations, depending on whether the component is meant to account for an emission or absorption feature, respectively.The broad reflection features are well accounted for.
Figure 52 shows the 1D posterior distributions resulting from the MCMC analysis of the independent observations for spin (left) and inclination (right).The majority of the spin posterior distributions agree on a high spin.However, fits to a few observations prefer lower spin values, likely due to an incomplete or improper characterization of the narrow spectral features present.Furthermore, partial obscuration of the inner disk would likely also artificially weigh the emission from the outer regions of the accretion disk more, leading to a lower spin measurement.The origin and shape of the narrow spectral features, together with their influence on the spin measurement will be properly quantified using spectra from the Resolve instrument on board the recently launched XRISM mission.Our combining algorithm finds a spin of a = 0.95 +0.04 −0.08 .Our measurement is consistent with archival spin measurements both through the relativistic reflection method on NuSTAR (see, e.g.Parker et al. 2015 who find a ≥ 0.97) and through measurements made using continuum fitting (see, e.g.Gou et al. 2014

V404 Cygnus
Given the mass estimate of M = 12 ± 3 M ⊙ (Shahbaz et al. 1994) and the distance of d = 2.39 ± 0.14 kpc (Miller-Jones et al. 2009), only two of the 11 available NuSTAR observations of V404 Cygnus happened while the source was at an Eddington ratio above 10 −3 : Ob-sIDs 90102007002 and 90102007003.The two observations were taken in succession and show significant flaring of the source, with count rates increasing by a factor of ∼ 5 − 10.For both observations, we split the light curves into "flare" and "nonflare" segments, and extracted spectral products for both states.We fit the resulting spectra independently.Figure 53 shows the flare (left) and nonflare (right) spectra extracted from ObsIDs 90102007002 (top) and 90102007003 (bottom), together with the residuals produced when being fit with our baseline model and with the best-performing reflection models.The source was in a hard spectral state, and the preferred reflection models did not require the presence of a diskbb component, which was subsequently removed in order to reduce the size of the parameter space.All spectra required strong ionized obscuration, which we accounted for using zxipcf.Despite the best-fit statistics being poor, especially for the nonflaring state, we note that the residuals produced by our best-performing models indicate that the reflection features have well been accounted for.At this stage, the remaining residuals come from calibration uncertainties in the NuSTAR band, which under normal circumstances are negligible, but given the high signal-to-noise ratios of the observations are strongly accentuated.Furthermore, King et al. (2015) found that the source showed a multitude of narrow spectral emission features which evolved on a timescale of a few kiloseconds.Our models do not account for any such features, as they are not resolvable given the NuSTAR energy resolution.Therefore, despite the formally poor statistic, we continue our analysis under the assumption that the relativistically broadened spectral features are properly accounted for.
Figure 54 shows the individual 1D posterior distributions for spin (left) and inclination (right) produced by the MCMC analysis of the individual fits.The spin is generally high, and the inclination measurements are in broad agreement.Our combining algorithm finds a = 0.935 +0.037 −0.075 and θ = 37 +9 −8 [ • ].The spin measurement is consistent with that of Walton et al. (2017), who report a ≥ 0.92, but our inclination measurement is lower than their reported value of θ ∼ 52 • .
Our models include significant obscuration of the central engine in this source, and both our results and the conclusion of Walton et al. (2017) suggest that the flaring episodes that this source experienced were in fact caused by a reduction in obscuration along the line of sight.If that is the case, it is likely that our spin inference should be treated as a lower limit only, as much of the emission from the innermost regions of the disk would have been attenuated, and larger contribution from outer regions of the disk would weigh more in the observed spectra, artificially biasing the measurement to lower values.

DISCUSSION
In this paper, we reanalyzed all the existing archival NuSTAR observations of sources that had previous spin measurement using the relativistic reflection method.This analysis was performed in a way fully consistent with that in Draghis et al. (2023b), adopting the same set of theoretical assumptions.In conjunction with the analysis in Draghis et al. (2023b), we have compiled a sample of spin measurements of considerable size.
For each source analyzed, we show the individual posterior probability distributions for the spin and inclination, obtained from the MCMC analysis of the spectra from individual observations.We report the values obtained by using a Bayesian combining algorithm to combine the individual measurements into a single result.This resulting combined distribution encapsulates the ignorance of the simplistic treatment of our models, which stems from using a single model to fit observations across various spectral states, at different accretion rates, when likely many more physical processes are present and relevant.Generally, the values that we report in this paper for the spin are in agreement with previous measurements, but the sizes of the uncertainties that we report are significantly larger.
While the spin measurements generally agree between observations, the inclination measurements often tend to disagree, especially in sources where there are many observations available.This is likely due to difficulties in constraining the blue wing of the relativistically broadened Fe K line, owing to uncertainties in constraining the underlying continuum and possible ionized outflows.This suggests that the inclination values obtained through relativistic reflection measurements from a single observation should be taken with caution, and not be treated as definitive.However, when the analysis is per-formed systematically on a large number of observations, the underlying true value can be better estimated.It is important to note that for a few sources, timing analysis has suggested the presence of the Lense-Thirring effect (see, e.g., Motta & Belloni 2023;Zhu et al. 2023), which could partly explain a change in the inferred viewing inclination of the inner disk.
To further probe the hypothesis regarding the effect of the ionized outflows, we ran a case study on ObsID 80402302008 of MAXI J1535-571 (Section 3.6).The best-fit solution obtained predicts a spin around a ∼ 0.89, and inclination of θ ∼ 5 • , with χ 2 r = 581/465.If we fix the spin to the value derived using all observations, a = 0.98, we obtain θ ∼ 7 • , with a slightly worse χ 2 r = 587/466.If instead we fix the inclination to the one derived by averaging all observations, θ = 45 • , we obtain a spin of a ∼ 0.97, but with worse χ 2 r = 643/466.However, adding an absorption gaussian component at 7keV while maintaining the inclination fixed at θ = 45 • maintains the spin high, in good agreement with the combined measurement, but the fit statistic does not represent a significant improvement (χ 2 r = 579/463) over the solution that does not fix the spin and inclination and does not include absorption features (χ 2 r = 581/465).However, the χ 2 value does formally improve, suggesting that in observations with higher signal-to-noise ratios, the absorption feature could have been significant and produce results in agreement with the rest of the observations.As the spin and inclination measurements might be somewhat correlated, and as absorption features from ioinized outflows produce narrow spectral features that overlap with the blue wing of the relativistically broadened Fe K line, improperly characterizing the impact of the ionized absorption could directly impact our ability to accurately constrain BH spins and inclinations.In the future, pairing NuSTAR broadband spectra with high energy resolution XRISM spectra will help to definitively address this issue.
We tested the impact of the ability to measure the viewing inclination on our ability to recover BH spin by reanalyzing all observations of GX 339-4 while fixing the inclination at the value determined through our combining algorithm, and comparing the results with the ones obtained with a free inclination.The top panel in Figure 2 shows the fit statistic produced when fitting the spectra with a free inclination (blue points) and with the inclination fixed at θ = 49 • (red points).Generally, the fits are very similar, or worse.Interestingly, in two cases, the statistic becomes slightly better, indicating that our automatic fitting algorithm might have not Figure 2. Top panel shows the reduced χ 2 obtained when fitting all the spectra of GX 339-4 with inclination free (blue) and fixed at 49 • (red).The central panel shows the spin measured through direct fitting when the inclination is fixed at 49 • (red), free (blue), and measured through the MCMC analysis when the inclination is free (green).Note that for the measurements obtained through χ 2 minimization in xspec, the uncertainties are not included.The bottom panel shows the inclinations measured through χ 2 minimization in xspec (blue) and through MCMC analysis (green).The horizontal red line and red shaded region show, respectively, the mode and ±1σ credible interval on the inclination measurement obtained by combining all the independent measurements.
fully reached a global statistic minimum, but became stuck in a nearby local minimum.
The bottom panel in Figure 2 shows the inclination measurements obtained when allowed to vary freely (in blue) by direct fitting within xspec, and the ones obtained following the MCMC analysis (in empty green circles with the associated error bars).The points showing the values from xspec fitting do not include error bars, as they are often uninformative, and as the errors reported in our work come from the MCMC analysis.The horizontal red line represents the value of 49 • to which the inclination was fixed during this experiment, and the red band shows the ±1σ uncertainty region on the combined measurement.The fact that the inclinations vary significantly despite the quality of the fits not being significantly different when compared to fits using a fixed inclination suggest that one needs to use caution when considering the inclination measurements obtained through this method, and again suggests that mischaracterizing disk winds which would imprint absorption features around 7keV could have a significant impact on our ability to recover the inclination.This hypothesis will be tested in the near future using XRISM data.
The central panel in Figure 2 highlights the effect of fixing the inclination on the recovered spins.The blue points represent the spins inferred from fitting the spectra within xspec when the inclination is allowed to vary freely, while the red points show the spins measured when the inclination is fixed at the value determined by combining all independent measurements.The green circles show the spins determined after running the MCMC analysis, starting from the best-fit values with free inclination (blue points).Generally, the spins measured when the inclination is free agree with those measured when the inclination is fixed at 49 • (note that uncertainties are not shown), suggesting the ability of the data to separate the spin and inclination constraints.The spin generally determines the red wing of the Fe line, while the inclination impacts the blue wing of the line.Nevertheless, through this experiment, we have reconfirmed that the inclination measurements need to be taken with caution, while the spin measurements appear to be robust, especially when combining multiple independent observations.However, a more interesting discrepancy, highlighted by Figure 2, is that obtained by running the MCMC analysis with free inclinations vs. results obtained from xspec χ 2 minimization.Despite the MCMC runs starting with the walkers initialized around the value obtained from χ 2 minimization, the posterior distributions sometimes cluster around a different value.This is likely because of the way in which the parameter space is explored by the walkers.When the quality of the data is reduced, there are often multiple families of solutions with similar statistics.Generally, when that happens, as highlighted previously in the paper, and as explained in Draghis et al. (2023a), solutions with a high spin and a high inner coronal emissivity index produce similar statistics with solutions with a low (even negative) spin and a low inner emissivity index (q 1 ).Due to the nature of the shape of the parameter space, walkers preferentially navigate from the high spin solution (narrow region of the parameter space) to the low spin solution (wider region of the parameter space), and have difficulty navigating back to the narrow parameter space around the high spin solution, despite it being slightly preferred in terms of statistics.This ability to explore the parameter space also explains why when fixing the inner disk radius at the ISCO, we sometimes find low or negative spins, while for the same data, when allowing the inner disk radius to vary, the low spin solutions disappear: the families of solutions become more statistically distinct, and the walkers navigate more easily to the high spin solutions.Additionally, when the inner disk radius is free, the size of the parameter space increases, making it less likely for walkers to wander outside of the starting solution and become stuck in low spin solutions.(See the discussion regarding the effects of allowing the inner disk radius to vary below.) The figures illustrating the posterior distributions obtained from our spectral analysis contain two different "combined" distributions.We report the mode and ±1σ credible interval of the black solid curves, which are obtained by averaging a large sample of beta distributions used to fit the individual measurements to a single distribution.This distribution encapsulates the information from all available measurements from individual observations.In contrast, the black dotted curves show the beta distributions obtained by the most frequently occurring shape of the beta function in the analysis, which combines the individual posterior distributions.While this value could be interpreted as the underlying, "true" spin value from which the measurements are drawn, it does not fully incorporate the size of the uncertainties originating from using incomplete theoretical models to fit the spectra.In other words, the dotted distribution assumes that our models are correct, and the solid curves allow for a more relaxed interpretation of the independent results, leading to larger uncertainties.
In order to ensure consistency between the spins measured in this paper with those in Draghis et al. (2023b), we ran the new version of the combining algorithm, which was established in Draghis et al. (2023a), on the 10 measurements in Draghis et al. (2023b).Additionally, we ran the same combining algorithm on the results from the analysis of single observation of Swift J1728.9-3613(Draghis et al. 2023c).We combined the measurements of BH spin and inclination from this paper with the results from Draghis et al. (2023b), the spin measurement of XTE J2012+381 (Draghis et al. 2023a), and Swift J1728.9-3613(Draghis et al. 2023c) to compile a list of 36 BH spin measurements performed in a systematic way.We show the results in Table 1.We present the original measurements for the BH spin and the inner disk inclination, together with our updated measurements.All our reported values have been obtained through our measurement combining algorithm, even when a single observation of the sources is available.This ensures that the same treatment is applied to

Walton et al. (2017)
Note-The sources are listed in order of increasing Right Ascension.All uncertainties reported are at the 1σ level.The values in the left columns are the measurements obtained using the pipeline highlighted in this paper, and represent the modes and ±1σ credible intervals of the combined posterior probability distributions obtained through the analysis (i.e. the solid black curves in the figures showing the posterior distributions).The values in the right columns show the previous measurements of spin and inclination for the sources, together with their respective references.† XTE J2012+381 does not have a spin measurement prior to Draghis et al. (2023a).‡ The spin of LMC X-3 was measured in Jana et al. 2021a using continuum fitting.* For LMC X-1 and LMC X-3, the inclination values were fixed to the orbital inclination.
all sources, and that the reported uncertainties encapsulate part of the inherent systematic uncertainty that comes with spin measurements made using the relativistic reflection method.
To better illustrate the change induced by our new methods, we plotted all the values in Table 1.The top panels in Figure 3 show the measurements updated using our pipeline versus the values in the literature for the spins (left) and inclinations (right).The bottom-left panel shows the complement to one of the updated and literature measurements, shown on a logarithmic scale.This panel is shown for visualization purpose only, and presents the same information as the top-left panel inverted across the origin.This was done to better highlight changes in the measured values at spins close to the maximal value.Generally, the measurements are consistent with their previous estimates.Three of the four spins with values less than a = 0.75 are now closer to maximal, indicated by their position above the diagonal line in the top-left panel of Figure 3.Many updated spin measurements now take smaller values than previously, but with uncertainties that encapsulate the previous measurements.For the inclinations, the updated measurements are again in good agreement, but with significantly larger uncertainties.The many points below the diagonal line in the top-right panel of Figure 3 indicate a decrease in the measured inclination compared to the literature values.1.
In order to further probe whether any trends are highlighted by our analysis, we plot the change in the measured spins and inclinations versus the literature measurements in the top panels of Figure 4.The lower-left panel in Figure 4 shows the same difference in spin, just focusing on values larger than a = 0.75, where most points are concentrated.Again, for spins, the largest change occurs for previously low values, while originally larger spins are now smaller, but with larger uncertainties.Similar to the previous plot, the trend shown in the top-right panel of Figure 4 suggests that lower inclinations are now measured to be slightly higher, while previously high inclinations are now slightly lower, all with increased uncertainties.Nevertheless, the measurements are still in good agreement with the previously reported values.Lastly, the bottom-right panel in Figure 4 shows the change in spin versus the change in inclination.As no trends are obvious, this suggests that the previous uncertain spin measurements were not generally due to uncertain inclination measurements.The markers represent the same sources as in Figure 3, with all values taken from Table 1.
Throughout our analysis, we only treated observations in which the sources were at an Eddington fraction between 10 −3 ≤ L/L Edd ≤ 0.3, under the assumption that for this range, the accretion disk extends to the ISCO and is not truncated.If applying the assumption that the inner disk extends to the ISCO while it is in fact truncated at larger radii, our measurements based on the particular spectra should be interpreted as lower limits.However, as the calculation of Eddington fraction can, in principle, be biased by erroneous existing BH mass or distance estimates, another observable that could point to the presence of disk truncation is the spectral hardness, with some studies indicating that during the hard state, the accretion disk should be truncated at large radii (see, e.g.Done et al. 2007).To test the robustness of our measurements, we recombined the measurements obtained when analyzing the spectra from four sources that have been observed multiple times, while focusing on observations occurring at high Eddington fractions and/or low hardness.Figure 5 shows the measured spins for all the treated observations of GRS 1915+105, Cygnus X-1, GX 339-4, and MAXI J1820+070 vs. the Eddington fraction during the observations (left, blue) and vs. the hardness ratio during the observations (right, red), defined as the ratio of the flux in the 10-70 keV band to the flux in the 3-5 keV band.The transparency of the points is proportional to the strength of reflection during the observations, which was used as weighting while combining the independent measurements into the reported values.The vertical dashed lines in the left and right sets of panels show the boundaries used for this experiment: we only combined observations for which 10 −2 ≤ L/L Edd ≤ 0.3 and for which Flux 10−79 keV /Flux 3−5 keV ≤ 5.In other words, we only combined the points to the right of the vertical dashed line in the left panels, and to the left of the vertical line in the right panels.While these choices are somewhat arbitrary, they clearly highlight the magnitude by which our measurements would be impacted by restricting the observations treated in this analysis to more explicitly agree with prior expectations of accretion disks extending to the ISCO.
When constraining the Eddington fraction, all observations of Cygnus X-1 fall while the source was below L/L Edd = 10 −2 , making this experiment irrelevant.For MAXI J1820+070, only one observation was cut out, for which the reflection strength was low and the spin was below the previously determined combined value, leading to a nearly unchanged measurement.For the other two sources (GRS 1915+105 and GX 339-4), by removing the observations for which L/L Edd ≤ 10 −2 , we filtered out the main contribution to the combined measurement that provides support for a low or negative spin, leading to a decrease in the uncertainty of our spin constraints.When only considering observations for which Flux 10−79 keV /Flux 3−5 keV ≤ 5, all four sources lose support for low spin, with the combined measurements resulting in narrower uncertainties on a higher BH spin.Still, all the measurements obtained through this experiment are in full agreement with the original measurements which incorporate all available data.
Additionally, we reanalyzed all observations of GX 339-4 to simultaneously allow the inner disk radius and the BH spin to vary.We chose the observations of GX 339-4 for this experiment, as this source spans the broadest range of Eddington fraction of the sources treated in this paper.We note that this treatment does not fix the spin to some value, as similar previous studies do.By allowing the inner disk radius and the spin to vary simultaneously, we properly infer the extent of the variation of the two parameters.If one were to fix the spin to some arbitrarily high value, the inner disk radius would be artificially biased to lower values.A maximally spinning BH with a disk truncated at 6 r g would produce similar gravitational distortions to a nonspinning BH with the disk extending to the ISCO.Therefore, if the spin is erroneously fixed at a high value, for the same spectral shape, the inferred disk truncation would be larger than when allowing the spin to vary freely.
Figure 6 shows the inner disk radius in units of the ISCO radius (top) and the BH spin (bottom) measured by running the MCMC analysis of all the treated observations of GX 339-4.The transparency of the points is proportional to the strength of reflection during the observation.The figure indicates that there is no clear sign of disk truncation throughout the entire range of Eddington fractions treated in our analysis.It is important to acknowledge that this experiment was performed for a single source, and that these conclusions do not necessarily translate directly to all sources.The spin of the BH, the degree of misalignment between the BH and binary system, and the mass and metallicity of the donor star, could cause the critical Eddington fraction to differ between sources.However, even when considering the same source, outbursts are observed to differ significantly; factors including the ratio of coronal and disk flux, mass loss in winds, and the rate of change in the mass accretion rate could cause the critical Eddington fraction to change between outbursts.The fact that no observations of GX 339-4 show strong evidence of truncation above L/L Edd ≥ 0.001 across several outbursts may signal that these considerations are minor.This analysis was performed in order to highlight that our initial assumption, of the accretion disk extending to the ISCO for Eddington fractions between 10 −3 ≤ L/L Edd ≤ 0.3, holds.While this assumption is likely to not be universally true for all sources, or even for all observations of a single source, it is a reasonable assumption given the observable quantities.
Based on the 36 spins reported in this paper, we can place constraints on the distribution of BH spins observed in XB systems.However, it is important to acknowledge that these values do not take into account any possible observational or selection biases, that could lead to a preferential detection of highly spinning BH systems.Within 1σ uncertainties, 86% of spins are consistent with a ≥ 0.95, 94% are consistent with a ≥ 0.9, and 100% are consistent with a ≥ 0.7, the theoretical upper limit for the spin of a neutron star (NS; Lo & Lin 2011).This suggests that despite possible selection effects, relativistic reflection could be a pragmatic tool for distinguishing between BH and NS XB systems.Given the reported values, 28% of the spins are definitively a ≥ 0.9 within 1σ uncertainty (value minus 1σ error ≥ 0.9).Conversely, looking at the lower limits on the distribution, we find limited support for low BH spins in the BHs in XBs, with ∼ 3% of measurements allowing values a ≤ 0 within 1σ, 5.5% allowing a ≤ 0.4, and ∼ 17% allowing a ≤ 0.6 within 1σ (value minus 1σ error ≤ 0.6).Again, these results strongly suggest that the BHs observed in XBs and those merging in BBH systems observed though GWs indeed come from different distributions, and that their formation paths are determined by the properties of the progenitor stars and by binary interactions both before and after BH formation.
We attempted to illustrate the importance that smallscale source variability between observations can have on inducing changes in the measured BH spin.The true value of the BH spin is unlikely to change on timescales comparable to the age of X-ray astronomy, meaning that the variability in BH spin that we see is entirely owed to our incomplete characterization of the physical processes that create the observed spectra.In the future, this uncertainty will likely be mediated by using more complex theoretical models, which can simultaneously explain the spectral, timing, and X-ray polarization measurements of our target sources.The era of high-resolution spectroscopy started by the launch of XRISM will shed light on narrow spectral features present in X-ray spectra that were previously unaccounted for, and the importance of those for spin measurements is yet to be determined.Analyzing XRISM spectra will likely require more complex theoretical models, which will incorporate a more physically accurate description of the mechanisms at play in XB systems.
Lastly, by systematically analyzing such a large sample of BH XB, we have compiled an extensive data set that illustrates the behavior of the models when applied to a variety of observations taken across an array of physical configurations of the sources.This entire data set will be treated in a future paper, in order to assess the existence of possible trends or degeneracies in the entire parameter space.That study will quantify the importance of correlations between parameters, and their impact on the ability to recover the BH spin.
The best fit for the spectra from both observations was achieved using the diskbb+powerlaw model, with Galactic absorption through TBabs not being statistically required by the data.For ObsID 90402367002, we obtain χ 2 /ν = 67.97/51,with a measured disk temperature of ∼ 0.5 keV and a power-law index Γ ∼ 3.9±0.2.For ObsID 90410345001, we obtain χ 2 /ν = 37.49/34, with a measured disk temperature of ∼ 0.6 keV and an unconstrained power-law index Γ.We conclude that no spin measurement of this source is possible given the available archival NuSTAR observations.

C. FIGURES
This appendix contains all the figures that show the fit residuals and the posterior probability distributions resulting from the MCMC analysis for all the observations of each source.

Figure 1 .
Figure 1.Source and background regions used for extracting spectra of T15 (left) and T37 (right).See text in Sections 3.15 and 3.16.
Figure 35.The middle and bottom panels of Figure 35 show the residuals produced when fitting the spectra with the TBabs*(diskbb+powerlaw) model and the best-performing reflection model, respectively, together with the statistic produced.The posterior probability distributions for BH spin and inclination of the inner accretion disk are shown by the blue curves in the left and right panels of Figure 36, respectively.Using only the posterior distributions from the MCMC analysis of this observation in our combining algorithm produces the black curves shown in Figure 36.We measure a = 0.884 +0.068 −0.109 and θ = 63 +10 −8 [ • ].These results are formally consistent with those of Mori et al. (2019), with larger uncertainties which better encapsulate the systematic uncertainties of the measurement.
, with the left panel showing the distributions for the spin parameter, and the right panel showing the distributions for the inclination parameter.By combining the individual posterior probability distributions using our combining algorithm, we measure a spin of a = 0.987 +0.007 −0.037 and an inclination of θ = 72 +6 −9 [ • ], in good agreement with the results of Feng et al. (2022a) and Coughenour et al. (2023).
who measured a spin of a = 0.967 ± 0.013 and an inclination of θ = 26.4 ± 0.5 [ • ].Due to the proximity of the source to the Sun during the observations, Mode 1 scientific data were unavailable.Therefore, we extracted spectra from the Mode 6 scientific data produced during the observation.The spectra taken during the two observations are shown in the top panels of Figure 43.Fits using the TBabs*(diskbb+powerlaw) model produce the residuals shown in the central panels of Figure 43.The residuals produced by the best-performing reflection models are shown in the bottom panels of Figure 43.

Figure 4 .
Figure 4. Top: change in spin (left) and inclination (right) when compared to the previously existing measurements.Bottom left: same as the top left, but focused on a old ≥ 0.75, for visual clarity.Bottom right: change in spin vs. change in inclination.The markers represent the same sources as in Figure3, with all values taken from Table1.

Figure 5 .
Figure 5. BH spin vs. Eddington fraction (left) and hardness (right), defined as the ratio of the flux in the 10-79 keV band to the flux in the 3-5 keV band, for GRS 1915+105, Cygnus X-1, GX 339-4, and MAXI J1820+070.The vertical dashed lines show an Eddington ratio of 10 −2 in the left panels, and a hardness of Flux 10−79 keV /Flux 10−79 keV = 5.The transparency of the points is proportional to the strength of reflection during the observations.

Figure 6 .
Figure6.Results of simultaneous fits to the inner disk radius (top) and BH spin (bottom) as a function of Eddington fraction for all the observations of GX 339-4 treated in this paper.The transparency of the points is proportional to the strength of reflection, which was used throughout the paper as weighting for combining independent measurements.No evidence for disk truncation is clear, throughout the entire regime of Eddington fraction.

Figure 7 .Figure 8 .
Figure 7.The three panels represent the three observations of AT 2019wey analyzed.The top panels show the unfolded spectra, with the NuSTAR FPMA spectra shown in blue and FPMB spectra in red.The reported best-fit models are shown by the solid lines, while the contributions of the relxill components are shown by the dotted lines and the contributions of the diskbb components are shown through the dashed lines.The middle and bottom panels show the residuals in terms of σ produced when fitting the spectra with TBabs*(diskbb+powerlaw) (middle) and the best-performing reflection models (bottom), respectively, together with the statistic produced.Figure discussed in Section 3.1.

Figure 9 .Figure 10 .
Figure 9. Unfolded spectra and fit residuals for the LMC X-3 data.Explanations are analogous to Figure 7. Figure discussed in Section 3.2.

Figure 18 .
Figure 18.Posterior distributions for the analysis of MAXI J1535-571 data.Explanations are analogous to Figure 8. Figure discussed in Section 3.6.

Figure 26 .
Figure 26.Posterior distributions for the analysis of GX 339-4 data.Explanations are analogous to Figure 8. Figure discussed in Section 3.10.

Figure 27 .
Figure 27.Fit residuals for the IGR J17091-3624 data.Explanations are analogous to Figure 7. Figure discussed in Section 3.11.

Figure 28 .
Figure 28.Posterior distributions for the analysis of IGR J17091-3624 data.Explanations are analogous to Figure 8. Figure discussed in Section 3.11.

Figure 49 .
Figure 49.Fit residuals for the GRS 1915+105 data.Explanations are analogous to Figure 7. Figure discussed in Section 3.22.

Figure 50 .Figure 51 .
Figure 50.Posterior distributions for the analysis of GRS 1915+105 data.Explanations are analogous to Figure 8. Figure discussed in Section 3.22.

Figure 52 .
Figure 52.Posterior distributions for the analysis of Cygnus X-1 data.Explanations are analogous to Figure 8. Figure discussed in Section 3.23.

Figure 53 .Figure 54 .
Figure 53.Unfolded spectra and fit residuals for the V404 Cygnus data.Explanations are analogous to Figure 7.The flaring state of the source is shown in the left panels, and is indicated by " f," while the non-flaring state is shown in the right panels, indicated by " n." Figure discussed in Section 3.24.

Table 1 .
All Current Spin Measurements in XBs Performed Using Our Pipeline