<figure><img src="image_1.png" /></figure>

# The Spin and Orientation of the Black Hole in XTE J1908+094

Paul A. Draghis¹, Jon M. Miller¹, Abderahmen Zoghbi¹, Elias S. Kammoun¹,², Mark T. Reynolds¹, and John A. Tomsick³

¹Department of Astronomy, University of Michigan, 1085 South University Avenue, Ann Arbor, MI 48109, USA; pdraghis@umich.edu

²IRAP, Université de Toulouse, CNRS, UPS, CNES 9, Avenue du Colonel Roche, BP 44346, F-31028, Toulouse Cedex 4, France

³Space Sciences Laboratory, 7 Gauss Way, University of California, Berkeley, CA 94720-7450, USA

Received 2021 May 28; revised 2021 June 29; accepted 2021 July 6; published 2021 October 18

## Abstract

NuSTAR observed the black hole candidate XTE J1908+094 during its 2013 and 2019 outbursts. We use relativistic reflection to measure the spin of the black hole through 19 different assumptions of relxill flavors and parameter combinations. The most favored model in terms of the Deviance Information Criterion (DIC) measures the spin of the black hole to be $a = 0.55_{-0.45}^{+0.29}$, and the inclination to be $\theta = 27_{-3}^{+2}$ degrees (1$\sigma$ statistical errors). We look at the effects of coronal geometry assumptions and the density of the accretion disk on the spin prediction. All 19 tested models provide consistent spin estimates. We discuss the evolution of spin measurement techniques using relativistic reflection in X-ray binaries and discuss the implications of this spin measurement in reconciling the distributions of stellar-mass black hole spin measurements made through X-ray and gravitational wave observations.

Unified Astronomy Thesaurus concepts: Accretion (14); Black hole physics (159); Low-mass x-ray binary stars (939); Stellar mass black holes (1611); X-ray binary stars (1811)

## 1. Introduction

The two main methods for measuring black hole spin in X-ray binaries are disk continuum fitting (see, e.g., McClintock et al. 2006; Morningstar et al. 2014; Steiner et al. 2014) and relativistic reflection (e.g., Miller 2007; Miller & Miller 2015; Reynolds 2020). Spin measurements made using continuum fitting are strongly dependent on prior knowledge about the mass of the black hole, the distance to the system, the inclination of the inner accretion disk, and the mass accretion rate, due to the degeneracy between these parameters in the model. By contrast, the relativistic reflection method requires no prior knowledge about the mass, distance, and mass accretion rate of the black hole in the system, with the inclination of the inner disk being treated as a free parameter. Currently, the main assumption of relativistic reflection models is the nature of the compact corona in the system, with some models assuming a “lamp post” coronal geometry, while others make no prior coronal geometry assumptions. Therefore, relativistic reflection can be applied in systems where, e.g., foreground extinction may inhibit binary mass constraints and is applicable in a larger number of sources.

The two main features of relativistic reflection are the Fe K fluorescence line and a broad flux excess between 20 and 40 keV known as the “Compton hump.” The Fe K spectral line is present at 6.4 keV for neutral gas and at progressively higher energies for ionized gas, up to 6.97 keV for H-like Fe XXVI. The line profile of the Fe K line originating from matter in the proximity of a black hole is “blurred” by relativistic Doppler shifts and gravitational redshifts. This method works under the assumption of an optically thick, geometrically thin accretion disk (Shakura & Sunyaev 1973) that extends all the way to the Innermost Stable Circular Orbit (ISCO). Numerical simulations show that for an Eddington fraction $\le 0.3$, the assumption of an accretion disk extending to the ISCO holds (Shafee et al. 2008) and that any gas at smaller radii is infalling onto the black hole and is optically thin (Reynolds & Fabian 2008). As the size of the ISCO is dependent on the black hole spin in the Kerr metric (Bardeen et al. 1972; Novikov & Thorne 1973), measuring the

blurring of the Fe K line allows an independent measure of the black hole spin.

In the emerging era of gravitational wave (GW) signals from binary black hole (BBH) mergers, the issue of breaking the degeneracy between the mass ratio $q = M_2/M_1$ in the binary and the effective spin parameter $\chi_{\text{eff}} = \frac{M_1\chi_1 + M_2\chi_2}{M_1 + M_2}$ (where $M_{1,2}$ are the masses of the components in the binary and $\chi_{1,2}$ are the components of the spin aligned with the orbital angular momentum) is prohibiting precise measurements of black hole spin (see, e.g., Pürrer et al. 2016; Tiwari et al. 2018). An additional degeneracy between the spins of the two black holes further complicates the spin measurement (Pürrer et al. 2016). As the posterior distributions of GW spin measurements are correlated to the assumed prior distribution (Vitale et al. 2017; Abbott et al. 2019b, 2021), it is important to have educated predictions for the prior distribution of spins (Zevin et al. 2020). The most pragmatic choice for an informative prior distribution comes from a robust measurement of the spins of stellar-mass black holes in X-ray binary systems.

The black hole system XTE J1908+094 was discovered in 2002 February using the Rossi X-ray Timing Explorer (RXTE) Proportional Counter Array (PCA; Woods et al. 2002). Early spectra indicated the presence of broadened Fe line emission, and Miller et al. (2009) measured a spin of $a = 0.75 \pm 0.09$ using BeppoSAX MECS spectra. Since its discovery, the source has undergone three outbursts, in 2003, 2013, and 2019. The 2013 and 2019 outbursts have been observed with NuSTAR (Harrison et al. 2013). During the 2013 outburst, Rushton et al. (2017) have detected resolved, expanding radio jets originating from XTE J1908+094 using the Very Long Baseline Array and the European VLBI Network.

Tao et al. (2015) and Zhang et al. (2015) have indicated the presence of a broadened Fe K line in the 2013 NuSTAR spectra of XTE J1908+094. Still, both papers concluded that spin measurements cannot be made using the data set. The 2019 NuSTAR observation was analyzed in Chatterjee et al. (2021), but no attempt at a spin measurement was made. Motivated by
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Light curves of the two analyzed observations in the 3–79 keV energy band. The light curves are binned so that each point represents 10 seconds. The left panels show the light curve of the entire observation, while the right panels show the light curves split according to their hardness, as they were considered in the extraction of the spectra that were used for the rest of the analysis. The top panels correspond to the 2013 observation, while the bottom panels represent the 2019 observation.</figcaption></figure>

the characteristics of NuSTAR such as its wide passband and its high sensitivity and by previous similar spin measurements using relativistic reflection (see, e.g., King et al. 2014; El-Batal et al. 2016; Draghis et al. 2020; Xu et al. 2020), we revisited the 2013 and 2019 observations of XTE J1908+094. In Section 2 we describe our data reduction process, in Section 3 we describe our analysis and results, and in Section 4 we discuss the implications of the result, connecting it to the future of spin measurements in X-ray binaries and BBH.

## 2. Observations and Data Reduction

NuSTAR observed XTE J1908+094 on 2013 November 8 under ObsID 80001014002 for a net exposure of 45 ks and on 2019 April 10 under ObsID 90501317002 for 41 ks. During both observations, the source was in the soft state (Negoro et al. 2013; Ludlam et al. 2019). The data were reduced using the routines in HEASOFT v6.27.1 through the NuSTARDAS pipeline v1.9.2 and CALDB v20200510. The source spectra were extracted from circular regions centered on the source position with radii of 100" in the two FPM NuSTAR sensors, and regions of the same size were used for extraction of background rates. The spectra were grouped using the “ftgrouppha” tool, through the optimal binning scheme described by Kaastra & Bleeker (2016). One potential issue with this binning scheme would be a small number of counts in some bins, particularly toward high energies, for which $\chi^2$ statistics would not be appropriate. Therefore, we compared the fits of our models on the spectra binned with the optimal scheme with the same spectra when grouped by oversampling the resolution by a factor of 3 and requiring a minimal signal-to-noise ratio of 6. The best-fit parameters were strongly consistent between binning schemes. The left panels in Figure 1 show the light curves of the two NuSTAR observations. For the 2013 observation (top left panel), the central part of the light curve shows a 20% increase in count rate over the 28 counts s⁻¹ continuum, while the 2019 observation (bottom left panel) shows variability on the order of 40%.

For both observations, we extracted the light curves in the 3–8 keV and 8–20 keV bands and compared the ratio of the count rates in the two energy bands to the total count rate in the 3–20 keV band. The lower limit of 3 keV was chosen as the

<figure><img src="image_3.png" /><figcaption>**Figure 2.** The hardness ratio throughout the two observations computed as the ratio of the 8–20 keV count rate to the 3–8 keV count rate, as a function of the total 3–20 keV count rate. The blue points represent the 2013 observation, while the yellow points represent the 2019 one. The horizontal dashed lines show the hardness cuts placed on the two observations to divide them into the “high hardness” and “low hardness” spectra. The light curves were binned so that each point represents 200 s.</figcaption></figure>

lower bound of the NuSTAR passband. The upper limit of 20 keV was chosen to be just under the Compton hump. The 8 keV break point was chosen to separate the effects of the Fe complex and the disk emission from coronal emission typical of X-ray binary systems. This hardness ratio is shown in Figure 2. For both observations, the instantaneous hardness appears to have a bimodal distribution. Therefore, we split the two observations into regions with hardness greater and lower than 0.16 for the 2019 observation and greater/lower than 0.1 for the 2013 observation. We generated GTI files using the XSELECT tool by filtering the intensity as described above, and re-extracted the spectra of the low and high hardness intervals of each observation. The right panels in Figure 1 show the intervals of the light curve corresponding to the high and low hardness part of the observations. For the 2013 observation, the high hardness interval corresponds to the central increase in overall flux, and results in the same division of the spectrum as that used by Tao et al. (2015) and Zhang et al. (2015). Similarly, for the 2019 observation the increase in hardness is associated with an increase in the total count rate of the source. The produced spectra are shown in panel (a) of Figure 3. We used the entire NuSTAR energy band (3–79 keV) throughout the spectral fits.

Upon the completion of the analysis, we tested the robustness of the split of the observations into hardness regimes by re-extracting the spectra after redefining the hardness ratio as the count rates in the 3–6 keV and 6–20 keV bands. We note that the final results are strongly consistent with the initial measurements, implying that in this case, the choice of definition of hardness ratio has little to no impact on the spin measurement.

## 3. Analysis and Results

The spectral fitting was carried out using XSPEC v12.11.0 m (Arnaud 1996), with the quality of the fit being evaluated using $\chi^2$ statistics. The CALDB version used accounts for the low energy excess of the spectra produced by the NuSTAR FPMA sensor due to the tear in its thermal blanket, described by Madsen et al. (2020). While a multiplicative constant is often allowed to vary between the spectra from the two FPM NuSTAR sensors, in this particular case the data did not require
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Panel (a) shows the spectrum of XTE J1908+094. The red, green, blue, and purple points represent the 2013 low hardness, 2013 high hardness, 2019 low hardness, and 2019 high hardness data groups, respectively. Panel (b) shows the residuals in terms of $\sigma$ for the $\text{TBabs}^*(\text{diskbb+powerlaw})$ model. Panel (c) shows the residuals of the $\text{TBabs}^*(\text{diskbb+relxilllp})$ model, which is the best performing model of the 19 tested in terms of the Deviance Information Criterion (model 6).</figcaption></figure>

it. By splitting the two observations into two hardness regimes, we obtain four data groups, each consisting of the FPMA/B spectra, fit jointly.

Visually inspecting the spectra in panel (a) of Figure 3 suggests that, on the zeroth order, models need at least two components to accommodate the apparent slope change at ~9 keV. We fit the spectra with the combination of a disk blackbody component `diskbb` (Mitsuda et al. 1984), characteristic of accretion disks around compact objects, and a power law. To account for interstellar absorption, we included the multiplicative component `TBabs` (Wilms et al. 2006), with the abundances computed according to Wilms et al. (2000) and photoionization cross sections computed by Verner et al. (1996). The residuals of this model are shown in panel (b) of Figure 3, with the fit producing $\chi^2/\nu = 2740.25/1621 = 1.69$. The main unmodeled feature is a strong emission line consistent with the Fe K line complex. The decrease in flux both below and above the Fe K line suggests that the fitting procedure altered the continuum shape of the disk blackbody in order to balance the negative and positive residuals, therefore impacting the shape of the continuum at low energy and inducing an apparent flux excess in the 3–4 keV range. Additionally, a broad excess is seen in the 20–40 keV range. These features are consistent with relativistic disk reflection.

For the measured disk temperature $kT_{\text{in}} \sim 0.7$ keV, the energy distribution peaks around ~2 keV, just below the lower

limit of the NuSTAR passband. As the high energy tail of the disk continuum easily extends into the NuSTAR band, the ability to constrain interstellar absorption and the temperature of the accretion disk would not be improved significantly through the addition of simultaneous observations using the Swift X-ray telescope (XRT; Burrows et al. 2007). Zhang et al. (2015) pointed out that when fitting the Swift and NuSTAR spectra from the same day simultaneously, the residuals show a difference between the spectra from the two instruments in the overlapping energy band. As Swift detectors suffer from pile-up, the reduced sensitivity obtained when extracting annuli does not permit an improvement of the constraints placed on parameters.

### 3.1. The relxill Family of Models

The relativistic reflection features were modeled by replacing the `powerlaw` component of the model with multiple models that are part of the `relxill` family of models (Dauser et al. 2014; García et al. 2014). We used `relxill` v1.4.0. Panel (c) of Figure 3 shows the residuals of one of such model, with the reflection component being modeled by `relxilllp` (described below), thus making the full model `TBabs*(diskbb+relxilllp)`. The quality of the fit is drastically improved, producing $\chi^2/\nu = 1707.44/1604 = 1.06$.

The `relxill` family of models is now a standard tool in modeling relativistic reflection. Variants of `relxill` allow for the probing of the effects on the spectrum of different physical assumptions, such as the geometry of the corona, density and ionization of the accretion disk, iron abundance, geometry of the accretion disk, and spin of the black hole. As the goal of this analysis is to measure the spin of the black hole, $a = cJ/GM^2$ ($-1 \le a \le 1$), we aimed to explore the entire range of parameter combinations and theoretical assumptions that the `relxill` family of models provides.

All `relxill` models include the spin of the black hole as a free parameter. The geometry of the accretion disk is constrained through the inner and outer radius of the disk ($r_{\text{in}}$ and $r_{\text{out}}$) and the inclination of the inner accretion region of the accretion disk. For our analysis, we fixed the inner radius of the accretion disk at the ISCO of the black hole, as this assumption is in agreement with MHD simulations such as those by Reynolds & Fabian (2008). The outer radius of the disk was fixed at a large value of $r_{\text{out}} = 990 r_g$. The inclination of the inner accretion disk was allowed to vary as a free parameter. Note that the inclination of the inner disk does not necessarily match that of the outer regions of the accretion disk and of the binary system, as a result of the Bardeen–Petterson effect (see, e.g., Bardeen & Petterson 1975; Nealon et al. 2015; Liska et al. 2020). Still, the orientation of the black hole spin and the inner accretion disk can precess due to the Lense–Thirring effect and align itself with the outer accretion disk (see, e.g., Banerjee et al. 2019).

The iron abundance of the reflector measured relative to the solar abundance $A_{\text{Fe}}$ and the ionization of the reflector ($\xi = L/nr^2$, with L being the luminosity of the emitting source, n the hydrogen number density of the reflector, and r the distance between source and reflector) are allowed to vary in all models. The illuminating flux is modeled as a power law with index $\Gamma$ and a high energy cutoff $E_{\text{cut}}$. The power-law index is a free parameter in all models, while the high energy cutoff is free in models having a fixed density of the accretion disk of $n = 10^{15} \text{ cm}^{-3}$ (`relxilllp` and `relxill`), and fixed by
construction of the models to $E_{\text{cut}} = 300$ keV in models that allow higher densities of the accretion disk (relxillD).

An alternative description of the emissivity of the corona is to model it as a thermal Comptonization continuum. The relxill family of models does offer this possibility through the relxillCp and relxilllpCp models. In order for black hole coronae to produce the observed X-ray fluxes while being powered by hot electrons, the size of coronae needs to be on the order of ~1000 $r_g$ (Merloni & Fabian 2001). A corona of this physical size cannot produce short-term variability on the timescales observed in some sources (see, e.g., Lee et al. 2000; Fabian et al. 2015). Merloni & Fabian (2001) suggested that the energy must be stored in the corona as magnetic fields generated by a sheared rotator. Another possibility is a combination of the two: the compact base of a jet surrounded by a larger cloud of hot electrons. As a test, we compare our best model (see below, Section 3.3) with the equivalent Comptonization model. The fit becomes worse by $\Delta\chi^2 \approx 10$ for no change in the number of degrees of freedom, while the best-fit parameters remain strongly consistent with the initial results. As the shape of the emissivity of the corona does not appear to influence the spin measurement, we continued to only analyze models that use a power-law emissivity.

In the relxilllp model, the geometry of the corona is modeled as a compact emitter at height $h$ above the accretion disk (the “lamp post” assumption). In relxill and relxillD, the geometry of the corona is unconstrained, but its emissivity is modeled as a broken power law with radius ($J \propto r^{-q_1}$ for $r < r_{\text{break}}$, and $J \propto r^{-q_2}$ for $r > r_{\text{break}}$). Additionally, all relxill models have a parameter representing the normalization of the flux of the component (norm_r) and a “reflection fraction” parameter (refl Frac) that represents the ratio of the intensity emitted toward the accretion disk and the intensity escaping to infinity in the frame of the primary source (Dauser et al. 2016a). For models assuming the lamp post geometry, this is calculated using ray tracing simulations as described by Dauser et al. (2014). In this version of the model, for the relxill variants that do not assume a coronal geometry, the reflection fraction is approximated by taking the ratio of the reflected flux to the direct flux in the 20–40 keV band (Dauser et al. 2016b).

The relxillpion variant of relxill replaces the constant ionization parameter with a ionization gradient modeled as a power law with radius. With this model, the solution converges to the same parameter combination as in the case with constant ionization, with the power-law index for the ionization gradient being consistent with zero. While this does not exclude the possibility of an ionization gradient across the accretion disk, this data set is not able to constrain any possible variation in ionization. Therefore, we chose not to explore models using variable ionization further. The lack of an apparent ionization gradient is consistent with a steep emissivity profile. A flatter emissivity would include more emission from large radii that could, in principle, have lower ionization. This would make the effects of an ionization gradient more apparent.

## 3.2. Markov Chain Monte Carlo Analysis

We explored a total of 19 variations of the above models. For each individual combination of model and parameter assumptions, we ran the default XSPEC fitting algorithm to find the “best-fit” parameter combination that minimizes the $\chi^2$ statistic given the

data. This set of parameters was used to construct a proposal distribution for running a Monte Carlo Markov Chain (MCMC). The proposal distribution was generated from a Gaussian distribution around the “best-fit” parameters. The chains were run with 200 walkers for a total of $2 \times 10^6$ steps, using the XSPEC EMCEE implementation written by A. Zoghbi$^4$. The number of walkers was chosen as an integer on the order of a few times the number of free parameters. Sokal (1996) suggested running MCMC chains with more steps than ~1000 × $\tau_f$, where $\tau_f$ represents the integrated autocorrelation time. We computed $\tau_f$ using the prescription given by Goodman & Weare (2010) and determined that running the chains for $2 \times 10^6$ steps is equivalent to 3–4000 × $\tau_f$, ensuring convergence of the chains. The first $5 \times 10^5$ steps of the chains were then assumed to be a “burn in” phase and therefore disregarded.

Using the MCMC posterior distribution, we computed the Deviance Information Criterion (DIC; Spiegelhalter et al. 2002) for each model. We used the DIC to quantify the goodness of the fit and to distinguish between the tested models. We used this approach instead of an F-test in order to properly account for how well the model fits the data and how complex the model is, all while not only looking at one single parameter combination of maximum likelihood, but at the entire distribution around the best-fit value. Table 3 in Appendix A presents the entire set of models tested, together with the relxill variant used, the assumptions made, the best-fit $\chi^2$ and number of degrees of freedom, reduced $\chi^2$, computed DIC, and model ranking based on the DIC.

Figure 4 shows the histograms of the posterior distributions of the MCMC chains for the black hole spin for the 19 tested models. The top panel shows the model variants that use relxilllp to model the reflection component, the middle panel shows the distributions for the variants using relxill, and the lower panel shows the models using relxillD. The models are presented below. The vertical lines represent the median of the posterior distribution for the spin parameter in each model. The highlighted green curve in the top panel and the blue curve in the bottom panel correspond to models 6 and 18, presented in detail below. The spin measurements are consistent among all models, suggesting a robustness of the measurement against the assumed relxill flavor used.

## 3.3. Best “Lamp Post” Model

Based on the computed DIC, the model that best fits the data in the lamp post assumption is that where the height of the corona is linked between the low and high hardness spectra from the same observation (model 6 in 3). Table 1 shows the median values across the posterior sample in the MCMC chains with $1\sigma$ uncertainties. For the black hole spin, the number in parentheses represent the value at which $\chi^2$ is minimized. The residuals of the fit using this model are shown in panel (c) of Figure 3. The number in parentheses in the $\chi^2/\nu$ row represents the median $\chi^2$ value across the posterior distribution from the MCMC chains, with its $1\sigma$ variation. The measured spin of the black hole for the best performing model in terms of DIC is $a = 0.55_{-0.45}^{+0.29}$. The posterior distribution of the spin parameter in the MCMC chain is shown in the top panel in Figure 4 by the highlighted green line. This measurement is consistent with all other measurements made using the other tested models.

<footnote>$^4$ https://zoghbi-a.github.io/xspec_emcee/</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Histogram of the posterior distributions of the black hole spin parameter in all 19 tested models. The top panel shows the models using *relxilllp*. The middle panel shows the models using *relxill*, and the bottom panel shows models using *relxillD*. The vertical lines represent the median value of the spin posterior distribution for each model. The highlighted green line in the top panel represents the best performing “lamp post” model in terms of DIC (model 6). The blue curve in the bottom panel shows the model performing best in terms of DIC from the *relxill* variants that adopt no assumptions about the geometry of the corona.</figcaption></figure>

Model 6 links the height of the compact corona in the *relxilllp* within the 2013 and 2019 observations. The measured values for the height of the corona are $h_{2013} = 14_{-2}^{+2} r_g$ for the 2013 observation and $h_{2019} = 6_{-1}^{+1} r_g$ for the 2019 observation. If instead we link the height of the corona between all four data groups, the measured value is the average of the previous two measurements, $h = 10_{-1}^{+1} r_g$. This corresponds to model 7 in Table 3, which produces the second-best DIC of the 19 models. Allowing the height of the corona to be free in all data groups (model 5) produces a similar $\chi^2$ value at the cost of two extra free parameters, resulting in a worse DIC, ranking sixth among the 19 tested models. The measurements for the height of the corona in this model are $h_{2013 low} = 14_{-2}^{+3} r_g$, $h_{2013 high} = 15_{-2}^{+3} r_g$, $h_{2019 low} = 6_{-1}^{+1} r_g$, and $h_{2019 high} = 6_{-1}^{+1} r_g$, consistent with the values obtained when the height of the corona is linked within each observation. This result suggests that in the lamp post geometry, the corona of XTE J1908+094 moves from ~14 $r_g$ in 2013 to ~6 $r_g$ in 2019.

The measured inclination of the inner region of the accretion disk through this model is $Incl = 27_{-3}^{+2}$ degrees. The measured Fe abundance is $A_{Fe} = 6_{-1}^{+2}$ of the solar value. Tomsick et al. (2018)

found that for Cygnus X-1, a source with a previously reported high Fe abundance ($A_{Fe} \sim 10$), fitting the data with a high density variant of the *reflionx* model (Ross & Fabian 2005) produces better fits than using *relxill* variants, even with a solar Fe abundance. For Cygnus X-1, the value of the disk density measured by Tomsick et al. (2018) was $n \sim 4 \times 10^{20} \text{ cm}^{-3}$, much higher than the range of densities that the *relxill* family of models can probe. We were unable to obtain a fit of comparable quality to that of the other tested models when using the model presented in Tomsick et al. (2018). If instead we fix $A_{Fe} = 1$ in our tested models, the fits become significantly worse, with an increase of $\chi^2 \ge 65$ for one fewer degree of freedom across all models. In all models, including the high density variants of *relxill*, the preferred Fe abundance is high.

As both relativistic Doppler shifts and gravitational redshifts have the effect of producing an asymmetry of spectral lines, the black hole spin parameter and the inclination of the inner accretion disk can be partially degenerate. The middle left panel in Figure 5 shows the $1\sigma$, $2\sigma$, and $3\sigma$ confidence intervals (in red, blue, and green) of the black hole spin parameter in the posterior MCMC distribution in relation to the inclination of the accretion disk in model 6. The top left and middle right panels show the 1D histograms of the posterior samples in the MCMC run for the spin and inclination of the inner disk. The solid red lines represent the median of the distributions, while the dashed red lines represent the $\pm 1\sigma$ confidence limits of the median. We note that no apparent degeneracy is present between the mentioned parameters in this measurement. The bottom left panel in Figure 5 shows the same confidence intervals for the spin measurement in relation to the $\chi^2$ of the specific parameter combination containing the given spin values. The bottom right panel shows the 1D histogram of the $\chi^2$ values in the MCMC chain. Note that the quality of the fit is not driven mainly by the black hole spin, indicating that error searches across a single parameter can be unreliable. This confirms the need to explore the effects of varying the entire parameter set simultaneously in an error analysis, in order to probe the full set of physical phenomena covered by the model. Apart from an apparent degeneracy between the temperature of the *diskbb* components and the normalization of the component, no other degeneracies are noticeable. For reference, Figure 7 in Appendix B shows the corner plot for a few parameters of interest in the analysis.

The emission from the accretion disk was modeled using *diskbb* (Mitsuda et al. 1984). The measured disk temperature in this model is $kT_{\text{in 2013 low}} = 0.709_{-0.003}^{+0.003}$ keV for the low hardness part of the 2013 spectrum, $kT_{\text{in 2013 high}} = 0.714_{-0.004}^{+0.004}$ keV for the high hardness part of the 2013 spectrum, $kT_{\text{in 2019 low}} = 0.598_{-0.004}^{+0.004}$ keV for the low hardness part of the 2019 spectrum, and $kT_{\text{in 2019 high}} = 0.606_{-0.005}^{+0.005}$ keV for the high hardness part of the 2019 spectrum. This is consistent with a disk of temperature $kT_{\text{in 2013}} \sim 0.71$ keV for the 2013 observation and $kT_{\text{in 2019}} \sim 0.60$ keV for the 2019 observation, indicating a constant disk temperature throughout the duration of the observation, but variable over longer periods of time.

The normalization of the *diskbb* component is defined as $(R_{\text{in}}/D_{10})^2 \cos\theta$ where $R_{\text{in}}$ is the apparent inner disk radius in km, $D_{10}$ is the distance to the source in units of 10 kpc, and $\theta$ is the inclination angle of the accretion disk. Assuming that the accretion disk extends to the ISCO for both observations, that the distance to the source is constant, and that the inclination does not vary in time, one would expect the normalization of
<table><caption>**Table 1**

Results from Model 6</caption><thead><tr><th>Parameter</th><th>2013 Low</th><th>2013 High</th><th>2019 Low</th><th>2019 High</th></tr></thead><tbody><tr><td>N<sub>H</sub> (×10<sup>22</sup> cm<sup>-2</sup>)</td><td></td><td>3.6<sup>+0.1</sup><sub>-0.1</sub></td><td></td><td></td></tr><tr><td>a</td><td></td><td>0.55<sup>+0.29</sup><sub>-0.45</sub> (0.71)</td><td></td><td></td></tr><tr><td>Incl (°)</td><td></td><td>27<sup>+2</sup><sub>-3</sub></td><td></td><td></td></tr><tr><td>A<sub>Fe</sub></td><td></td><td>6<sup>+2</sup><sub>-1</sub></td><td></td><td></td></tr><tr><td>h (r<sub>g</sub>)</td><td></td><td>14<sup>+2</sup><sub>-2</sub></td><td></td><td>6<sup>+1</sup><sub>-1</sub></td></tr><tr><td>kT<sub>in</sub> (keV)</td><td>0.709<sup>+0.003</sup><sub>-0.003</sub></td><td>0.714<sup>+0.003</sup><sub>-0.004</sub></td><td>0.598<sup>+0.004</sup><sub>-0.004</sub></td><td>0.606<sup>+0.005</sup><sub>-0.005</sub></td></tr><tr><td>norm<sub>d</sub> (×10<sup>3</sup>)</td><td>1.25<sup>+0.04</sup><sub>-0.04</sub></td><td>1.17<sup>+0.04</sup><sub>-0.04</sub></td><td>1.07<sup>+0.06</sup><sub>-0.05</sub></td><td>0.98<sup>+0.06</sup><sub>-0.05</sub></td></tr><tr><td>Γ</td><td>2.08<sup>+0.02</sup><sub>-0.03</sub></td><td>2.21<sup>+0.02</sup><sub>-0.02</sub></td><td>2.12<sup>+0.02</sup><sub>-0.03</sub></td><td>2.16<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>Log(ξ)</td><td>3.8<sup>+0.1</sup><sub>-0.1</sub></td><td>4.3<sup>+0.1</sup><sub>-0.1</sub></td><td>3.9<sup>+0.1</sup><sub>-0.1</sub></td><td>4.0<sup>+0.1</sup><sub>-0.1</sub></td></tr><tr><td>E<sub>cut</sub> (keV)</td><td>500<sup>+200</sup><sub>-200</sub></td><td>700<sup>+200</sup><sub>-200</sub></td><td>600<sup>+200</sup><sub>-200</sub></td><td>700<sup>+200</sup><sub>-200</sub></td></tr><tr><td>refl<sub>frac</sub></td><td>0.6<sup>+0.1</sup><sub>-0.1</sub></td><td>0.7<sup>+0.1</sup><sub>-0.1</sub></td><td>0.9<sup>+0.2</sup><sub>-0.2</sub></td><td>0.8<sup>+0.2</sup><sub>-0.2</sub></td></tr><tr><td>norm<sub>r</sub> (×10<sup>-4</sup>)</td><td>13<sup>+1</sup><sub>-1</sub></td><td>21<sup>+2</sup><sub>-2</sub></td><td>11<sup>+2</sup><sub>-2</sub></td><td>17<sup>+3</sup><sub>-3</sub></td></tr><tr><td>χ<sup>2</sup>/bins</td><td>531.42/432</td><td>398.82/417</td><td>376.17/388</td><td>401.03/401</td></tr><tr><td>χ<sup>2</sup>/ν</td><td colspan="4">1707.44/1604 = 1.06 (1735<sup>+8</sup><sub>-7</sub>)</td></tr></tbody></table>

**Note.** In the row presenting the black hole spin, the number in parentheses represents the value that produced the lowest $\chi^2$ in the analysis. In the $\chi^2/\nu$ row, the number in parentheses represents the median $\chi^2$ value in the MCMC analysis with associated $\pm 1\sigma$ variation.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Histogram of the spin–inclination parameter space for model 6 (middle left panel) and of the spin–$\chi^2$ space (bottom left panel). The $1\sigma$, $2\sigma$, and $3\sigma$ confidence intervals are represented by the red, blue, and green contours in the two aforementioned panels. The top left, middle right, and bottom right panels show the 1D histograms of the posterior distribution in the MCMC chain of the spin, inclination, and $\chi^2$. The solid red lines represent the median of the distributions and the dashed red lines represent the $\pm 1\sigma$ uncertainty around the median.</figcaption></figure>

the *diskbb* component to be constant. Models 1, 2, and 3 fix the normalization of the *diskbb* component across all data groups under the same assumption about the corona height as models 5, 6, and 7, while maintaining the same reflection component (*relxillp*). Additionally, model 4 fixes the normalization and temperature of *diskbb* and the height of the corona within one observation. Models 1–4 produce overall worse fits than models 5–7 where the normalization of the accretion disk component is allowed to vary, indicating that the

*diskbb* model does not fully capture the range of physical effects occurring in the accretion disk. If we characterize the emission from the accretion disk using *kerrbb* (Li et al. 2005), which models the radiation from a multi-temperature blackbody accretion disk around a Kerr black hole, the fits are similar in quality as when using the *diskbb* component. Still, the addition of multiple free parameters drives the DIC to much higher values, and the degeneracy between the black hole mass, distance to the source, mass accretion rate, and disk inclination cannot be broken without prior knowledge of these parameters. Therefore, we note that while *diskbb* does not fully describe all physical effects on the emission from accretion disk, its simplicity makes it the preferred method over more complicated models. In practice, the emergent X-ray spectrum from an accretion disk depends on the spectral hardening factor (Shimura & Takahara 1995), which must depend on the atmospheric properties of the disk, such as gas density, temperature, ionization state, and magnetization (Salvesen & Miller 2021). Therefore, it is possible for the disk radius to remain constant, while the inferred radius and disk temperature be observed to vary.

In all *relxill* models, the illuminating flux is modeled as a power law with index $\Gamma$ and high energy cutoff $E_{\text{cut}}$. This model is unable to constrain the high energy cutoff, returning values between $E_{\text{cut}} \sim 300-900$ keV. The estimated power-law index is $\Gamma_{2013 \text{ low}} = 2.08_{-0.03}^{+0.02}$ for the low hardness part of the 2013 spectrum, $\Gamma_{2013 \text{ high}} = 2.21_{-0.02}^{+0.02}$ for the high hardness part of the 2013 spectrum, $\Gamma_{2019 \text{ low}} = 2.12_{-0.03}^{+0.02}$ for the low hardness part of the 2019 spectrum, and $\Gamma_{2019 \text{ high}} = 2.16_{-0.02}^{+0.02}$ for the high hardness part of the 2019 spectrum. We note that the power-law index $\Gamma$ increases for both observations between the low and high hardness regimes. While this may be contrary to expectations, it is not uncommon, with similar behavior shown by Sobczak et al. (2000) in their monitoring of XTE J1550 564.

The measured ionization of the accretion disk (Log $\xi$) takes high values $\ge 3.7$. Similarly to the case of the power-law index, the ionization parameter increases during both observations
<table><caption>**Table 2**

Results from Model 18</caption><thead><tr><th>Parameter</th><th>2013 Low</th><th>2013 High</th><th>2019 Low</th><th>2019 High</th></tr></thead><tbody><tr><td>N<sub>H</sub> (×10<sup>22</sup> cm<sup>-2</sup>)</td><td></td><td>3.7<sup>+0.1</sup><sub>-0.1</sub></td><td></td><td></td></tr><tr><td>a</td><td></td><td></td><td>0.5<sup>+0.3</sup><sub>-0.6</sub> (0.74)</td><td></td></tr><tr><td>Incl (<sup>1</sup>&deg;))</td><td></td><td></td><td>33<sup>+7</sup><sub>-3</sub></td><td></td></tr><tr><td>A<sub>Fe</sub></td><td></td><td></td><td>6<sup>+1</sup><sub>-1</sub></td><td></td></tr><tr><td>Log(n) (cm<sup>-3</sup>)</td><td></td><td></td><td>19<sup>*</sup></td><td></td></tr><tr><td>R<sub>br</sub> (r<sub>g</sub>)</td><td></td><td></td><td>12<sup>*</sup></td><td></td></tr><tr><td>q<sub>1</sub></td><td>6<sup>+3</sup><sub>-3</sub></td><td>6<sup>+2</sup><sub>-2</sub></td><td>5<sup>+2</sup><sub>-2</sub></td><td>6<sup>+2</sup><sub>-2</sub></td></tr><tr><td>q<sub>2</sub></td><td>2<sup>+2</sup><sub>-1</sub></td><td>3<sup>+1</sup><sub>-1</sub></td><td>3<sup>+3</sup><sub>-2</sub></td><td>3<sup>+3</sup><sub>-2</sub></td></tr><tr><td>kT<sub>in</sub> (keV)</td><td>0.704<sup>+0.004</sup><sub>-0.004</sub></td><td>0.710<sup>+0.004</sup><sub>-0.004</sub></td><td>0.597<sup>+0.005</sup><sub>-0.005</sub></td><td>0.603<sup>+0.005</sup><sub>-0.005</sub></td></tr><tr><td>norm<sub>d</sub> (×10<sup>3</sup>)</td><td>1.31<sup>+0.05</sup><sub>-0.06</sub></td><td>1.23<sup>+0.05</sup><sub>-0.05</sub></td><td>1.10<sup>+0.07</sup><sub>-0.06</sub></td><td>1.02<sup>+0.07</sup><sub>-0.06</sub></td></tr><tr><td>Γ</td><td>2.03<sup>+0.02</sup><sub>-0.02</sub></td><td>2.15<sup>+0.02</sup><sub>-0.02</sub></td><td>2.04<sup>+0.02</sup><sub>-0.02</sub></td><td>2.09<sup>+0.02</sup><sub>-0.02</sub></td></tr><tr><td>Log(ξ)</td><td>3.60<sup>+0.09</sup><sub>-0.09</sub></td><td>3.94<sup>+0.09</sup><sub>-0.10</sub></td><td>3.70<sup>+0.10</sup><sub>-0.10</sub></td><td>3.71<sup>+0.09</sup><sub>-0.10</sub></td></tr><tr><td>refl<sub>frac</sub></td><td>0.25<sup>+0.05</sup><sub>-0.04</sub></td><td>0.31<sup>+0.07</sup><sub>-0.06</sub></td><td>0.31<sup>+0.07</sup><sub>-0.05</sub></td><td>0.29<sup>+0.08</sup><sub>-0.05</sub></td></tr><tr><td>norm<sub>r</sub> (×10<sup>-4</sup>)</td><td>9.0<sup>+0.7</sup><sub>-0.6</sub></td><td>13<sup>+2</sup><sub>-1</sub></td><td>4.5<sup>+0.4</sup><sub>-0.4</sub></td><td>6.8<sup>+0.7</sup><sub>-0.7</sub></td></tr><tr><td>χ<sup>2</sup>/bins</td><td>531.58/432</td><td>400.80/417</td><td>374.67/388</td><td>398.39/401</td></tr><tr><td>χ<sup>2</sup>/ν</td><td colspan="4">1705.44/1602 = 1.06 (1741<sup>+9</sup><sub>-8</sub>)</td></tr></tbody></table>

Note.\* represents a parameter fixed in the fit. In the row presenting the black hole spin, the number in parentheses represents the value that produced the lowest χ² in the analysis. In the χ²/ν row, the number in parentheses represents the median χ² value in the MCMC analysis with associated ± 1σ variation.

between the low and high hardness parts of the spectrum. Lastly, the measured reflection fraction in relxill is broadly consistent with a value ~0.7 and does not significantly vary between observations.

### 3.4. Best Geometry-free Model

While relxillp assumes a lamp post coronal geometry, the geometry of the corona is unspecified in relxill and relxillD. Therefore, models 8–19 in Table 3 test the impact of the coronal geometry assumption on the spin measurement. Model 8 replaces the relxillp component with relxill. The parameters describing the emissivity of the corona ($q_1$, $q_2$, and $R_{br}$) are allowed to vary freely. We note that for this model, the breaking radius $R_{br}$ is unconstrained, so model 9 fixes this to a value of $R_{br} = 10 r_g$. Model 10 allows the breaking radius to vary, but links the normalization of the diskbb component due to the arguments mentioned above. Component relxill allows the high energy cutoff $E_{cut}$ of the power law describing the incident flux to vary, but works under the assumption of an accretion disk density of $Log(n) = 15$. To probe the effects of the disk density on the fits, we replaced the relxill component in the model with relxillD. Models 11–13 adopt the same assumptions about coronal emissivity and diskbb normalization as models 8–10, respectively, but with the disk density fixed to $Log(n) = 15$. Note that this differs from models 8–10 through the fact that relxillD fixed the high energy cutoff to $E_{cut} = 300$ keV, while relxill allows it to vary. Models 14–16 and 17–19 are the analogs of models 11–13 for $Log(n) = 17$ and $Log(n) = 19$, respectively.

Of the models that make no assumption about the geometry of the corona, model 18 performs best in terms of DIC and ranks third among all 19 models. Model 18 uses the relxillD component, fixes the density of the accretion disk to $Log(n) = 19$, and $R_{br} = 12 r_g$. We note that, by fixing $R_{br} = 12 r_g$, models 18, 12, 15, and 9 perform better overall than the other variants of relxill/D ranking 3rd, 4th, 5th, and 7th of the 19 models in terms of DIC. All

of these models produce black hole spin measurements that are consistent within the 1σ uncertainties, as shown by Figure 4. In particular, model 18 predicts a black hole spin of $a = 0.5_{-0.6}^{+0.3}$, in good agreement with the result of the best performing lamp post model. The entire set of parameters of model 18 together with their 1σ uncertainties are shown in Table 2, and the histogram of the posterior distribution of the spin parameter in the MCMC chain is shown by the highlighted blue curve in the bottom panel of Figure 4.

This model produces similar constraints on the parameters as the best performing lamp post model. The measured inclination of the accretion disk for model 18 is $Incl = 33_{-3}^{+7}$ degrees, slightly larger than the one measured by model 6. The measured Fe abundance is consistent with the previous result, $A_{Fe} = 6_{-1}^{+1}$, despite the higher disk density. The temperature and normalization of the diskbb component are once again similar to the ones measured by model 6. The measurements of the power-law index Γ and ionization parameter of the accretion disk are once again slightly lower, but show the same trend of increasing in the hard part of the spectra when compared to the soft part, in both observations. The reflection fraction parameter in the model also decreases to a value of ~0.3, once again consistent between the two observations.

It was argued by Wilkins & Fabian (2012) that the emissivity profiles in active galactic nuclei (AGNs) have a steep power-law index with radius between 6 and 8 for the inner regions of the accretion disk, which flattens for intermediate distances to almost a constant, and then tends to be a constant power-law index of 3 over the outer regions of the disk. While current theoretical models do not yet have the capability of describing the emissivity of the illuminating source as a three component power law, the results of this model on the emissivity of the corona in an X-ray binary follow the same general trend: a steep value $q_1 = 6 \pm 3$ for $r < 12 r_g$ and a lower value $q_2 = 3 \pm 2$ for $r > 12 r_g$. Models that link $q_1$ to $q_2$ for this data set perform poorly in terms of χ² and have therefore not been explored further.
## 4. Discussion

We revisited the NuSTAR observations of XTE J1908+094 during its 2013 and 2019 outburst. We divided each observation into two regimes, corresponding to an increase in hardness. Fitting the four data groups together using a set of 19 models that explore the entire set of physical assumptions that the *relxill* family of models covers returned consistent values of the black hole spin. The best performing model in terms of DIC (model 6) measures the black hole spin to have an intermediate value, $a = 0.55_{-0.45}^{+0.29}$. This model assumes a lamp post coronal geometry and fixes the height of the corona within the 2013 and 2019 observations. The corona appears to shift from $h_{2013} = 14 \pm 2 r_g$ in 2013 to $h = 6 \pm 1 r_g$ in 2019. The inclination of the inner accretion disk is low, $\text{Incl} = 27_{-3}^{+2}$ degrees, the measured Fe abundance is $A_{\text{Fe}} = 6_{-1}^{+2}$ times the solar value, and the ionization of the accretion disk is high (Log $\xi \sim 4$) during both observations, and increases with an increase in hardness.

For the 2013 observation of XTE J1908+094, the measured unabsorbed flux in the low and high hardness states was $F_{2013,low} = 1.20(6.44) \times 10^{-9} \text{erg cm}^{-2} \text{s}^{-1}$ and $F_{2013,high} = 1.36(6.77) \times 10^{-9} \text{erg cm}^{-2} \text{s}^{-1}$ in the 3–79 (0.5–100) keV band. For the 2019 observation, the measured fluxes are $F_{2019,low} = 0.42(2.78) \times 10^{-9} \text{erg cm}^{-2} \text{s}^{-1}$ and $F_{2019,high} = 0.51(2.91) \times 10^{-9} \text{erg cm}^{-2} \text{s}^{-1}$. Zand et al. (2002) estimate the distance to XTE J1908+094 to be greater than 3 kpc based on its X-ray flux, while Chaty et al. 2006 constrain the distance to be between 3 and 10 kpc based on optical measurements. Using the flux measurement at the soft-to-hard X-ray state transition, Curran et al. (2015) estimated the source distance to be between 4.8–8.3 kpc for a 3 $M_{\odot}$ black hole and between 7.8 and 13.6 kpc for an 8 $M_{\odot}$ black hole. Using these constraints on the distance to XTE J1908+098 and the mass of the central black hole, even during the highest flux phase (high hardness state of 2013 observation), the luminosity of the system is at most ~2.9% of its Eddington luminosity when computed using the 3–79 keV flux or ~14.4% when using the 0.5–100 keV flux. This maximal fraction is obtained for a 3 $M_{\odot}$ black hole at a distance of 8.3 kpc or a 8 $M_{\odot}$ black hole at 13.6 kpc. Even for the most conservative parameter combination that yields the maximum Eddington fraction, the value is well below the requirement of numerical simulations of $L/L_{\text{Edd}} \le 30\%$ for which the disk obeys the test particle ISCO (see, e.g., Fragile et al. 2018). This indicates that the assumption of inner accretion disk extending all the way to the size of the ISCO radius is reasonable and that the inner disk regions are gravity dominated. During the faintest part of the observations (low hardness state of 2019 observation), the luminosity of the system is greater than ~0.3% $L_{\text{Edd}}$ using the 3–79 keV flux or ~1.9% when using the 0.5–100 keV flux. These lower limits were obtained for a 3 $M_{\odot}$ black hole at a distance of 4.8 kpc or a 8 $M_{\odot}$ black hole at 7.8 kpc. These lower limits are indicative of a soft-intermediate state, transitioning into a hard state.

The 2013 data were previously analyzed by Tao et al. (2015) and Zhang et al. (2015). Both papers indicated the presence of a relativistic broadened Fe K line, but the attempts to measure the black hole spin using relativistic reflection were unsuccessful. By fitting the 2013 and 2019 observations jointly, the effects of relativistic reflection become more discernible from the disk continuum radiation. Previous works rebinned the spectra requiring at least 50 counts per bin, while in this work, we used the optimal binning algorithm described by Kaastra & Bleeker (2016). This binning scheme maintained the high signal-to-noise ratio at low

energies, while allowing multiple bins at high energies. Therefore, for relatively low flux observations such as the ones analyzed in this paper, the optimal binning scheme allows the models to both capture the details of the strong, broadened Fe K line and accurately track the shape of the Compton hump at energies above ~20 keV. The recent improvement in the NuSTAR spectra extraction pipeline now accounts for the tear in the FPMA thermal blanket that used to lead to a flux difference between the two NuSTAR sensors at low energies. All these factors in addition to the evolution of the *relxill* family of models in the recent times have encouraged pushing the potential of NuSTAR observations to new frontiers and obtaining spin measurements from spectra that were previously considered unable to allow the placement of constraints on the black hole spin.

The *relxill* variants adopting the lamp post geometry produce consistent parameter constraints with the models that make no prior coronal geometry assumption, while maintaining a similar quality of the fit. For the spin measurement of EXO 1846-031, Draghis et al. (2020) found that *relxilllp* fails to measure the black hole spin due to the model requiring a large coronal height. The observation of EXO 1846-031 on which the measurement was made occurred when the source had a flux ~10 times the flux of XTE J1908+094 during the 2013 observation and ~25 times larger than during the 2019 observation. For lower signal-to-noise spectra, the lamp post models are able to produce similar quality fits as other model variants, and are preferred in terms of DIC due to their relative simplicity owed to a smaller number of free parameters. This analysis shows that despite the difference in assumption about the corona geometry, in relatively low signal-to-noise spectra the lamp post models are able to place similar spin constraints, despite not fully covering the same space of physical properties. For higher quality spectra such as the one presented in Draghis et al. (2020), the data are additionally able to place constraints on the geometry of the corona and differentiate between lamp post and non-lamp post models.

Using BeppoSAX spectra, Miller et al. (2009) estimated the spin of XTE J1908+094 to be $a = 0.75 \pm 0.09$. This analysis predicts the spin of XTE J1908+094 to be $a = 0.55_{-0.45}^{+0.29}$. The two are formally consistent, but it is worth examining why the agreement is not better and why the errors are larger with better data. First, with the increase in energy resolution E/dE with NuSTAR over BeppoSAX MECS by a factor of nearly 2 at 6 keV and a factor of 4 at 10 keV, both the line profile and the shape of the continuum around the Fe K line can be more accurately constrained, reducing the possibility that the fit becomes stuck in a local minimum of $\chi^2$ in the parameter space as opposed to the global minimum, with falsely small errors. Second, in our fits to the NuSTAR observations, we attempted to understand possible sources of systematic uncertainty by exploring the entire set of physical phenomena covered by the *relxill* family of models. We concluded that the choice of model does not have an important impact on the spin measurement, as the predictions of all models are consistent with each other. However, the advances made in reflection models to include a wider range of physical effects has necessarily led to an increase in the number of free parameters in the models, leading to uncertainties that could not be probed in the past. By including a broader set of physics, the *accuracy* of the measurements has likely increased. It should also be noted that the error scans in this work are more complete and conservative. By running an MCMC analysis on the entire set
of parameters, the effects of all parameters are fully explored and incorporated into the uncertainty of the spin measurement. The larger and more physical uncertainty in the spin measurement has the potential to aid in eventually reconciling the spin distribution in X-ray binaries and BBH.

The spin distribution of black holes in low-mass X-ray binaries (LMXBs) determined from relativistic reflection measurements has been estimated to have a mean around $a = 0.66$ (Miller & Miller 2015). Based on the first two observing runs of Advanced LIGO and Advanced Virgo, it was estimated that half of the black holes in BBH systems have spins less than 0.27, and 90% have spins smaller than 0.55 (Abbott et al. 2019a). While the posterior distribution of spin measurements from GW signals from BBH mergers is still strongly influenced by the assumed prior distribution, a better understanding of the spin distribution of black holes in X-ray binaries will significantly aid the effort of breaking the degeneracies between parameters.

While works such as Kammoun et al. (2018) have shown the robustness of supermassive black hole (SMBH) spin measurements in AGNs using relativistic reflection, it is important to consider the effects of possible observational biases. One example is highlighted by Vasudevan et al. (2016), who argue that due to a selection effect, high spin SMBHs are more likely to be detected and therefore measured. This effect is likely to also translate to black holes in X-ray binaries. Additionally, it was argued by Jonker et al. (2021) that the mass range in the observed Galactic X-ray binaries is biased toward low-mass black holes. The spins of black holes of different masses are likely to both be formed at different values and be influenced differently by evolution effects. Additionally, using numerical GR MHD simulations of accretion, Król & Janiuk (2021) have shown that the accretion rate and black hole mass and spin evolution are strongly influenced by rotation dynamics of the infalling gas flow. Nevertheless, a black hole with an initial spin of $a=0$ must accrete at least half its mass to reach a spin parameter of $a = 0.84$ (Bardeen 1970), and since stellar-mass black holes in LMXBs are more massive than their companion stars, it is unlikely that stellar-mass black holes can significantly alter their spin parameters on relatively short timescales.

The best performing model measures an inclination of the inner accretion disk of XTE J1908+094 of $27 \pm 3$ degrees, with all models predicting inclinations <40°. By contrast, Rushton et al. (2017) measured an inclination larger than 79° by analyzing the proper motion of expanding radio jets in this source under the assumption of two symmetric jets propagating ballistically outwards, with subsequent evolution dictated by light travel time delays and relativistic boosting. Panel (a) in Figure 6 shows the unfolded spectrum when fit with model 6, but requiring the inclination to be larger than 79°. The dashed lines represent the contributions of the *diskbb* component of the model, while the dotted line shows the contribution of the *relxilllp* component. The $\chi^2$ becomes worse by ~85 and the residuals (panel (b) of Figure 6) show that the large inclination model predicts an emission excess on the blue wing of the Fe K line.

The residuals of the high inclination model shown in panel (b) in Figure 6 could be interpreted as the absorption from a high velocity disk wind. If we introduce a Gaussian component to our model with a negative amplitude and link its parameters within the 2013 and 2019 observation, the best-fit $\chi^2$ improves,

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Panel (a) shows the unfolded spectrum of XTE J1908+094 when fit with model 6, requiring an inclination >79°. The colors of the spectra match those described in Figure 3. The dashed line shows the contribution of the *diskbb* component of the model, while the dotted line shows the contribution of the *relxilllp* component. Panels (b) and (c) show the residuals in terms of $\sigma$ of the fits using the high inclination and free inclination models, respectively. The $\chi^2$ of the model requiring a high inclination is ~85 worse than the model that allows a free inclination. Note the increased residuals at ~7.4 keV (indicated by the vertical gray line) in panel (b) when compared to panel (c), suggesting that the high inclination model overpredicts the shape of the blue wing of the Fe K line. Panel (e) is a copy of panel (d) in Figure 3.</figcaption></figure>

but is still worse than the case of the free inclination model by 21, for six additional parameters. In order for the 6.9 keV absorption line of H-like Fe XXVI to be blueshifted to the measured line center of the Gaussian absorption component at 7.4 keV, the disk wind would have to travel at $v = 0.07$ c. An outflow with that velocity would shift the 8.3 keV line of He-like Fe XXV to 8.9 keV. The absence of this feature in the residuals of the model would be indicative of very high ionization of the gas.

Disk winds are typically observed in sources that are viewed at a high inclination. Sharp flux dips are even more typical of sources viewed at high inclination. The *lack* of dips in XTE J1908+094 and the requirement of a disk wind (of the same high velocity during the 2013 and 2019 outbursts) that perfectly balances the blue wing of the relativistic Fe K line represent a fortuitous combination of physical effects needed to explain a statistically disfavoured model. The low inclination model is then preferred both in terms of simplicity of the model and in terms of statistical significance, indicating that the large inclination determined by Rushton et al. (2017) was likely influenced by interactions of the radio jets with the interstellar medium. Rushton et al. (2017) note that the velocities implied by the jets do not match the flux differences between them, and
that the techniques used to measure the inclination of jets in sources like GRS 1915+105 may not hold in XTE J1908+094. In the future, simultaneous measurements of disk winds and relativistic reflection will be possible by combining NuSTAR observations with XRISM (Tashiro et al. 2018).

We thank James Miller-Jones and Elena Gallo for helpful discussions on the intriguing radio observations of XTE J1908+094. We thank Fiona Harrison, Karl Forster, and the NuSTAR team for executing observations of this source. The authors also acknowledge Brandon Case and Linda Hudson for computer support that made this work possible during COVID restrictions. We also thank the anonymous referee for their helpful comments and suggestions, which have clarified and

improved this paper. JAT acknowledges partial support under NASA NuSTAR Guest Investigator grant 80NSSC20K0644. ESK acknowledges financial support from the Centre National d'Etudes Spatiales (CNES).

## Appendix A
Tested Models

Table 3 shows the complete list of the 19 models tested in this analysis, together with the *relxill* variant they adopt, a brief summary of their assumptions, the best-fit $\chi^2$, number of degrees of freedom $\nu$, reduced $\chi^2$, computed DIC based on the MCMC runs, and the model rank in terms of lowest DIC.

<table><caption>**Table 3**
List of the 19 Tested Models</caption><thead><tr><th>Model</th><th>relxill Component</th><th>Assumptions</th><th>χ<sup>2</sup></th><th>ν</th><th>χ<sup>2</sup>/ν</th><th>DIC</th><th>Rank (by DIC)</th></tr></thead><tbody><tr><td>1</td><td>relxilllp</td><td>linked diskbb norm, free h</td><td>1726.39</td><td>1605</td><td>1.076</td><td>1810.25</td><td>13</td></tr><tr><td>2</td><td>relxilllp</td><td>linked diskbb norm, linked h within observation</td><td>1731.16</td><td>1607</td><td>1.077</td><td>1787.72</td><td>8</td></tr><tr><td>3</td><td>relxilllp</td><td>linked diskbb norm, linked all h</td><td>1779.16</td><td>1608</td><td>1.106</td><td>1846.02</td><td>17</td></tr><tr><td>4</td><td>relxilllp</td><td>linked diskbb norm, T, and h within observation</td><td>1728.36</td><td>1608</td><td>1.075</td><td>1797.4</td><td>11</td></tr><tr><td>5</td><td>relxilllp</td><td>free diskbb norm, free h</td><td>1707.56</td><td>1602</td><td>1.066</td><td>1784.94</td><td>6</td></tr><tr><td>6*</td><td>relxilllp</td><td>free diskbb norm, linked h within observation</td><td>1707.44</td><td>1604</td><td>1.064</td><td>1762.7</td><td>1</td></tr><tr><td>7</td><td>relxilllp</td><td>free diskbb norm, linked all h</td><td>1720.16</td><td>1605</td><td>1.072</td><td>1774.06</td><td>2</td></tr><tr><td>8</td><td>relxill</td><td>free diskbb norm, free q1,q1,r_br</td><td>1709.62</td><td>1594</td><td>1.073</td><td>1811.62</td><td>14</td></tr><tr><td>9</td><td>relxill</td><td>r_br fixed at 12</td><td>1706.26</td><td>1598</td><td>1.068</td><td>1786.95</td><td>7</td></tr><tr><td>10</td><td>relxill</td><td>linked diskbb norm, free q1,q2,r_br</td><td>1714.83</td><td>1597</td><td>1.074</td><td>1821.33</td><td>16</td></tr><tr><td>11</td><td>relxillD</td><td>N = 15 free diskbb norm, free q1,q2,r_br</td><td>1708.45</td><td>1598</td><td>1.069</td><td>1800.79</td><td>12</td></tr><tr><td>12</td><td>relxillD</td><td>N = 15, r_br = 12</td><td>1710.35</td><td>1602</td><td>1.068</td><td>1781.03</td><td>4</td></tr><tr><td>13</td><td>relxillD</td><td>N = 15 linked diskbb norm, free q1,q2,r_br</td><td>1724.13</td><td>1601</td><td>1.077</td><td>1866.92</td><td>18</td></tr><tr><td>14</td><td>relxillD</td><td>N = 17 free diskbb norm, free q1,q2,r_br</td><td>1708.59</td><td>1598</td><td>1.069</td><td>1818.26</td><td>15</td></tr><tr><td>15</td><td>relxillD</td><td>N = 17, r_br = 12</td><td>1710.53</td><td>1602</td><td>1.068</td><td>1783.98</td><td>5</td></tr><tr><td>16</td><td>relxillD</td><td>N = 17 linked diskbb norm, free q1,q2,r_br</td><td>1723.25</td><td>1601</td><td>1.076</td><td>1880.98</td><td>19</td></tr><tr><td>17</td><td>relxillD</td><td>N = 19 free diskbb norm, free q1,q2,r_br</td><td>1705.02</td><td>1598</td><td>1.067</td><td>1789.01</td><td>9</td></tr><tr><td>18*</td><td>relxillD</td><td>N = 19, r_br = 12</td><td>1705.44</td><td>1602</td><td>1.065</td><td>1778.81</td><td>3</td></tr><tr><td>19</td><td>relxillD</td><td>N = 19 linked diskbb norm, free q1,q2,r_br</td><td>1711.81</td><td>1601</td><td>1.069</td><td>1789.51</td><td>10</td></tr></tbody></table>

Note. Models 6 and 18 (marked by *) are the best performing models in terms of DIC among those adopting the lamp post geometry and those making no prior coronal geometry assumptions, respectively. The two models are discussed in detail in this work. Model 18 performs best out of the models using *relxill* or *relxilllp*, but ranks only third among all tested models.
**Appendix B**
**Corner Plot for Model 6**

Figure 7 shows the corner plot of a few parameters of interest in model 6. The parameters presented are the equivalent hydrogen column along the line of sight $N_H$ in

units of $10^{22}$ cm$^{-2}$, the height of the corona in the lamp post geometry in the 2013 and 2019 observations in units of $r_g$, the black hole spin $a$, inclination of the inner accretion disk in degrees, the Fe abundance $A_{\text{Fe}}$ in solar units, and the $\chi^2$ distribution.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Corner plot of some parameters of interest in model 6. The parameters shown are the interstellar hydrogen column $N_H$ in units of $10^{22}$ cm$^{-2}$, the height of the corona in the lamp post geometry in the 2013 and 2019 observations in units of $r_g$, the black hole spin $a$, inclination of the inner accretion disk in degrees, the Fe abundance $A_{\text{Fe}}$ in solar units, and the $\chi^2$ distribution.</figcaption></figure>
ORCID iDs

Paul A. Draghis https://orcid.org/0000-0002-2218-2306

Abderahmen Zoghbi https://orcid.org/0000-0002-0572-9613

Elias S. Kammoun https://orcid.org/0000-0002-0273-218X

Mark T. Reynolds https://orcid.org/0000-0003-1621-9392

John A. Tomsick https://orcid.org/0000-0001-5506-9855

References

Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019a, ApJL, 882, L24

Abbott, B. P. & LIGO Scientific Collaboration 2019b, PhRvX, 9, 031040

Abbott, R. & LIGO Scientific Collaboration 2021, PhRvX, 11, 021053

Arnaud, K. A. 1996, in ASP Conf. Ser., 101, Astronomical Data Analysis Software and Systems V, ed. G. H. Jacoby & J. Barnes (San Francisco, CA: ASP), 17

Banerjee, S., Chakraborty, C., & Bhattacharyya, S. 2019, MNRAS, 487, 3488

Bardeen, J. M. 1970, Natur, 226, 64

Bardeen, J. M., & Petterson, J. A. 1975, ApJL, 195, L65

Bardeen, J. M., Press, W. H., & Teukolsky, S. A. 1972, ApJ, 178, 347

Burrows, D. N., Kennea, J. A., Abbey, A. F., et al. 2007, Proc. SPIE, 6686, 668607

Chatterjee, D., Jana, A., Chatterjee, K., et al. 2021, Galax., 9, 25

Chaty, S., Mignani, R. P., & Israel, G. L. 2006, MNRAS, 365, 1387

Curran, P. A., Miller-Jones, J. C. A., Rushton, A. P., et al. 2015, MNRAS, 451, 3975

Dauser, T., Garcia, J., Parker, M. L., Fabian, A. C., & Wilms, J. 2014, MNRAS, 444, L100

Dauser, T., García, J., Walton, D. J., et al. 2016a, A&A, 590, A76

Dauser, T., García, J., & Wilms, J. 2016b, AN, 337, 362

Draghis, P. A., Miller, J. M., Cackett, E. M., et al. 2020, ApJ, 900, 78

El-Batal, A. M., Miller, J. M., Reynolds, M. T., et al. 2016, ApJL, 826, L12

Fabian, A. C., Lohfink, A., Kara, E., et al. 2015, MNRAS, 451, 4375

Fragile, P. C., Etheridge, S. M., Anninos, P., Mishra, B., & Klužniak, W. 2018, ApJ, 857, 1

García, J., Dauser, T., Lohfink, A., et al. 2014, ApJ, 782, 76

Goodman, J., & Weare, J. 2010, CAMCS, 5, 65

Harrison, F. A., Craig, W. W., Christensen, F. E., et al. 2013, ApJ, 770, 103

in't Zand, J. J. M., Miller, J. M., Oosterbroek, T., & Parmar, A. N. 2002, A&A, 394, 553

Jonker, P. G., Kaur, K., Stone, N., & Torres, M. A. P. 2021, arXiv:2104.03596

Kaastra, J. S., & Bleeker, J. A. M. 2016, A&A, 587, A151

Kammoun, E. S., Nardini, E., & Risaliti, G. 2018, A&A, 614, A44

King, A. L., Walton, D. J., Miller, J. M., et al. 2014, ApJ, 784, L2

Król, D. Ł., & Janiuk, A. 2021, ApJ, 912, 132

Lee, J. C., Fabian, A. C., Reynolds, C. S., Brandt, W. N., & Iwasawa, K. 2000, MNRAS, 318, 857

Li, L.-X., Zimmerman, E. R., Narayan, R., & McClintock, J. E. 2005, ApJS, 157, 335

Liska, M., Hesp, C., Tchekhovskoy, A., et al. 2020, MNRAS, 507, 983

Ludlam, R. M., Remillard, R., Homan, J., et al. 2019, ATel, 12652, 1

Madsen, K. K., Grefenstette, B. W., Pike, S., et al. 2020, arXiv:2005.00569

McClintock, J. E., Shafee, R., Narayan, R., et al. 2006, ApJ, 652, 518

Merloni, A., & Fabian, A. C. 2001, MNRAS, 321, 549

Miller, J. M. 2007, ARA&A, 45, 441

Miller, J. M., Reynolds, C. S., Fabian, A. C., Miniutti, G., & Gallo, L. C. 2009, ApJ, 697, 900

Miller, M. C., & Miller, J. M. 2015, PhR, 548, 1

Mitsuda, K., Inoue, H., Koyama, K., et al. 1984, PASJ, 36, 741

Morningstar, W. R., Miller, J. M., Reis, R. C., & Ebisawa, K. 2014, ApJL, 784, L18

Nealon, R., Price, D. J., & Nixon, C. J. 2015, MNRAS, 448, 1526

Negoro, H., Suzuki, K., Ueno, J. S. S., et al. 2013, ATel, 5549, 1

Novikov, I. D., & Thorne, K. S. 1973, in Black Holes (Les Astres Occlus), ed. C. DeWitt & B. DeWitt (New York: Gordon & Breach), 343

Pürrer, M., Hannam, M., & Ohme, F. 2016, PhRvD, 93, 084042

Reynolds, C. S. 2021, 43rd COSPAR Scientific Assembly (Sydney, Australia)

Reynolds, C. S., & Fabian, A. C. 2008, ApJ, 675, 1048

Ross, R. R., & Fabian, A. C. 2005, MNRAS, 358, 211

Rushton, A. P., Miller-Jones, J. C. A., Curran, P. A., et al. 2017, MNRAS, 468, 2788

Salvesen, G., & Miller, J. M. 2021, MNRAS, 500, 3640

Shafee, R., Narayan, R., & McClintock, J. E. 2008, ApJ, 676, 549

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33

Shimura, T., & Takahara, F. 1995, ApJ, 445, 780

Sobczak, G. J., McClintock, J. E., Remillard, R. A., et al. 2000, ApJ, 544, 993

Sokal, A. 1996, in Monte Carlo Methods in Statistical Mechanics: Foundations and New Algorithms, ed. C. DeWitt-Morette, P. Cartier, & A. Folacci, Vol. 361 (Boston, MA: Springer)

Spiegelhalter, D. J., Best, N. G., Carlin, B. P., & Van Der Linde, A. 2002, J. R. Stat. Soc. Series B Stat. Methodol, 64, 583

Steiner, J. F., McClintock, J. E., Orosz, J. A., et al. 2014, ApJL, 793, L29

Tao, L., Tomsick, J. A., Walton, D. J., et al. 2015, ApJ, 811, 51

Tashiro, G., Maejima, H., Toda, K., et al. 2018, Proc. SPIE, 10699, 1069922

Tiwari, V., Fairhurst, S., & Hannam, M. 2018, ApJ, 868, 140

Tomsick, J. A., Parker, M. L., García, J. A., et al. 2018, ApJ, 855, 3

Vasudevan, R. V., Fabian, A. C., Reynolds, C. S., et al. 2016, MNRAS, 458, 2012

Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996, ApJ, 465, 487

Vitale, S., Gerosa, D., Haster, C.-J., Chatziioannou, K., & Zimmerman, A. 2017, PhRvL, 119, 251103

Wilkins, D. R., & Fabian, A. C. 2012, MNRAS, 424, 1284

Wilms, J., Allen, A., & McCray, R. 2000, ApJ, 542, 914

Wilms, J., Juett, A., Schulz, N., & Nowak, M. 2006, BAAS, 38, 376

Woods, P. M., Kouveliotou, C., Finger, M. H., et al. 2002, IAU Circ.No. 7856

Xu, Y., Harrison, F. A., Tomsick, J. A., et al. 2020, ApJ, 893, 30

Zevin, M., Berry, C. P. L., Coughlin, S., Chatziioannou, K., & Vitale, S. 2020, ApJL, 899, L17

Zhang, L., Chen, L., Qu, J.-l., Bu, Q.-c., & Zhang, W. 2015, ApJ, 813, 90