A Dependable Distance Estimator to Black Hole Low-Mass X-ray Binaries

Black Hole Low Mass X-ray Binaries (BH-LMXBs) are excellent observational laboratories for studying many open questions in accretion physics. However, determining the physical properties of BH-LMXBs necessitates knowing their distances. With the increased discovery rate of BH-LMXBs, many canonical methods cannot produce accurate distance estimates at the desired pace. In this study, we develop a versatile statistical framework to obtain robust distance estimates soon after discovery. Our framework builds on previous methods where the soft spectral state and the soft-to-hard spectral state transitions, typically present in an outbursting BH-LMXB, are used to place constraints on mass and distance. We further develop the traditional framework by incorporating general relativistic corrections, accounting for spectral/physical parameter uncertainties, and employing assumptions grounded in current theoretical and observational knowledge. We tested our framework by analyzing a sample of 50 BH-LMXB sources using X-ray spectral data from the \swift/XRT, \maxi/GSC, and \rxte/PCA missions. By modeling their spectra, we applied our framework to 26 sources from the 50. Comparison of our estimated distances to previous distance estimates indicates that our findings are dependable and in agreement with the accurate estimates obtained through parallax and H\,{\sc i} absorption methods. Investigating the accuracy of our constraints, we have found that estimates obtained using both the soft and transition spectral information have a median uncertainty (1$\sigma$) of 20\%, while estimates obtained using only the soft spectral state spectrum have a median uncertainty (1$\sigma$) of around 50\%. Furthermore, we have found no instrument-specific biases.


INTRODUCTION
It is unequivocal that the accretion mechanisms around black holes are not yet fully understood.Black hole X-ray binaries (BHXBs) provide ideal observational laboratories to study the accretion process.In particular, transient low mass X-ray binaries (BH-LMXBs) show variations in their accretion rate at a time scale of days to a few years (Fender & Gallo 2014;Tetarenko et al. 2016).This renders investigations into the underlying physics of accretion much more feasible compared to supermassive black hole accretion systems present in Active Galactic Nuclei (AGN; Silk & Rees 1998).
Depending on the mass of the companion star and the variability of the emission, BHXBs are divided into (1) Persistent, typically highmass X-ray binaries (BH-HMXBs) where the companion star mass is comparable to or greater than the primary BH mass and mass is always being transferred, (2) Transient, which are found to be mostly low mass X-ray binaries (BH-LMXBs), where the companion star mass is ≤  ⊙ and the mass accretion rate onto the black hole is typically low.However, these transient systems undergo recurrent outbursts in which an increased accretion rate gives rise to a bright emission followed by an exponential decay back to their usual faint state, typically known as quiescence (Corral-Santana et al. 2016;Tetarenko et al. 2016).
★ E-mail: youssefabdulghani@montana.edu There are, however, notable exceptions to this picture of black hole X-ray binaries, one such exception is 4U 1957+11.This source has been studied multiple times since its discovery in 1974 (Giacconi et al. 1974).Evidence from most studies strongly suggests that it is a persistent BH-LMXB observed to be in the soft state (e.g., Wĳnands et al. 2002;Gomez et al. 2015;Barillier et al. 2023) The transient behavior of BH-LMXBs is thought to be driven by instabilities in the accretion disk (Dubus et al. 2001;Remillard & McClintock 2006) leading to large fluctuations of the accretion rate.During a typical outburst, a rapid increase in the accretion rate (timescale of days to months) causes characteristic changes to the accretion geometry and emitted radiation (e.g.Remillard & McClintock 2006;Fender et al. 2004).
Furthermore, transient BH-LMXBs are observed to evolve through a series of different spectral states during typical outbursts.This evolution of spectral states is thought to be explained by a hot inner advection-dominated accretion flow (ADAF) and a truncated accretion disk (Esin et al. 1997).In this picture, the accretion rate is low initially before the outburst.During that time, the accretion disk does not extend all the way to the black hole but is recessed.As a result, the radiation emitted by the accretion disk during this time is weak.At low accretion rates, the "hard" X-ray emission from the inner ADAF dominates over the disk emission, and this state of the BH-LMXB is therefore referred to as the hard state.As the accretion rate increases, the accretion disk starts to move closer to the black hole.The disk contribution to the X-ray spectrum increases.The BH-LMXB moves into the intermediate state.Then, as the accretion rate gets closer to the maximum possible rate, the disk extends until the innermost stable orbit (ISCO), and the accretion disk spectrum completely dominates the X-ray spectrum.This high luminosity state is called the soft state, where the low energy (soft X-ray) disk black body emission dominates.Finally, as the episodic increased accretion event comes to an end, the accreting BH-LMXB fades via the intermediate state back into the quiescence (for review see Done et al. 2007, and reference therin).Thus, in this framework, the observed radiation is produced primarily via two processes: 1) the thermal emission from the accretion disk itself (Shakura & Sunyaev 1973;Mitsuda et al. 1984) and 2) the up-scattered emission from the inner ADAF.This up-scattered emission is created by some of the primary disk photons scattering off the hot electrons in the inner ADAF and gaining energy, resulting in a power-law-like spectrum towards higher (harder) X-ray energies (Thorne & Price 1975).
Although the current models for accretion can explain some of the observed behaviour of BH-LMXBs, much remains unknown.Moreover, intrinsic observables such as accretion geometry, jet power, luminosity, and magnetic field strength are needed to refine existing frameworks, such as the inner ADAF and the truncated disk model, and to test new theoretical models.Crucially, getting constraints on such intrinsic properties is often hampered by the poor quality of the distance estimate of these systems.
The most accurate distances to BH-LMXBs can obtained from a model-independent parallax measurement to the companion star/jet (e.g.Miller-Jones et al. 2009;Gandhi et al. 2019;Atri et al. 2020;Arnason et al. 2021).However, this approach cannot be used if the companion star/jet is too faint, too far or there are not enough observations available to permit the parallax calculation.Other methods include using interstellar extinction (Jonker & Nelemans 2004), analyzing the spectrum of the companion star (Jonker & Nelemans 2004), using the 21-cm line of neutral hydrogen (e.g.Dhawan et al. 2000;Chauhan et al. 2021), modelling the X-ray dust scattering halo radial profiles (e.g.Trümper & Schönfelder 1973;Thompson & Rothschild 2009;Kalemci et al. 2018), or measuring the proper motion of jet ejecta (e.g.Mirabel & Rodríguez 1994).All these methods also suffer from limitations, such as the requirement that the companion star is bright enough for spectral analysis, the existence of bipolar jets and/or that the inclination of the source is known.Their constraints also often suffer from high uncertainties (Jonker & Nelemans 2004).
As of December 2023, the catalogue of stellar-mass black holes in X-ray binaries 1 produced by Corral-Santana et al. (2016) (Black-CAT), lists a total of 70 observed transient sources where the compact object is either a dynamically confirmed BH or a strong BH candidate.From these 70 sources, only 28 have good distance estimates, and 14 have poorly constrained distances (error > 50% or no lower/upper bounds).While 28 have no distance estimate to the best of our knowledge.It is estimated that the number of BH-LMXBs in our Galaxy is of the order of thousands (e.g.Romani 1998;Kiel & Hurley 2006).Moreover, there has been an exponential increase in the detection of transient sources in the past two decades (Corral-Santana et al. 2016), driven by recent missions/instruments such as the Rossi X-ray Timing Explorer (RXTE; Bradt et al. 1993), Neil Gehrels Swift Observatory/Burst Alert Telescope (Swift/BAT; Barthelmy et al. 2005) and the Monitor of All-sky X-ray Image (MAXI; Matsuoka et al. 2009).This rate is only expected to keep increasing.Keeping up with the ever-growing BH-LMXBs discoveries will require a simple and generally applicable distance estimation method.
While various BH-LMXB systems (e.g., GRO J0422+32, XTE J1118+480, some outbursts of GX 339-4) have been found to exhibit only the hard state (Alabarta et al. 2021), the majority of the current sample of BH-LMXBs (about 70% of the observed outbursts, according to Alabarta et al. 2021) are detected in the X-ray band when they undergo a "canonical" outburst.Here, we define a "canonical" outburst as one where the hard, soft, and intermediate states are observed.Consequently, a strong candidate for a distance estimation method is one that utilizes the modeling results of the X-ray spectra from specific spectral states typically present in any canonical BH-LMXB outburst, namely, the soft state and the soft-to-hard transition (e.g., Remillard & McClintock 2006;Kalemci et al. 2022).Moreover, throughout the transient BH-LMXB literature, this method has been previously used (e.g.Nakahira et al. 2012Nakahira et al. , 2014;;Oda et al. 2019;Tominaga et al. 2020) however, no systematic study of the robustness of this method was ever conducted.
In this study, we further develop this method, quantify its uncertainties, and implement it in a framework that allows for its application to any BH-LMXB exhibiting the typical outburst pattern.Our paper is organized as follows: Sections 2 and 3 describe how we selected the sample, the spectral data, and its analysis.Subsequently, in Section 4, we detail the spectral models used and the results of modeling the spectral data.In Section 5, we establish the statistical framework used to obtain our distance and mass estimates, presenting its results.Section 6 assesses the accuracy of our method and discusses the overall resulting mass and distance distributions, along with limitations and biases.Finally, Section 7 provides our conclusions.

SAMPLE & DATA SELECTION
Out of the 70 transient black hole sources observed as of December 2023, we selected all sources that underwent at least one typical outburst according to the catalogues of Tetarenko et al. (2016), Corral-Santana et al. (2016) 2 , and Alabarta et al. (2021).This yields an initial sample of 50 potential BH-LMXB sources.For this initial sample of targets, we created separate light curves of every catalogued outburst using available archival data from MAXI/GSC, Swift/XRT and/or RXTE/PCA.The light curve creation and analysis are described below (Sections 3.1 and 3.2).Using the light curve analysis, we could further limit our sample to sources with an identifiable soft state in at least one of the observed outbursts.
Moreover, we excluded sources from our sample whose spectra were later found to be of low quality for our analysis.Namely, sources that had less than 10 data points after binning to a signal-to-noise ratio of 3 or sources that were unable to obtain an X-ray spectral description with 2 red < 1.4 using our models.After applying these filtration criteria, we were left with a final sample of 26 unique sources for our analysis, with some sources having multiple outbursts or the same outburst detected in multiple instruments.
As we will see in Section 5, our distance estimation method is ideally based on the spectral modelling of the soft and the softto-hard transition state spectra during a typical outburst from the same instrument.However, a lower-quality distance estimate can be obtained from the soft state alone.In our final sample of 26 sources, the 18 sources with both soft state and transient state spectral information can be found in Table 1, and the 8 sources with only soft spectral information are listed in Table 2.
We have chosen to consider data from MAXI/GSC, Swift/XRT and RXTE/PCA for our study as these instruments have observed the majority of BH-LMXB outbursts.Additionally, MAXI/GSC and Swift/XRT continue to make new discoveries.We did not simultaneously model data or combine estimates coming from these different instruments for the same source/outburst.This allows us to assess any biases arising from using a specific detector since we have used data of the same source/outburst from multiple detectors whenever possible.

DATA REDUCTION & PREPARATION
This section describes how the MAXI/GSC, Swift/XRT, and RXTE/PCA data for this study have been reduced and how light curves and spectra of certain points in the hardness intensity diagram were created.
All locally downloaded data were reduced and analyzed using HEASoft v6.32.1 (Nasa High Energy Astrophysics Science Archive Research Center (Heasarc) 2014) with the latest CALDB files (as of Dec 2023) for each of the used instruments.

MAXI/GSC
For MAXI/GSC, data were obtained using the mxdownload_wget command and reduced using the mxproduct pipeline both provided in HEASoft from the MAXI team3 .Source and background regions were determined using the MAXI on-demand4 web tool (Nakahira et al. 2012).For each source, daily average spectra were then produced using the tool.Daily light curves in two bands, 2-4 keV and 4-20 keV, were created from the reduced daily spectra.

Swift/XRT
For those outbursts covered by Swift/XRT, light curves were obtained using the Swift/XRT web tool (Evans et al. 2007(Evans et al. , 2009)), which used HEASoft v6.29 at the time we produced the light curves.For each outburst, light curves in the energy range of 0.3-4 keV and 4-10 keV were produced with a binning by observation.

RXTE/PCA
Ready-analyzed RXTE/PCA light curves were downloaded using the HEAVENS5 tool, which is part of the INTEGRAL Science Data Centre (ISDC; Lubiński 2009;Walter et al. 2010).The 2-10 keV and 10-20 keV light curves were binned in 1-day bins and limited to counts from the Proportional Counter Unit (PCU) 2 detector and the top layer.

Spectral State Identification
After obtaining the light curves for each instrument, we calculate the hardness ratio (HR) with time for individual outbursts by dividing the number of counts in the harder band light curve by the softer band light curve counts in a given time bin.The soft and transition periods are then identified manually using the values of the daily HR.In particular, following the creation of the HR for each outburst and instrument, we first identify the entirety of the soft state period, the time when there is a consistently lower HR relative to the HR at the start of the outburst.The soft-to-hard transition period is then determined by identifying the first few days when the HR recovers from the consistently low values.This selection is performed by visual inspection of the HR plot.
In order to validate this manual identification method, we also ran a K-means clustering algorithm where we use the HR and intensity (total flux) as input features similar to a recent study by Sreehari & Nandi (2021).We use a total of 4 clusters corresponding to the hard state, hard-to-soft intermediate state, soft state, and the soft-to-hard intermediate state.We then let the algorithm assign each data point to one of these clusters.For most instances, there was a consensus between the manual identification process and the results derived from K-means clustering.However, there were a few cases where the manual identification did not align with the K-means clustering outcome.In these situations, we relied on the K-means clustering results to identify the appropriate periods of the soft state and the soft-to-hard transition.
Once the soft state and ideally the transition state dates are identified (Table A1), we then produce corresponding spectra for further analysis as described in the next sub-section.

MAXI/GSC
For MAXI, we have already produced daily spectra to obtain a light curve, and we utilize those for our soft and transition state spectra.
However, most of the individual daily spectra did not have the necessary statistics to provide meaningful spectral parameters.Thus, spectra with similar hardness ratios (Photons s −1 m −2 of 4-20 keV to 2-4 keV) were grouped together such that the combined spectrum has at least 2000 net (background-subtracted) photons.Similarly, we required a minimum of 2000 net counts for the soft-to-hard transition period.The final combined spectra were then binned to a signal-tonoise (S/N) ratio of 3.

Swift/XRT
For Swift/XRT, we downloaded the Windowed Timing (WT) mode observations matching the soft and transition state times and reduced the data locally using the xrtpipeline.To create a spectrum, the source region for non-piled-up observations was taken to be a circle with a radius   of 50 pixels centered on the source coordinates.While the background region was taken to be an annulus with an inner radius  ,in of 80 pixels and an outer radius  ,out of 120 pixels.Following the guidelines by the Swift team6 , both the source and background had their BACKSCAL keyword edited to 2  and  ,out −  ,in − 1, respectively.
For sources with 100 or more counts per second in the resulting spectrum, a pile-up correction was performed similar to method 2 described in (Mineo et al. 2006) but with a simpler approach.In our approach, we specify an annulus region centered on the source center.The starting inner radius is arbitrarily chosen, while an outer radius is chosen to be 100 arcseconds for all observations.We then increase or decrease the inner radius by several arcseconds until we have a spectrum with a count rate of less than 100.The background region shape for the piled-up sources is then taken to be the same as the source annulus but sufficiently away from the center.The spectra of the individual observations are much more statistically significant than the individual observations from MAXI, so no single observations were combined.All spectra were binned to an S/N ratio of 3.

RXTE/PCA
Data reduction of the observations matching the soft/transition states was done locally by first processing the data into usable Standard2 spectral products using the PCA Standard2 extraction pipeline provided by the RXTE team7 and included in HEASoft v6.32.1.We used the latest CALDB file released on May 1, 2023, which automatically adds a systematic error that is recommended by the RXTE team.The pipeline runs 4 steps: (1) data preparation, where basic filtering and data processing happen, (2) merging observations if using more than one, (3) making a good time filter which excludes bad data points, (4) extracting the final spectral product which is a "Standard2" data product.The Standard2 data products contain the energy and temporal information required for spectral modelling.The spectra were also binned to a S/N ratio of 3 and no observations were combined.

SPECTRAL ANALYSIS
To ultimately obtain a constraint for the distance to a BH-LMXB from the X-ray spectrum during an outburst.First, we need to model the spectrum during the soft state and, if possible, the soft-to-hard transition.
All spectral analyses presented here were performed in XSPEC (v12.13.1;Arnaud 1996).Solar abundance values were set to wilm (Wilms et al. 2000) while the cross-section table was set to vern (Verner et al. 1996).The built-in Markov Chain Monte Carlo (MCMC) was used to find the uncertainties in all variable spectral parameters.All MCMC chains used the Goodman-Weare (GW) algorithm and had a total length of 2 × 10 6 with a burn-in of 1 × 10 5 and 40 walkers.
In BH-LMXBs, the soft state and the soft-to-hard transition Xray spectrum in the energy range 0.5-20 keV can be described well to first order by combining two models: a multi-color disk black body model accounts for the accretion disk emission and a Compton up-scattered emission component (Remillard & McClintock 2006), which a power law can approximate.This intrinsic X-ray spectral model is then further altered by photoelectric absorption along the line of sight.Thus, in XSPEC we implemented this base model as: TBabs * (powerlaw+ezdiskbb) (Model S1), where the TBabs (Wilms et al. 2000) model accounts for the interstellar absorption, a phenomenological powerlaw models the Compton up-scattered photons, and the ezdiskbb (Zimmerman et al. 2005) models the disk's radiation.In this model set-up, we always fit for the photon index (Γ), power law normalization, the maximum temperature of the disk ( max ; assumed to be > 0.1 keV), and the disk normalization ( ezdiskbb ).
For the MAXI/GSC and RXTE/PCA spectra, the energy ranges used for spectral modelling were 2-20 keV and 3-20 keV, respectively.This leads to the data being unable to constrain the interstellar absorption column, ( H ), and the parameter was that kept fixed8 to the literature values for each source (see Table A1 for the used absorption columns).On the other hand, the energy range used for Swift/XRT spectra was 0.7-10 keV, which allowed for  H to be a free 9 parameter.
For the soft state spectra, we expect that the spectrum is very steep, and consequently, Γ was only permitted to vary between 1.7 and 3.0.The parameter relevant for the distance estimation from the spectral parameters of this state is the normalization of the disk component.On the other hand, we expect the spectrum to be harder for the transition period, and the photon index was confined to the lower range of 1.0-2.2.The parameter that will be relevant for distance estimation from the transition state spectrum is the 0.5-200 keV power-law flux, as discussed in the next section.
We note that for some RXTE/PCA sources, the spectrum contained a significant 6.4 keV iron K line that was affecting the spectral parameters.Therefore, for those sources, we added a Gaussian component (gauss) centered at 6.4 keV and a maximum width of 1 keV.The full model used for these sources is implemented in XSPEC as: TBabs * (powerlaw+ezdiskbb+gauss) (Model S2).
We require the unabsorbed power-law flux of the transition state in the 0.5-200 keV energy band for the distance estimation.The models above were modified to obtain a robust constraint of this flux from the convolution model cflux.We extended the response matrix range using the energies command, so the model calculation covers the entire 0.5-200 keV band.As the cflux flux value and the power law normalization would be degenerate in the modeling, we adjust the modeling set-up for the transition state.Namely, we freeze the normalization of the unabsorbed model power-law component and apply it as the following: TBabs * (cflux * powerlaw+ezdiskbb) (Model T1) or, if this does not yield to a satisfactory description, TBabs * (cflux * powerlaw+ezdiskbb+gauss) (Model T2).
Using the version of Model 1 or Model 2 corresponding to the state, as described above, we modelled a total of 66 spectra.More specifically, we fitted 39 soft-state spectra and 27 soft-to-hard transition spectra for a total number of 34 unique outbursts.
Our basic models 1 and 2, were able to describe these spectra well with  2 red ≤ 1.4 in all cases and the majority of fits having  2 red very close to 1.A summary of our spectral modeling results can be found in Table A1 in the appendix.
From our spectral modeling, we obtain the posterior distribution of the normalization of ezdiskbb ( ezdiskbb ) by performing MCMC analysis with Gaussian priors based on the best-fit model parameters for each of the soft state spectra.Similarly, for soft-to-hard transition fits, we obtain the posterior distribution of the power law flux in the 0.5-200 keV range by utilizing the cflux component.For computational efficiency, these resulting distributions were down-sampled from the initial length 2 × 10 6 (length of the MCMC chain) to 2 × 10 4 by averaging every 100 adjacent values.The MCMC distributions were then used as inputs for the statistical framework described in Section 5.3 to produce distance estimates for each source in our final sample.

DISTANCE CALCULATION
In the following, we discuss how the modelling results from the soft and the transition state spectra can be used to obtain a high-quality distance constraint.

The Basic Idea
A fairly simple distance estimation method based on only information from X-ray spectral modelling was used multiple times throughout the literature (see Nakahira et al. 2012Nakahira et al. , 2014;;Oda et al. 2019;Tominaga et al. 2020) to produce a quick constraint on the distance to a BH-LMXB in outburst.
This method is based on the same fundamental idea as many astronomical distance estimation methods: if we measure the flux and know what the intrinsic luminosity should be, we can find the distance.In the soft state, the X-ray spectrum of black hole X-ray binaries can be assumed to be disk-dominated, with only a small portion of the flux potentially originating from the power-law component.Therefore, we can get very accurate information about the total observed flux in this state from the accretion disk flux.In practice, the spectral modeling of the soft state spectrum does not return the accretion disk flux directly, but the accretion disk spectrum is modeled and then the disk normalization is used to derive the fluxdistance relation.In particular, the normalization of the multi-color disk black body accretion disk model, ezdiskbb, that we have used in our modeling is defined as (Zimmerman et al. 2005): where  col = 1.7 (see Section 6.4 for discussion about this value) is the color (spectral) hardening factor (Shimura & Takahara 1995),  is the distance to the source in kpc, and  is the disk inclination as observed in the sky (e.g. = 0 • is face-on).As evident from equation ( 1), the observed disk spectrum depends not only on the distance to the source but also on the inclination and inner radius of the accretion disk.Observational evidence indicates that this inner radius is nearly constant during the soft state and is thought to be equal to the theoretical innermost stable circular orbit (ISCO; e.g.Steiner et al. 2010).The location of the ISCO generally depends on the BH mass and spin.However, if we assume a non-rotating BH, then it only depends on mass, where  ISCO = 6 / 2 , with  being the black hole mass,  is the gravitational constant and  is the speed of light.As a result, the constraint of the normalization of the ezdiskbb model obtained from the spectral modelling can be used to produce a relation between the mass, the inclination, and the distance to the source.Ideally, we need to constrain the parameter space further to break some of the degeneracy between the distance, the mass, and the inclination for a given spin value.To do this, we make use of an empirical relation first found by Maccarone, T. J. ( 2003 2019), in which the bolometric luminosity of the state transition from the soft state to the hard state was typically observed to happen at 1-4 % of the Eddington luminosity ( Edd ).More specifically, Maccarone, T. J. (2003) andrecently Vahdat Motlagh et al. (2019) found that the (∼ 0.5 − 200 keV) power law flux in the soft-to-hard spectral state transition has been consistently observed to have a mean of 2% of the Eddington luminosity.A similar study by Tetarenko et al. (2016), where they considered not only the power law flux but also the emission of the disk, also found a small, mean percentage of ∼3%.From modeling the X-ray spectrum at the transition time, we are able to obtain the flux.Since the flux is related to the luminosity via the squared distance and the Eddington luminosity depends on mass, we obtain an additional relation between distance and mass.For instance, using the (1-4 %) range found by Maccarone, T. J. (2003): where  trans is the bolometric flux at the soft-to-hard transition,  Edd is the Eddington fraction which is 0.01-0.04,and we have used  Edd = 1.3 × 10 38   ⊙ erg/s and  = 4 2  (isotropic emission).So theoretically, one can obtain a constraint on the distance from the spectral modeling results and equations 1 and/or 2 after making assumptions for  and .However, this simple approach assumes a non-rotating black hole and ignores the specific statistical distribution of the spectral parameters and other relativistic effects, as discussed below.

Limb Darkening and Relativistic Corrections
In the previous section, limb-darkening or relativistic effects were not considered, yet these are important effects that influence the accuracy of this type of distance estimate.Despite their importance, applications of the method outlined in the previous section have mostly ignored these effects.In this study, we take these corrections into account by following the work done by Salvesen & Miller (2021) (see also Zhang et al. 1997) and re-writing the normalization of ezdiskbb as: where Υ() = 1 2 + 3 4 cos  is the limb darkening law assuming a disk atmosphere that is plane-parallel, semi-infinite and electron scattered (Chandrasekhar 1960),   (, ) is a correction factor that takes into account the relativistic effects on photon propagation, and    () is a correction factor that takes into account the relativistic radial structure of the accretion disk.We note that   (, ) depends on the black hole spin, , and the disk inclination, , while    () only depends on the black hole spin.These dependencies allow us to pre-calculate these correction factors for possible values of spin and inclination regardless of the properties of any particular source.In fact, these corrections were tabulated by Salvesen & Miller (2021) and made available in a public repository10 .In their work, the authors utilize the kerrbb model (Li et al. 2005) and numerical integration to calculate the flux, taking into account the relativistic effects.They then determine the correction factors by calculating the ratio of fluxes with and without these effects.In our work, we use their tabulated correction values to account for the relativistic effects.

Obtaining Distance/Mass Probability Distributions
Rather than just making assumptions using a few discrete values for , , , we use an approach that considers all of the uncertainties entering from the calculations and yields a probability distribution of the distance values for each source.
It starts with the MCMC spectral modeling error analysis of the soft state spectrum of a source, which produces a probability distribution of  ezdiskbb .We then obtain a joint mass/distance probability density function  , (, |soft) by calculating the distance using equation (3) over sampled , , and  values.Similarly, the spectral modelling followed by an MCMC error analysis of the source when it goes through the soft-to-hard transition produces a probability distribution of  trans .Once again, using an assumed distribution of black hole mass, as well as an Eddington fraction  Edd , we find another mass/distance probability density function  , (, |trans) using equation ( 2).Finally, we assume that the best estimate for the true mass/distance joint probability distribution is: where  is some normalization factor that makes the combined probability density function  , (, ) integrated over all possible distance and mass values equal to one.
We can find the marginal probability distribution of distance or mass by integrating over all possible values of the other variable.For example, to obtain the marginal probability density for the distance we integrate over all possible masses: Furthermore, equation ( 4) shows that unless  , (, |soft) and  , (, |trans) are both non-zero at all distances , we in fact gain new information on the mass probability density function.So, we obtain a refined mass distribution by requiring that the only valid distance values of a certain source are where both constraints have non-zero probability density values.This approach is analogous to Bayesian inference tools, where we obtain a posterior distribution based on the prior guess.
If only a constraint from the soft state spectrum exists, we follow the same procedure until we obtain a mass and distance probability distribution.As we do not have a second mass and distance probability, the distribution of the combination step does not take place.An overview diagram of the process is shown in Fig. 1.
As mentioned before, we need to use probability distributions for the inclination, black hole spin, Eddington fraction, and black hole mass.For the inclination, we assume an isotropic distribution in the range (0°, 84°) where the probability of observing a given inclination is given by   () = sin .We chose this specific inclination range since it is very unlikely to observe the disk spectrum at very high inclinations (it will be very faint).Additionally, the correction values   (, ) and    () in equation ( 3) are calculated by utilizing the kerrbb model in XSPEC which only allows for angles 0°to 85°.We note here that this inclination prior could be improved in two ways: universally, by accounting for observational selection effects (a topic beyond the present paper's scope) or, specifically, for individual sources, by considering a probable inclination distribution for the specific source.Notably, we draw attention to a method proposed by Muñoz-Darias et al. (2013) in which the shape of the HID track can be used as an estimator for inclination.This method is of particular relevance given that both relativistic effects and limb-darkening significantly alter the observed pattern.In future work, we aim to conduct a detailed systematic study where we use these HID properties and develop them into a method yielding inclination constraints.
For the black hole spin parameter, we assume a uniform distribution in the range 0.0 to 0.998, where all the spin values in that range are equally likely:   () = 1/(0.998− 0).This range of spins was chosen because of the rarity of BH-LMXBs that are estimated to have retrograde spin (Reynolds 2021).However, we investigate the impact of allowing for retrograde spin on our results in Section 6.4.
For the Eddington fraction at the transition point, we assume a uniform distribution in the range 0.01 to 0.04 according to Maccarone, T. J. ( 2003):   Edd (  Edd ) = 1/(0.04− 0.01).We chose to use a range of 1-4% as it captures most of the observed variations in the transition luminosities.We also note that we use only the power law flux following both Maccarone, T. J. (2003) andVahdat Motlagh et al. (2019).
Finally, for the black hole mass distribution, we assume a Gaussian prior centered at 7.8  ⊙ with a standard deviation of 1.2  ⊙ .This particular Gaussian distribution was selected based on the work of Özel et al. (2010), Farr et al. (2011), and Kreidberg et al. (2012).
In all our calculations of  , (, ), we use these distributions to sample 101, 100, 101, 1000 values for the inclination, spin, Eddington fraction, and black hole mass, respectively.As shown in the example in Fig. 2, we plot the marginal probability distributions of distance from both constraints (red/dashed for soft, blue/dotted for transition) and then obtain the best estimate of the marginal probability by multiplying both distributions (black/solid).Similarly, for the joint mass/distance distribution, we plot the 2D contours of the soft constraint (red/dashed contours) and the transition constraint (blue/dotted contours).We then multiply both of the joint probabilities to obtain the best estimate of the joint mass/distance probability distribution (black/solid contours).

Sources with Multiple Outbursts
Some sources in our sample have been observed during multiple outbursts.To obtain the best estimate for the true mass/distance joint probability distribution, we first find mass/distance probability distri- The joint mass/distance probability distribution.The red region corresponds to the probability distribution from only the soft state constraint, while the blue region shows the probability distribution from only the transition luminosity.The grey region visualizes the combined probability distribution that results when both distributions are assumed to be true, equation ( 4).The contour levels shown are 1, 2, and 3.Right: The marginalized distance probability distribution.The red line corresponds to the probability distribution from only the soft state constraint, while the blue line indicates the probability distribution from only the transition luminosity.The black line is the probability distribution that results by combining both constraints.bution,  (  ) , (, ), for each outburst  according to equation (4), we then find the overall mass/distance probability distribution by multiplying them and normalizing again: where  is the total number of outbursts and K is another normalization factor.To produce the overall estimate of only the soft state or the soft-to-hard transition mass-distance joint probability distribution, we use equation ( 6) but with only the soft state or the soft-to-hard transition probabilities for each outburst (see for example the joint probability distribution for GX 334-9 in Figs.B1, and B2).

Results
Using the results from our spectral modeling and the method described above, we successfully estimated the distance to 26 sources.Two of the sources, GX 339-4 and IGR J17091-3624, had multiple observed outbursts, so we combined the individual outburst results to an overall distance/mass estimate for these sources as per the method specified in Section 5.3.1.18 of our sources had constraints from both, the soft and the soft-to-hard transition.The distance marginal probability distributions and the joint distance/mass distributions for these 18 sources are presented in Figs. 2, B1, and B2 (note that some outbursts was observed by multiple instruments).From these sources, the majority show a large overlap in the distance marginal probability distribution between the constraints from the two states (red line: soft state, blue line: transition point), resulting in somewhat broad combined probability distributions (black line).We note that in most cases, the constraint from the transition flux suggests larger distances than the constraint from the soft state disk normalization.flux constraints.In particular, for IGR J17091-3624, the soft state constraint suggests a distance larger than 20 kpc, while the transition flux constraints from the same outbursts require distances in the 12 to 17 kpc range.The combined distance distribution has a narrow peak around 14 kpc, following the better-constrained transition fluxes.For the case of H 1743-322, the discrepancy is much smaller, with the soft state peak slightly higher than the transition flux peak.
In addition to analyzing the eighteen sources for which we have constraints from both states, we were also able to find the proba-Table 1. Summary of distance and mass estimates when using both the soft state and the soft-to-hard transition spectral information.For each source, we show the median distance/mass along with their 1 statistical uncertainties (see Section 6.4 for a discussion of systematic errors).The probability distributions are shown in Figs.B1, and B2 -Notes. "Combined" means the resulting distance or mass estimate is obtained though the combination of probability distributions from all the outbursts using equation 6.
Quoted errors are statistical-only.
bility distributions for eight additional sources, using only the soft constraint.In Figure B3, we plot the probability distribution of these eight sources.
Six of these sources (Swift J1727.8-1613,MAXI J1727-203, XTE J1720-318, XTE J1752-223, XTE J1817-330, XTE J1818-245, and XTE J1859+226) appear to have their distance probability distributions peaks in the lower distance range between 1 and 5 kpc.While the rest show higher distance distributions.We find that the spread of the soft-only distributions is in general wider than the distributions found by combining the probability distributions from both the soft and transition constraints, as one might expect.For example, the narrowest distribution of the source XTE J1752-223 as observed by MAXI has a 1 interval of (1.02, 2.66) kpc which is considerably larger compared to the narrowest distribution we get for distance using both constraints for the source MAXI J1348-630 also as observed by MAXI which has a 1 interval of (2.26, 2.62) kpc.As was the case for the two constraint sources, it is evident that the soft state constraint typically necessitates shorter distances.
Additionally, a few of the sources from the 18 that we found both constraints for, had observations from other outbursts/instruments that we were only able to obtain a distance estimate using the soft state.So, we added their results to Fig. B3.
A summary of the resulting distance constraints for each source is shown in Tables 1 and 2. Table 1 shows the results for sources where both the soft state and soft-to-hard transition were utilized, and Table 2 shows the estimates when using the soft state only.In these tables, we quote the median distance/mass and the 68% (1) credible interval (statistical errors).The median of the distances using both the soft state and soft-to-hard transition information is 7.47 +3.494.02 kpc, the closest source is MAXI J1348-630 with a distance 2.42 +0.2 −0.16 kpc and the farthest is IGR J17091-3624 with distance 14.16 +2.03 −1.94 kpc (combined outbursts).Furthermore, we find that the mass distribution for some of the sources was updated.Using the updated overall mass estimate, we report that the median of all masses is 8.46 +1.41 1.45  ⊙ , with lowest mass 7.36±0.86 ⊙ (IGR J17091-3624, combined outbursts) and highest mass 12 +0.77−0.68  ⊙ (XTE J1752-223).Figs. 3 and 4 show the overall distributions of the distances and masses of all sources in our sample, respectively.Both of these distributions are derived by summing all the individual source probabilities and then re-normalizing such that the total area under curves equals one.Note that for the overall distance distribution (Fig. 3), we have included the soft-only estimates for sources that did not have both constraints.We also avoid double-counting sources by taking the narrowest probability distribution (i.e.estimates with lowest errors) of the common sources across instruments.The resulting distance distribution reveals narrow peaks at 2.4 kpc and 7.1 kpc, and slightly  wider peaks at 4 kpc and 8.1 kpc.We also observe a heavy tail that rabidly decays beyond 15 kpc.The final mass distribution shows a shift in the observed peak mass value from the prior of 7.8 to 8.7  ⊙ .Additionally, there is an increase in the probability of finding masses in the 10 to 12.5  ⊙ range.Furthermore, using the tightest (lowest 1 errors) median distance estimate we obtained for each of the sources in our sample as well as their RA and Dec, we produced the Milky Way spatial distribution of the sources in our sample.We have used the astropy.coordinatesmethod to transform the (ICRS frame) RA, dec, and distance into galactocentric coordinates.We then plotted the possible positions of the sources in the projection of the Milky Way as shown in Fig. 5.The plot shows an abundance of sources in the region near or at  Bulge.We also see some distant sources observed by Swift/XRT, RXTE/PCA with the majority of sources likely to be associated with the spiral arms.

DISCUSSION
In this study, we have used data from MAXI/GSC, Swift/XRT and RXTE/PCA to develop a robust distance estimation method for black hole low-mass X-ray binaries.Utilizing our distance estimation method and the spectral information of the soft spectral state and the soft-to-hard spectral state transition, we have provided high-quality distance constraints for 26 unique BH-LMXB sources, often exceeding the accuracy of existing constraints (see Section 6.2).From our results, we constructed a new overall distribution of low-mass black hole X-ray binaries' distances from Earth.Moreover, we have obtained an updated overall black hole mass distribution by combining the distance/mass distribution from both the soft state and soft-tohard transition state.The python pipeline used to calculate the probability distributions can be found in a public GitHub repository11 .
We have also implemented a simplified web version of our improved method for the purpose of providing quick distance estimates for newly discovered BH-LMXBs12 .

Comparison to Previous Estimates
First, we want to assess the agreement of our newly determined distance estimates with existing distance constraints.To facilitate an easy comparison, we plot the estimated distance from each instrument against the most accurate literature distance in Fig. 6.We limited ourselves to those sources, where we were able to use both spectral states in our analysis.
From the graph, it is evident that distances derived from MAXI and RXTE data generally show good agreement with previous estimates.They are either consistent within errors or almost equal within the error bars.Most notably, our estimates for MAXI J1820+070 (red point) and MAXI J1348-630 (blue point) using MAXI/GSC data are in perfect agreement with the highly accurate Gaia parallax distance (Gaia Collaboration et al. 2021;Tetarenko et al. 2021) and H i absorption distance (Chauhan et al. 2021), respectively.
Distances derived from Swift data also show good agreement for three sources, which include MAXI J1820+070.On the other hand, two sources at distances ∼ 7 kpc appear to have considerably higher distance estimates from our method than the previous measurements.These two sources are XTE J1908+094 (2013 outburst, orange point) and GRS 1739-278 (2014 outburst, green point).Below, we take a closer look at these sources and their previous distance constraints.
XTE J1908+094 was first discovered in February 2002 (Woods et al. 2002).The source later went into two other outbursts in 2013 and more recently in 2019 (Krimm et al. 2013;Rodriguez et al. 2019).The first attempt to estimate the distance to this source was made soon after its discovery in 2002, where in 't Zand et al. (2002) used the outburst peak luminosity to derive a lower limit of 3 kpc.Later, in 2006, Chaty et al. ( 2006) conducted a detailed astrometric analysis and identified a likely candidate for the companion star, which they called "object 1".Using a  − ′ ,   ′ , absolute colour-magnitude diagram, they determined that "object 1" is an intermediate or late-type mainsequence star.With spectral type between A and K. Consequently, they determined that its distance should be in the range of 3-10 kpc.This estimate is what we compare to our result.Using our method, we measure a distance of 11.52 +1.93  −2.74 kpc.Although we estimate a median distance considerably higher than the average of their lower and upper limits, the two estimates are still in agreement within the error bars.Moreover, we point to the fact that according to Jonker & Nelemans (2004), using template non-rotating stars to determine the spectral type could lead to underestimation of distance.So if the distance was underestimated by Chaty et al. (2006), that makes our estimates consistent.It is also worth noting that using the same softto-hard transition luminosity empirical relation by Maccarone, T. J. ( 2003) that we used, Curran et al. (2015) found an almost identical range (7.8 -13.6 kpc) to what we found.However, our transition state probability distribution of XTE J1908+094, which is shown in Fig. B1, indicates that the possible distances using only the transition state flux are significantly higher than they found.This is evident from the broad peak of the distribution that appears to be in the range of 11-24 kpc.This discrepancy with the estimation of Curran et al. (2015) is most likely due to our different determination of the state transition date, which led them to use a higher transition flux than we used.In their opinion, the state transition was not observed by Swift/XRT due to its proximity to the sun at that period and that in the earliest observation after this pause, the source had already returned to its hard state.In contrast, our HR analysis determined that the source was still transitioning into the hard spectral state after this pause.In their paper, the authors state that the source spectrum was still hardening until MJD 56730 (our transition state spectrum date is MJD 56703), which is consistent with our determination.It is unclear why the authors considered the source already in the hard state and estimated a higher transition flux and closer distances.
On the other hand, GRS 1739-278 was discovered initially in its 1996 outburst (Paul et al. 1996;Vargas et al. 1997).The source went into another outburst in 2014 and was observed by Swift (Krimm et al. 2014).The main 2014 outburst lasted over a year, and the softto-hard transition was not observed due to the Sun's obscuration.Following this main outburst, two mini-outbursts were also observed (Yan & Yu 2017).In the second mini-outburst, observations during the soft-to-hard transition were obtained (Yan & Yu 2017).As the soft-to-hard transition was observed during that mini outburst, we have used it to estimate a distance of 12.74 +3.39  −2.42 kpc.In contrast, using interstellar extinction Greiner et al. (1996) found the distance to the source in the range of 6-8.5 kpc.This estimate was later revisited by Yan & Yu (2017), where they used an  H from Swift/XRT spectral modelling as well as an updated extinction map to estimate a distance of 7.5±2 kpc.While their upper limit of 9.5 kpc is very close to our lower bound of 10.32 kpc (1), it is still slightly higher.Once again, according to Jonker & Nelemans (2004) systematic bias in estimating  H through X-ray observations leads to underestimation of distance.This could possibly explain this discrepancy.However, by looking at the distance marginalized distribution in Fig. B1 as well as the softonly estimates in Table 2 and Fig. 7, we notice that the higher distance values in our derived combined distance distribution are mainly due to the transition constraint.Due to this fact, we think that another possible explanation for the discrepancy between our estimates is that this source might have its transition luminosity occur at a value lower than the 1-4% range we used.If we use a lower value of the Eddington fraction, then the transition constraint would shift toward lower distances at the same masses and, in turn, make the combined distribution peak at lower distances.This is entirely possible since we have used a mini-outburst transition flux.The mini-outburst sources observed up till now are a handful, and no systematic studies were made on whether they also follow the empirical correlations found in major outbursts (Yan & Yu 2017).However, it is worth noting that preliminary analysis shows that the same mechanisms drive them, and thus, we expect them to behave similarly (Yan & Yu 2017).
While it is clear that we can obtain better constraints from our method if we information from both states, we also want assess the robustness of our distances derived from the soft state spectra alone.In Fig. 7, the same comparison as above is made, but for the outbursts/sources where we could only find the soft constraint.In the graph, data points from MAXI/GSC and RXTE/PCA suggest that the distance is systematically underestimated by using the soft-only constraint.On the other hand, estimates made using data from Swift/XRT appear to be less affected by the loss of information from the transition state and agree very well with previous estimates from the literature.Nevertheless, when looking at the same source (GX 339-4, black point) that was observed by all instruments, we find consistent evidence that the soft state distance constraints are systematically lower than literature distances.To investigate this further, we also plotted a comparison between our estimates that were found using both con-straints and those that were found using the soft-only constraint in Fig. 8.The plot shows that apart from six data points at distances larger than 5 kpc, all other data points indicate that the soft state constraint estimates are systematically lower than estimates found when using both constraints.This is also corroborated by the distance distributions of each source as shown in Fig B1 (compare red lines to blue lines).Inspection of the differences between the two estimates (our soft-only estimates versus estimates using both constraints) reveals that distances obtained through the soft-only constraint might be un- Comparison of the distance constraints obtained using information from both the soft and transition spectra and those obtained using the soft spectra only.A trend is apparent where the soft-only constraints yield lower distances than those constraints using the spectral information from both states.
derestimated by up to 70%.However, exceptions exist for sources located further than 5 kpc, where the soft-state constraint suggests larger distances are needed.This discrepancy could stem from a systematic bias associated with the disk spectrum model we employed.
We delve into potential biases with greater detail in Section 6.4.

Consistency Between Instruments
To investigate whether there are any biases arising from the use of any particular detector, we examined the distance estimates determined from outbursts that observed in multiple instruments.We looked at estimates obtained by using both constraints and through using the soft-only constraint, as the number of overlapping outbursts were small.We found that all estimates are in agreement within their uncertainties.Our findings indicate that there is no detector-specific bias when applying our method to data from MAXI/GSC, Swift/XRT, or RXTE/PCA.This is in line with the expected results since the energy bands used for modelling the data from each instrument are nearly identical.Thus, spectral modelling results are not expected to change significantly.

Fractional Uncertainty and Dependence on Flux
So far, we have considered the absolute errors of our distance estimates when determining the quality of the constraints.Using the absolute error, however, ignores that large distances might carry a larger uncertainty as the source flux from these more distant sources would be lower.To be able to compare the quality of measurements better, we first find the average fractional uncertainties of each measurement by calculating the average 1 upper/lower errors, and then we divide this average error by the median distance estimate.Once the fractional uncertainties have been found, we investigate the relationship of the fractional uncertainties with flux.Given that our dataset is small, we calculated the Kendall's tau correlation coefficient between the fractional uncertainty and the log 10 of each of the soft and the transition fluxes.Both Kendall's tau correlation coefficients were close to zero and had -values of 1.0 and 0.86 for the soft and transition fluxes, respectively.This suggests little to no evidence of an association between the fractional uncertainty and the fluxes.We also tested these correlations for each instrument separately and found similar -values.Having determined that the fractional uncertainties have no detectable dependence on the observed flux.Using all estimates we find a median fractional uncertainty of 0.2, when using spectral information from both the soft and transition states.
Repeating the above tests for the fractional uncertainties obtained through the soft-only constraint, we calculated the Kendall's tau correlation coefficients separately for each instrument.The Kendall's tau correlation coefficients for Swift and RXTE indicated there is no association with the soft state log flux.Since the coefficients were close to zero and their -values were equal to 0.84 and 0.28, respectively.On the other hand, the test provided weak evidence of a relationship between MAXI's fractional uncertainties and the soft state flux with a coefficient value of -0.45 (-value = 0.05).So, we conclude that the uncertainties in the distance estimates from MAXI decrease, albeit weakly, as the source's log 10 flux increases when only the soft state constraint is found.While for Swift and RXTE, there is no detectable dependence on the source's flux.We then calculate that the median fractional uncertainty is 0.48 when considering all estimates found using the soft-only constraint.Although more data are needed to provide more robust conclusions regarding the relationship between fractional uncertainty and flux, our analysis here provides a preliminary investigation.

Improved Distances
Our method significantly improved the quality of the distance constraint for several sources.In this section, we take a closer look at sources where we found much more precise distance constraints, for sources which previously only had upper or lower limits on their distances.
We start by looking at Swift J1910.2-0546, which is a peculiar source because it exhibited two state transitions13 (Nakahira et al. 2014;Reis et al. 2013).Applying our method to data from the 2012 outburst of Swift J1910.2-0546obtained by MAXI/GSC, we were able to find a precise distance estimate to the source of 3.8 +0.23  −0.25 kpc.The only available estimate for this source in the literature was made by Nakahira et al. (2014) using the simple approach of our method, which assumes a non-rotating black hole and enforces discrete estimates on inclination.The authors only reported a lower bound of 1.7 kpc.Since the authors used the same approach as our method, the consistency between our measurements is not surprising.Another source is XTE J1652-453; the distance had been previously estimated to be ≥ 5 kpc (Reynolds et al. 2009), using the spectrum of the companion star.Applying our method to data from RXTE/PCA provides the most precise constraint on the distance of XTE J1652-453, as we calculated a distance of 6.69 +0.54  −0.26 kpc.Moreover, using our method on data from RXTE/PCA of the 1998 outburst of XTE J1748-288 and the 2008 outbursts of Swift J1842.5-1124,we confined their distances to 8.66 +1.98  −1.84 kpc and 8.12 +1.67 −2.2 kpc, respectively.The prior estimate of XTE J1748-288 was ≥8 kpc (Hjellming et al. 1998), which was obtained using an H i absorption measurement.While the best estimate of Swift J1842.5-1124 of > 5 kpc (Zhang et al. 2022) was found by utilizing radio and X-ray correlations.
Additionally, using the MAXI/GSC data of Swift J1727.8-1613 in its ongoing 2023 outburst, we have utilized the soft state spectra to estimate a distance of 1.52 +0.85  −0.61 kpc.Although the source is currently still in outburst, the HR data suggest that the source is in the soft state (Bollemeĳer et al. 2023).This estimate perfectly agrees with the current only estimate in the literature of 1.5 kpc using flux scaling arguments (Veledina et al. 2023).We point that our estimate is obtained through the soft state information only, so the distance might be underestimated by up to ∼ 70% (see Section 6.1.1).
In general, we have improved the distance estimates for 10 sources, achieving up to a 45% decrease in fractional uncertainty for Swift J1910.2-0546compared to previous estimates (see Table C1).
However, it is important to note that, while our pragmatic framework has found much tighter constraints than many previously estimated distances, we must emphasize that our errors are statistical and may be subject to systematic biases.In our discussion in Section 6.4, we will attempt to identify these biases and suggest ways to refine our framework to accommodate them in future work or in form of adjustments to the pipeline.

Distance and Spatial Distributions
In Fig. 3, we plotted the overall distance distribution obtained by adding the probability distributions from each source in our sample, visualizing the probability of finding BH-LMXBs as a function of distance.It is tempting to correlate the multi-peak structure of the observed distribution with the Milky Way's spiral structure.However, since these distances are measured directly from the sun to the respective sources, they do not necessarily provide insights into the Galaxy's overall structure.Thus, we cannot infer any relationship between Milky Way's structure and this distance distribution.Furthermore, this distribution was obtained considering only the statistical uncertainties.Including systematic errors (see Section 6.4) would possibly smear the spiky shape reminiscent of the spiral arms over-densities.
On the other hand, the spatial distribution of sources in our sample, as shown in Fig. 5, shows that the sources are likely to be located near spiral arms or the Bulge.Comparing this observed distribution to a similar observed distribution as in Corral-Santana et al. (2016); Gandhi et al. (2019), we found that our inferred distribution is consistent with their distribution.In particular, all distributions find that BH-LMXBs are concentrated towards the direction of the Galactic Bulge, with distances indicating positions at or very near the Galactic center.We also find a small number of sources located further away from the Galactic center.Careful investigations on the selection effects affecting our observed distributions will be the subject of one of our future paper.

Mass
On the other hand, looking at the posterior mass distribution in Fig. 4. We find that the probability density peak appears to have moved slightly higher (to 8.7  ⊙ from 7.8  ⊙ ) when compared to our prior guess.Thus, this indicates that we have gained new information on the most frequent or likely mass values, which are now slightly higher than before.Additionally, there's a notable increase in probability density within the 10 to 12.5  ⊙ range.This observed mass distribution appears to reaffirm the observed lower mass gap between the 3 and 5  ⊙ (Özel et al. 2010;Farr et al. 2011;Kreidberg et al. 2012).This is expected since we use the distribution found by Özel et al. (2010), which is our prior mass guess.But since there was no a priori constraint on the direction of the information gain in the posterior mass distribution (i.e.we could have found sources where both constraints only agree at much lower masses), we argue that our results contribute to the evidence of the observed lower mass gap.Nonetheless, we cannot exclude that it is due to an observational bias.Further detection of new BH-LMXBs and determining their masses would pave the way to understanding these possible biases.

Limitations and Biases
In the previous sections, we have demonstrated that our framework has the potential to provide robust estimates for most BH-LMXBs.However, we must discuss a number of limitations and possible biases that may arise from the assumptions we made.
First, in our disk model, we have assumed a color correction factor  col = 1.7.This correction factor is needed because the emission from the disk is not a perfect blackbody radiation (Ebisuzaki et al. 1984).Shimura & Takahara (1995) argued that the spectra produced could be approximated by a blackbody with maximum temperature shifted to higher temperatures by  col .By numerically solving the equations of the vertical structure of the disk, Shimura & Takahara (1995) found this factor should be in the range 1.5-1.9, with 1.7 being the most acceptable approximation.Since  col enters in the distance/mass relation found by modelling the soft state (equation 3), it affects our distance estimations.Several studies (Merloni et al. 2000;Davis et al. 2005) have shown that there is still large uncertainty in the ranges that  col can take.Nonetheless, most theoretical models suggest it should be in the narrow range of 1.4-2.0(Davis et al. 2005;Davis & El-Abd 2019).Our assumption of  col = 1.7 was motivated by empirical evidence, which found that most systems are consistent with the value of 1.7 (Reynolds & Miller 2013).We also note that for the disk-dominated state (soft spectral state) with disk temperatures of < 1 keV (typical of all our spectral fits), the latest theoretical models predicated that  col should be 1.7-1.8(Davis et al. 2005;Reynolds & Miller 2013;Davis & El-Abd 2019).Furthermore, the higher values  col are only predicted for systems, which are accreting at high Eddington fractions (> 30%) (Davis & El-Abd 2019).Observational evidence hints that these systems are rare (Davis & El-Abd 2019).Additionally, relativistic disk models14 already utilize the Novikov-Thorne disk model (Novikov & Thorne 1973), which assumes thin disks with Eddington ratios less than 30%.Therefore, considering higher values seems unwarranted as it is not consistent with the models we use.Notwithstanding, we identify this as a possible bias source.In preliminary analysis, where we allowed  col to vary in the range of 1.4 to 2.0 for MAXI J1820+070 using Swift/XRT data, we have found minimal change in the distance from 2.64 +0.71 −0.56 kpc to 2.63 +0.71 −0.55 kpc.This very minimal change is almost certainly associated with the tight soft-to-hard transition constraint of this source (see Fig. B2) which means that the overall estimate is not highly sensitive to changes in  col .Thus, being able to draw conclusions on this matter would require analysis of multiple sources with varying degrees of narrowness for the soft-to-hard transition constraint.Since including the uncertainty of  col in addition to other certainties in our current implementation of the method renders the probability calculation very computationally intensive, we have reserved full investigation of the effect of  col on our estimations for future work.
Another source of bias is the use of the Özel et al. (2010) mass prior.We argue that our motivation was purely observational since several studies have reinforced this prior (Farr et al. 2011;Kreidberg et al. 2012).But we also note that there are some possible physical mechanisms that can explain the observed distribution, especially for the low mass gap at 2-5 ⊙ , such as natal kicks that favour the formation of X-ray binaries in a certain mass range (Fryer & Kalogera 2001;Fryer et al. 2012).We also note that using both the soft and transition constraints decreases the effect of this bias.Because agreement between both constraints is required for updating the prior, as shown in Fig. 4.
The use of the 1-4% Edd assumption of the soft-to-hard transition luminosity (Maccarone, T. J. 2003) is also another source of bias.It has been suggested that the state transition luminosity might depend only on the dimensionless viscosity parameter, and the observational evidence suggests that this parameter does not appear to change significantly between different sources (Maccarone, T. J. 2003;Vahdat Motlagh et al. 2019).As we mentioned in the methods section, several studies found observational evidence for the same range proposed by Maccarone, T. J. (2003).Most recently, a study by Wang et al. (2023) investigated the transition luminosities again using data from Swift/XRT and NICER/XTI.The authors found that the power law flux is tightly constrained around 1.3%.Moreover, the authors point to the importance of the methodology for determining the state transition.They argue that the use of disk recession (when disk normalization begin to increase) would be simpler to verify when analyzing spectra that have been observed in the soft bands (0.5-10 keV).Given that we used the common approach of determining the state transition through HR analysis, this is yet another source of bias because the definitions of HR can vary between studies.Nevertheless, we believe that the 1-4% range we considered is large enough to include the uncertainties in determining the state transition dates.
A limitation of our method is the presence of insufficient spectral data at the state transition.We note that we could not obtain a distance estimate for most of the sources from our initial sample because they suffered the absence of viable observations at the state transition.Additionally, our most robust measurements are obtained using both constraints.Fortunately, in the future, continuous monitoring from missions such as MAXI/GSC and Swift/BAT would reduce cases of fewer observations in the transition state.
Finally, a systematic bias in our estimates could be from the assumption of prograde spins only ( > 0).Strong observational evidence shows that the majority of discovered BH-LMXBs have prograde spins (Reynolds 2021).However, we have to note that there are some exceptions where it has been suggested that a retrograde spin is most likely (e.g.Reis et al. 2013;Rout et al. 2020).This would affect the soft state constraint and would increase the probability density in the direction of higher distances15 .To further investigate what impact the inclusion of retrograde spin would have on the calculated distances, we recomputed the distances for six sources, employing a flat prior for  ranging from -0.999 to -0.998.The sources analyzed include GRS 1739-278, with a recalculated distance of 14.05 +3.2   −3.4   kpc; MAXI J1348-630, at 2.52  1 reveals that, except for Swift J1910.2-0546, the median distances for all sources have increased.Additionally, the fractional uncertainty in most of the distances has increased.This is a consequence of allowing for retrograde spin, which widens the soft state constraint and shifts its peak to larger distances.This alteration further results in a greater overlap between the soft and transition distributions at larger distances, leading to an increase in the 1 statistical errors of the combined distribution.
As a final note, we point to the fact that our aim was to use the best current theoretical and observational knowledge to obtain reliable estimates of the distances to BH-LMXBs.Thus, we have made a number of pragmatic assumptions to be able to achieve this goal.We also emphasise that our versatile framework could be easily modified to address some of these possible biases if more observational evidence surfaces against the assumptions we used in our work.So, the framework's robustness lies not only in the current implementation but also in its ability to evolve with the latest observational evidence and theoretical understanding.For instance, further modification to our method would include the possibility of accounting for  col uncertainties, retrograde spin, or a spin prior based on the observed BH-LMXB population.

CONCLUSION
We have analyzed the light curves of 50 transient low-mass X-ray binaries selected from the catalogues of Tetarenko et al. (2016), Corral-Santana et al. (2016), and Alabarta et al. (2021) using archival data from MAXI/GSC, Swift/XRT and/or RXTE/PCA.We further selected sources that had an identifiable soft state and excluded sources with very poor spectral fits ( red > 1.4).Consequently, this led to a final sample of 26 sources.We performed X-ray spectral analysis of 39 soft-state spectra and 27 soft-to-hard transition spectra for a total number of 34 unique BH-LMXB outbursts.The main results from our study can be summarized as follows: (i) We have developed a robust framework to constrain the black hole BH-LMXB distances from their X-ray spectra, considering general relativistic effects and other physical parameter uncertainties.
(ii) Using this framework, we constrained the distances to 26 sources, often improving their existing constraints.
(iii) Comparison to previous distance estimates shows good agreement.In particular, our estimates are in good agreement with the highly accurate estimates obtained through Gaia parallax measurements and H i absorption methods.
(iv) We have not detected any instrument-specific biases in the distance estimates obtained using data from MAXI/GSC, Swift/XRT, or RXTE/PCA.
(v) We have found that estimates obtained using both the soft and transition spectral information have a median uncertainty of 20% while estimates obtained using only the soft state spectrum have a median uncertainty of about 50%.
(vi) We have created and published a python pipeline that can be used to calculate the distance/mass probability densities and has implemented an online tool based on the results of this study that can be readily used to quickly calculate the distance estimates.
(vii) Using our distance estimates, we found an overall probability distribution of distances (Fig. 3) from the sun with peaks at 2.4, 4, 7.05, and 8.1 kpc.Investigating the spatial distribution of the sources in the Milky Way plane projection (Fig. 5), we found that most of the sources are directed toward the Bulge and near the spiral arms.
(viii) The posterior mass distribution found by requiring the estimates from the soft state and the transition state must agree on the true value of distance, contributes to the evidence for the 3-5 ⊙ mass gap (e.g.Özel et al. 2010;Fryer et al. 2012).
We note that our distance estimation framework can be applied on data from other X-ray detectors such as NICER.Nevertheless, due to the low amount of outbursts observed with NICER compared to other the instruments we used, we have reserved analysis using NICER for future projects when more data is available.
ac.uk/user_objects/.Observational data from the Monitor of All-sky X-ray Image/Gas Slit Camera (MAXI/GSC) were also employed in this study.The MAXI/GSC data are freely accessible via the MAXI team's website hosted by the RIKEN, JAXA, and the MAXI team http://maxi.riken.jp/mxondem/and part of the HEAsoft package https://heasarc.gsfc.nasa.gov/docs/software/heasoft/.Additionally, data from the Rossi X-ray Timing Explorer/Proportional Counter Array (RXTE/PCA) were used.These data are available in the public domain and can be retrieved from the NASA's HEASARC database: https://heasarc.gsfc.nasa.gov/cgi-bin/W3Browse/w3browse.pl for the spectra and https://www.isdc.unige.ch/integral/heavensfor light curves.Furthermore, the pipeline used to generate the distance/mass probability densities is available on GitHub at https: //github.com/ysabdulghani/lmxbd,and the scripts for data reduction and preprocessing can be provided upon request.
Table A1.This table contains the information of each spectrum modelled in this work.We show the source name, outburst year, and instrument.We also show the selected spectra dates for the soft state spectrum and, if available, the soft-to-hard transition spectrum (See Section 3.2 for more details on selection method).We list the  H and its reference, if fixed.The  2 red for the fitted soft state spectrum and, if applicable, the soft-to-hard transition spectrum are also listed.Finally we specify the model used in the fitting (see Section 4 for the definition of models).The sources are ordered in alphabetical order from top to bottom.

Source
Outburst  This figure shows the marginalized distance probability distribution for all of the outbursts in our sample that both constraints.The red line shows the probability distribution from only the soft state constraint, while the blue line shows the probability distribution from only the transition luminosity.The black line shows the probability distribution that results by combining both constraints such that both are assumed to be true, equation (( 4)).For sources with multiple outbursts, we only use the combined probability distribution.* See Table 1 for outbursts dates used to obtain the combined estimate.
MNRAS 000, 1-19 (2024) Estimating Distances to BH-LMXBs 21 This figure shows the joint mass-distance probability distribution for all of the outbursts in our sample using both constraints.The red region shows probability distribution from only the soft state constraint while the blue region show the probability distribution from only the transition luminosity.The gray region shows the combined probability distribution that results when both distributions are assumed to be true, equation (4).All contours lines are 1, 2, and 3 regions.For sources with multiple outbursts we only use the combined estimated probability distribution.* See Table 1 for outbursts dates used to obtain the combined estimate. Quoted errors are statistical (1).Refer to Section 6.4 for discussion of systematic errors.

Figure 1 .
Figure 1.This diagram summarizes the statistical framework, which produces mass/distance probability distributions.The framework is discussed in detail in section 5.3.

Figure 2 .
Figure2.This figure shows the probability distributions resulting from applying our method to the RXTE/PCA observations of Swift J1539.2-6227 during its 2008/2009 outburst.Left: The joint mass/distance probability distribution.The red region corresponds to the probability distribution from only the soft state constraint, while the blue region shows the probability distribution from only the transition luminosity.The grey region visualizes the combined probability distribution that results when both distributions are assumed to be true, equation (4).The contour levels shown are 1, 2, and 3.Right: The marginalized distance probability distribution.The red line corresponds to the probability distribution from only the soft state constraint, while the blue line indicates the probability distribution from only the transition luminosity.The black line is the probability distribution that results by combining both constraints.

Figure 3 .
Figure3.This graph shows the overall heliocentric distance probability distribution, obtained by summing the individual distributions, for the 26 sources in our sample.The plot shows peaks at 2.4, 4.0, 7.1, and 8.1 kpc.We avoid double-counting sources by taking the narrowest probability distribution for the sources observed by multiple instruments.
"Combined" means the resulting distance estimate is obtained though the combination of probability distributions from all the outbursts. Quoted errors are statistical.

Figure 4 .Figure 5 .
Figure 4.This plot shows the observed overall mass probability distribution (black) and the prior mass distribution (red line).The observed overall distribution shows that the peak shifted from 7.8 to 8.67  ⊙ when compared to the peak of the prior.We also see a higher probability density in the range of 10-12.5  ⊙ compared to the prior Gaussian distribution N(7.8,1.2).

Figure 6 .Figure 7 .
Figure 6.This plot compares the distance constraints obtained from MAXI/GSC, Swift/XRT, and RXTE/PCA data (this work) and distance estimates found in the literature.The dashed line is  = .The error bars in the y-axis correspond to 1 while the error bars in the x-axis correspond to the quoted errors in each of the respective literature works used, where arrows indicate either lower or upper limits.References for the literature distance can be found in Table C1.The colors correspond to the method(s) used to estimate the literature distance.Methods.A: Gaia parallax, B: H i absorption, C: Reflection & continuum-fitting, D: Interstellar extinction, E: Companion star, F: X-ray decay time scale, G: X-ray luminosity & disk temperature, H: Hard-to-soft transition luminosity, I: Proper Motion (Jet), J: Soft state and soft-to-hard transition without corrections, K: Radio & X-ray correlation, L: Assumed to be at Galactic Bulge.M: Flux Scaling, O: Soft state & peak luminosity.
Figure8.Comparison of the distance constraints obtained using information from both the soft and transition spectra and those obtained using the soft spectra only.A trend is apparent where the soft-only constraints yield lower distances than those constraints using the spectral information from both states.
Figure B1.This figure shows the marginalized distance probability distribution for all of the outbursts in our sample that both constraints.The red line shows the probability distribution from only the soft state constraint, while the blue line shows the probability distribution from only the transition luminosity.The black line shows the probability distribution that results by combining both constraints such that both are assumed to be true, equation ((4)).For sources with multiple outbursts, we only use the combined probability distribution.* See Table1for outbursts dates used to obtain the combined estimate.
Figure B2.This figure shows the joint mass-distance probability distribution for all of the outbursts in our sample using both constraints.The red region shows probability distribution from only the soft state constraint while the blue region show the probability distribution from only the transition luminosity.The gray region shows the combined probability distribution that results when both distributions are assumed to be true, equation (4).All contours lines are 1, 2, and 3 regions.For sources with multiple outbursts we only use the combined estimated probability distribution.* See Table1for outbursts dates used to obtain the combined estimate.

Figure B3 .
Figure B3.This figure shows the distance probability distribution for sources/outbursts from our sample with only spectral information from the soft state.
parallax, B: H i absorption, C: Reflection & continuum-fitting, D: Interstellar extinction, E: Companion star, F: X-ray decay time scale, G: X-ray luminosity & disk temperature, H: Hard-to-soft transition luminosity, I: Proper Motion (Jet), J: Soft state and soft-to-hard transition without corrections, K: Radio & X-ray correlation, L: Assumed to be at Galactic Bulge.M: Flux Scaling, O: Soft state & peak luminosity. .

Table 2 .
Summary of distance estimates using the soft state spectral information only for all sources used in this work.For each source, we show the median distance along with 1 errors statistical uncertainties (see Section 6.4 for a discussion of systematic errors).The probability distribution of distance are shown in Fig.B3.

Table C1 .
In this table we show the best estimated distance (this work) for each source in our sample as well as the most accurate distance from the literature.We also show the distance method used to obtain the literature estimate.Sources are ordered alphabetically from top to bottom.