High-density reflection spectroscopy of black hole X-ray binaries in the hard state

We present a high-density relativistic reflection analysis of 21 spectra of six black hole X-ray binaries in the hard state with data from \textit{NuSTAR} and \textit{Swift}. We find that 76\% of the observations in our sample require a disk density higher than the 10$^{15}$~cm$^{-3}$ assumed in the previous reflection analysis. Compared with the measurements from active galactic nuclei, stellar mass black holes have higher disk densities. Our fits indicate that the inner disk radius is close to the innermost stable circular orbit in the hard state. The coronal temperatures are significantly lower than the prediction of a purely thermal plasma, which can be explained with a hybrid plasma model. If the disk density is fixed at 10$^{15}$~cm$^{-3}$, the disk ionization parameter would be overestimated while the inner disk radius is unaffected.


INTRODUCTION
Accreting black holes are efficient at converting gravitational energy into electromagnetic radiation (e.g. Shakura & Sunyaev 1973;Thorne 1974). The accretion process is believed to be responsible for many high energy phenomena in the Universe, i.e., active galactic nuclei (AGNs) and X-ray binaries (XRBs). Moreover, accretion is an important ingredient in understanding the growth of supermassive black holes (SMBHs) and galaxies (e.g. King 2003;Alexander & Hickox 2012;Fabian 2012). The X-ray emission is a powerful probe to study the innermost regions of accreting black holes (e.g. Mushotzky et al. 1993;McHardy et al. 2006).
One of the primary components in the X-ray spectra of black holes is the power-law continuum with a high-energy cutoff. This component is thought to originate from inverse Compton scattering of seed photons by a hot plasma (the corona) near the black hole (e.g. Shapiro et al. 1976;Haardt & Maraschi 1993). This coronal emission can illuminate the optically thick accretion disk and produce the reflected emission that is featured with fluorescent emission, photoelectric absorption and the Compton scattering hump (e.g. George & Fabian 1991;García & Kallman 2010). Due to strong relativistic effects near the black hole, the observed reflected features are significantly skewed (e.g. Fabian et al. 1989;Dauser et al. 2010;Bambi 2017a) and contain rich information about the spacetime (e.g. Reynolds 2014; Bambi 2017b; Bambi et al. 2021;Reynolds 2021) and disk-corona geometry (e.g. Fabian et al. 2009Fabian et al. , 2012. This technique has been successfully applied to many black hole XRBs and AGNs (e.g. García et al. 2015;Walton et al. 2016;Jiang et al. 2018;Liu et al. 2019) Previous reflection models assume the disk electron density to be log(n e /cm −3 ) = 15. This value is appropriate for very massive (e.g. M BH > 10 8 M ) black holes in AGNs (e.g. Jiang et al. 2019b) but the standard disk model (Shakura & Sunyaev 1973) predicts a higher density for less massive black holes. For AGNs, reflection-based analysis has shown that a disk density larger than log(n e /cm −3 ) = 15 is sometimes required for SMBHs with M BH < 10 7 M and the measured densities are consistent with the prediction of a radiation pressure-dominated disk (e.g. Jiang et al. 2019c;Mallick et al. 2022). The high-density reflection model is able to explain the "soft excess" feature found in many Seyfert galaxies (e.g. Arnaud et al. 1985) without including an additional component (e.g. García et al. 2019). The soft X-ray reverberation signature also supports the reflection origin of the soft excess (e.g. De Marco et al. 2013;Kara et al. 2013).
Moreover, a higher disk density can also relieve the puzzle of super-solar iron abundance in some reflection spectral modelling (e.g. Tomsick et al. 2018;Jiang et al. 2019a). Since the standard disk model predicts an anticorrelation between the disk density and black hole mass times the square of the Eddington-scaled accretion rate (M BHṁ 2 ) for a radiation pressure-dominated disk (e.g. Svensson & Zdziarski 1994, and Eq. 1), we expect an even higher density for the accretion disk of XRBs (especially in the hard state). This has been confirmed with an analysis of individual sources but the measured disk densities are lower than theoretical predictions of a radiation pressure or gas pressure-dominated disk (see Fig. 8 of Jiang et al. 2019c).
Previous high-density reflection modelling of black hole XRBs has only been conducted for a few sources (e.g. Reis et al. 2008Reis et al. , 2009Reis et al. , 2011Steiner et al. 2011;Reis et al. 2012b;Steiner et al. 2012b;Walton et al. 2012;Chiang et al. 2012;Tomsick et al. 2018;Jiang et al. 2019aJiang et al. , 2020Connors et al. 2021;Chakraborty et al. 2021). Moreover, on the log(n e )-log(M BHṁ 2 ) diagram, there is still an obvious gap between the XRBs and AGNs sample (see Fig. 8 of Jiang et al. 2019c). In this work, we conduct a systematic analysis on the broadband X-ray spectra of six black hole XRBs in the hard state that have not been studied with high-density reflection models. The data selection and reduction are described in Sec. 2. Sec. 3 explains the spectral analysis. In Sec. 4 we discuss the results.

OBSERVATIONS AND DATA REDUCTION
We select six black hole XRBs from the BlackCAT catalog 11 (Corral-Santana et al. 2016). All the sources have measurements of their distances, enabling the estimation of their bolometric luminosities through broadband spectral fitting. We also require the sources to have been observed by NuSTAR so the Compton hump is captured. Tab. 1 shows information about the selected sources. For those sources that do not have measurements of their black hole mass, we assume a value of 10±1 M throughout this work to calculate the Eddington-scaled luminosity. For each source, we select the NuSTAR observations in the hard state that show relativistic reflection features. One of the main features of the reflection by high-density accretion disk is the enhanced soft X-ray emission (Jiang et al. 2019c). Therefore, we also include contemporaneous (on the same day) Swift data if possible. Details of the selected observations are listed in Tab. 2. Note that a few other sources may also match our selection criteria but they have already been studied with reflection models of variable electron density (e.g. GX 339-4, Jiang et al. 2019a;GRS 1716-249, Jiang et al. 20204U 1630-47, King et al. 2014, Connors et al. 2021MAXI J1348-630, Chakraborty et al. 2021). We do not include those sources in the analysis of this work. MAXI J1820+070 also matches the selection criteria, but there is still debate on its accretion geometry (e.g. Buisson et al. 2019;Zdziarski et al. 2022). We will explore this source in a future publication.

NuSTAR
We produce cleaned event files for both FPMA and FPMB with the tool nupipeline v0.4.9 and the calibration version 20220301. The source spectra are then extracted from circular regions centered on the sources using the nuproducts task. The background spectra are extracted from source-free areas with polygon regions created with ds9. Note that for V404 Cyg, five fluxresolved spectra are extracted from the two observations (see details in Sec. B.4). The Swift/XRT data are all in the Window Timing (WT) mode and are free from pile-up effects. The cleaned events files are produced using xrtpipeline v0.13.7 and the last calibration files as of September 2021. We extract the source spectra from circular regions centered on the source with a radius of 100 arcsec. The background regions are chosen to be annuli with an inner radius of 110 arcsec and an outer radius of 200 arcsec. We only include events with grade 0.

SPECTRAL ANALYSIS
Spectral fittings are conducted with XSPEC v12.12.1 (Arnaud 1996). The NuSTAR data are used in the 3-79 keV band and the Swift data in the 1-10 keV band. All data are grouped to ensure a minimum of 30 counts per bin. We implement element abundances of Wilms et al. (2000) and cross-sections of Verner et al. (1996). χ 2 statistics is used to find the best-fit values and uncertainties (at 90% confidence level unless otherwise specified) of the model parameters.
We first fit the broadband spectra from each source with a simple absorbed continuum model: constant * tbabs * ( diskbb + nthcomp ). The constant model is to fit the cross-normalization between instruments. The absorption by Galactic inter-stellar medium is modelled by tbabs (Wilms et al. 2000). The diskbb (Mitsuda et al. 1984) component is to fit the thermal emission from the multicolor accretion disk. The Comptonization model nthcomp (Zdziarski et al. 1996;Życki et al. 1999) fits the coronal emission. The data to model ratios for this test are shown in Fig. 1. Common features for the plots in Fig. 1 are a broad line feature around 6-7 keV and a hump above 20 keV. These features are known as indications of a relativistic reflection component in the spectra (e.g. Ross et al. 1999). In some cases, e.g. GRS 1739-278, we also see a strong tail above 50 keV. This high energy excess may result from extra emission of the jets. Another explanation could be that the corona is a hybrid plasma with both thermal and non-thermal particles (e.g. Parker et al. 2015, Jiang et al. submitted). The non-thermal component can reduce the coronal temperature and produce more hard X-ray photons compared to the pure-thermal model (e.g. Coppi 1999;Fabian et al. 2017).
Then we model the reflection features with relativistic reflection models: constant * tbabs * ( cflux * nthcomp + cflux * relconv * reflionx HD). In this model, reflionx HD 12 is a rest-frame disk reflection model calculated with the reflionx code (Ross & Fabian 2005). The incident spectrum to the accretion disk is assumed to be described by nthcomp. Therefore, we link the photon index (Γ) and the coronal temperature (kT e ) parameters between nthcomp and reflionx HD. Other free parameters of the reflection model include the ionization parameter ξ = L/nr 2 (where L is the ionizing luminosity from the primary source, n is the density and r is the distance from the ionizing source), the iron abundance (A Fe ) and the electron density that can vary between log(n e /cm −3 ) = 15-22. The convolution kernel relconv (Dauser et al. 2010 is required to include  Greiner et al. (1996) Khargharia et al. (2010). Note that there is still debate on the distance of GS 1354-64 (also known as BW Cir). A smaller value has been reported by Gandhi et al. (2019).
the relativistic effects. In this way, the model calculates the angle-averaged spectrum with only some minor bias (smaller than the statistical errors) in the estimates of some parameters for low disk viewing angles . This model has parameters like the black hole spin (a * ), the disk inclination angle (i), the inner disk radius (R in ) and the emissivity profile. In this work, we fix the spin parameter in order to constrain the R in . The spin is fixed at a * = 0.2 for H 1743-322 (Steiner et al. 2012a), a * = 0 for IGR J17091-3624 (see Sec. B.3) and at the maximum value (0.998) for the other sources (see Sec. B for details). The inclination angle and the inner disk radius are left free during the fit. For the emissivity profile, we implement a broken power-law profile (i.e., ∝ 1/r qin for R in < r < R br and ∝ 1/r qout for R br < r < R out where R br is the breaking radius). In cases that q out and R br are not constrained, a powerlaw emissivity (q in = q out ) is used instead. We include a cflux model on each additive component to calculate the flux in the 0.1-100 keV band. Note that in some cases, a distant reflection component, a disk thermal emission component or additional absorption components are required to fit the data (see Fig. 10 and Sec. B for details).
To illustrate the impact on the measurements of spectral parameters after including a variable electron density, we also fit each spectrum with log(n e ) setting to 15, which is a value adopted by traditional reflection models. The best-fit parameters of the two sets of fittings are shown in Tab. 3 and Tab. 4. Moreover, we also test the relxill model (García et al. 2014) with log(n e ) = 15 (see Sec. A and Tab. 5). Fig. 2 shows an example of high-density reflection of GRS 1739-278. The best-fit electron density for this source is log(n e /cm −3 ) = 19.0 ± 0.4. Allowing the electron density to be a free parameter improves the χ 2 by 60 with one more degree of freedom (see Tab. 3 and Tab. 4). If fixing the density at log(n e /cm −3 ) = 15 (without refitting the spectrum), the reflection component would be depressed in the soft X-ray band, and there appears to be an excess in the residuals below 2 keV (see the third panel of Fig. 2). On the contrary, if increasing the density to log(n e /cm −3 ) = 22, the soft X-ray emission is apparently strengthened. This is a result of stronger free-free process on the disk when the density gets higher, which could also increase the disk surface temperature (see left panel of Fig. 6).

Comparisons with the low density model
In this section, we discuss the impact of the highdensity reflection model on spectral parameters compared to the traditional low density model. In Fig. 3, we show the measurements of the disk ionization parameter, the iron abundance and the reflection fraction. It is shown that when assuming a disk density of log(n e /cm −3 ) = 15, the ionization parameter is systematically overestimated. This is likely because increasing the ionization parameter produces more soft X-ray emissions in the reflected spectrum (see Fig. 3 of García et al. 2013), which mimics the effect of high-density reflection (see Fig. 4 of García et al. 2016).
When fitting the relativistic reflection models to X-ray spectra of XRBs and AGNs, a super-solar iron abundance is commonly required (e.g. Jiang et al. 2018;García et al. 2018a). Some studies have shown that highdensity models could lower the measured iron abundance (e.g. Tomsick et al. 2018;Jiang et al. 2019a). We do not see this effect in our sample (see the middle panel of Fig. 3) possibly because most of our observations do not need a super-solar iron abundance even in the case of log(n e /cm −3 ) = 15.
We define the reflection strength parameter (R str ) to be the energy flux ratio between the reflected and coronal emission in the 0.1-100 keV band. Fig. 3 shows that, except for V404 Cyg, the model with variable n e always returns a higher reflection strength. As evident in Fig. 2, additional thermalized emission appearing in the soft bands contributes additional emission which may explain the correlated increase in reflection strength with density.
It is of crucial importance to see if the assumption of the electron density would affect measurements of the inner disk radius. One of the reasons is that the black spin measurements with X-ray reflection spectroscopy rely on accurate determination of the ISCO size (e.g. Reynolds 2014). The constraints of R in of each observation are shown in Fig. 4 by plotting the distribution of ∆χ 2 as a function of the inner disk radius. It shows that allowing n e to be free to vary provides consistent measurements on R in compared to the case of fixed n e . This is expected  The red dashed and magenta dotted lines represent the cases when the density is set to 10 15 cm −3 and 10 22 cm −3 respectively without changing the other parameters. Lower three panels: The data to model ratios for the three cases in the upper panel. The blue, black and red data represent Swift-XRT, NuSTAR-FPMA and NuSTAR-FPMB respectively. since R in is mainly constrained by the red wing of the broad iron line as a result of the gravitational redshift effect, which does not depend on the disk density.

Disk densities
In the left panel of Fig. 5, we present our measurements of the disk density (log(n e )) versus log(M BHṁ 2 ). Previous studies in the literature of XRBs (the left cluster) and AGNs (the right cluster), which mainly cover the parameter space of log(M BHṁ 2 ) < 0 and log(M BHṁ 2 ) > 4, are plotted in grey. Our data fit into the gap 0 < log(M BHṁ 2 ) < 3. Overall, we find that XRBs require a disk density significantly higher than AGNs. Moreover, 16 out of our 21 spectra show a disk density higher than log(n e ) = 15.
According to Svensson & Zdziarski (1994), the density of radiation pressure-dominated standard disk (Shakura & Sunyaev 1973) follows: (1) where α = 0.1 is adopted for the viscosity parameter, σ T = 6.64 × 10 −25 cm 2 is the Thomson cross-section, R s is the Schwarzschild radius, R is the disc radius in units of R s ,ṁ is the dimensionless mass accretion rate defined asṁ = L bol /ηL Edd (η is the accretion efficiency that can reach 0.32 for a * = 0.998 and 0.057 for a * = 0, Thorne 1974), ξ is the conversion factor in the radiative diffusion equation that is chosen to be unity by Shakura & Sunyaev (1973) and f is the fraction of power transported from the disc to the corona. The solutions for different values of f are also plotted in the left panel of Fig. 5. We can see that most of the AGN data (the lower right cluster) can be explained by these solutions, with a f varying between 0.0 and 0.9 assuming ξ = 1. However, except for a few observations with the highest accretion rates, most disk density measurements of black hole XRBs are below the theoretical predictions.
We note that a large fraction (13/21) of our samples are in the range of L bol /L Edd < 10%. In this case, the gas pressure may play a significant role on the disk, especially when f is large (Svensson & Zdziarski 1994). According to Svensson & Zdziarski (1994) the radius at which the radiation pressure (P rad ) equals the gas pressure (P gas ) is determined by: (2) The left-hand side of this equation reaches minimum when R ≈ 5.7. Therefore, given certain values for M BH , f andṁ, the left-hand side can be always larger than the right-hand side in which case the disk is dominated by gas pressure. In the right panel of Fig. 5, we plot the thresholdṁ as a function of M BH for a few values of f . By fitting the AGN sample, Mallick et al. (2022) found f ≈ 0.7. We can see that the AGN sample is well above the threshold line for f = 0.7 and the radiation pressure should be important. However, this is not the case for our XRB sample, which indicates the important role of gas pressure. In the left panel of Fig. 5 we are also showing the solution for a gas pressure-dominated disk, which still predicts a disk density larger than our measurements by two orders of magnitude. One of the explanations could be that the reflection model only measures the density of the disk surface. The current reflection models assume a uniform density in the vertical direction, which may not be the case in reality. Developing reflection models that consider the vertical density structure of the disk would be important to understand this problem.
In Mallick et al. (2022), the authors find a balance between the incident radiation pressure (P rad = L/4πr 2 c = ξn e /4πc where L is the corona luminosity) and the disk gas thermal pressure (P th = n e k B T ). We explore this possibility with our XRBs sample. The radiation pressure can be directly calculated from the measured spectral parameters. To calculate the gas thermal pressure, we run the reflionx code to obtain the temperature at the Thomson depth τ = 1 for log(n e ) = 16 − 22 (see the left panel of Fig. 6). The inputs for the reflection calculation are chosen to be the average on our sample, i.e., Γ = 1.7, E cut =100 keV, ξ=100 erg cm s −1 and A Fe =1. The results are shown in the right panel of Fig. 6, where the dashed line shows P rad = P th . We can see that there is indeed a balance between the incident radiation pressure and the gas thermal pressure. Note that in this comparison we are neglecting the role of the magnetic field and the radiation pressure from the disk.

Disk inner radius
The inner disk radius (R in ) is an important parameter to understand the accretion process. The disk truncation     (Steiner et al. 2012a), a * = 0 for IGR J17091-3624 (see Sec. B.3) and a * = 0.998 for the other sources. The horizontal grey lines represent the case ∆χ 2 = 2.706 (90% confidence level for one parameter of interest). Data are offset in the vertical direction for visual clarity. The regions between ∆χ 2 -R in and ∆χ 2 = 2.706 are shaded with magenta (ne free) and grey (ne = 10 15 ) colors. The grey areas result after a mirror symmetry transformation with respect to the horizontal lines. For sources with multiple observations, the sequence follows Tab. 2 scenario has been commonly used to explain the distinctive states of XRBs (e.g. Esin et al. 1997). It is generally believed that the standard cold disk is truncated at large radius in the low hard state when the source luminosity is low (e.g. Tomsick et al. 2009) and it reaches the innermost stable circular orbit (ISCO) in the high luminosity soft state (e.g. Gierliński & Done 2004;Steiner et al. 2010). However, it is still under debate whether the disk is truncated in the bright hard state (e.g. Dzie lak et al. 2019; Mahmoud et al. 2019). We plot our reflectionbased measurement of the inner disk radius versus the Eddington-scaled luminosity in Fig. 7 along with previous measurements in the hard state of GX 339-4, which is one of the best-studied sources on R in (see also Fig. 8 of Wang-Ji et al. 2018). For GX 339-4, the measurements of R in from timing mode data of XMM-Newton EPIC/pn are systematically larger than those from NuSTAR or RXTE data. This may be due to the complex pile-up effects that are hard to eliminate (see the discussion in Wang-Ji et al. 2018).
Our measurements of R in from six black hole XRBs in the hard state are broadly consistent with the trend of GX 339-4 and we extend this tendency to the regime where L/L Edd > 20%. We can see that a small R in (< 10 R g ) is found when L/L Edd > 1% even though all the data in Fig. 7 are from the hard state. For L/L Edd > 10%, almost all measurements are consistent with R in being smaller than twice of the ISCO size (assuming a * = 0.998). These results support the idea that the inner disk radius can reach the ISCO in the bright hard state. Although it appears that the H 1743-322 disk is truncated at a large radius in the hard state, the black hole in this system has a low spin measured by the independent continuum-fitting method (a * ∼ 0.2, Steiner et al. 2012a) and thus the ISCO size is large (R ISCO = 5.3 R g , see the magenta dashed line in Fig. 7). Things are more complicated for IGR J17091-3624 since its black hole spin parameter is still uncertain (see Sec. B.3).

Inclination angle
In this analysis, we always fit the inclination as a free parameter. It allows us to compare the inclination angle measured from reflection modelling with that from other methods (often for the jet inclination or the binary inclination). For GS 1354-64, only an upper limit of 79 • has been obtained for its binary inclination through the absence of X-ray eclipses (Casares et al. 2009). There is a tight constraint on the binary inclination of V404 Cyg ((67 +3 −1 ) • , Khargharia et al. 2010) while the reflection analysis gives a lower value.
Similar discrepancies also have been found in other sources, e.g. for Cygnus X-1 the reflection spectra always require an inclination around 40 • (e.g. Tomsick et al. 2014;Walton et al. 2016;Liu et al. 2019) while its binary inclination is 27.51 +0.77 −0.57 (Miller-Jones et al. 2021). We note that the relativistic reflection spectra are only sensitive to the inner part of the accretion disk while the inner disk inclination may not align with the orbital inclination, which leads to a warped disk (e.g. Pringle 1996). Simulations have shown that the misalignment between the black hole spin axis and orbital angular momentum can be common for black hole XRBs (e.g. Brandt & Podsiadlowski 1994;Fragos et al. 2010). Observational indications of a warped disk have been found in a few sources, e.g. MAXI J1535-571 (Miller et al. 2018) and MAXI J1820+070 (Poutanen et al. 2022;Thomas et al. 2022). The warped disk may have an impact on the relativistic reflection spectra (e.g. Abarr & Krawczynski 2021) but a detailed analysis along this aspect is out of the scope of this work.
For H 1743-322, Steiner et al. (2012a) measured an inclination angle of 75 • ± 3 • for its large-scale ballistic jets, which are supposed to be aligned with the black hole spin axis. The inner disk region should also align with the spin axis due to the Bardeen-Petterson effect (Bardeen & Petterson 1975) while the reflection model gives a lower measurement of 30 • −40 • . We note that the inclination angle measurements with reflection analysis can be affected by systematic uncertainties, e.g. the lack of knowledge of the corona geometry or simplifications in the model calculations (see discssuions in Bambi et al. 2021). In some cases, the systematic uncertainty can be as large as ∼ 30 • (e.g. García et al. 2018b;Connors et al. 2022). This might explain the discrepancy we are seeing.

Reflection strength
Note that the reflection strength (R str ) defined in this work is not the reflection fraction (R f ) parameter in Dauser et al. (2016) which is the ratio between the corona intensity that shines on the disk and that reaches infinity. R f is determined by the accretion geometry while R str could be affected by inclination between the line-ofsight and the black hole spin axis (Dauser et al. 2016). We note that, in our sample, most sources are consistent with an inclination in the range of 20 • -30 • (except for MAXI J1535-571). Therefore, R str can still provide insight into the disk-corona geometry. In the left panel of Fig. 8, we are showing the relation between R str and the Eddington-scaled X-ray luminosity from our analysis. There is no strong correlation between the two parameters with a Pearson correlation coefficient (r) of 0.32 (1 − p = 75%). We also test the correlation between the photon index and the reflection strength (see the middle panel of Fig. 8) but still find no strong correlation (r = 0.22, 1 − p = 66%).

Photon index
It is well known in the literature that there is a positive statistical correlation between the Eddington ratio (λ Edd = L bol /L Edd ) and X-ray photon index (Γ) at high accretion rates in AGNs (e.g. Shemmer et al. 2008;Risaliti et al. 2009;Brightman et al. 2013) and XRBs (e.g. Kubota & Makishima 2004;Wu & Gu 2008;You et al. 2023). At low accretion rates (e.g. λ Edd < 1%, Yuan et al. 2007;Constantin et al. 2009), a negative correlation is often found (e.g. Gu & Cao 2009;Díaz et al. 2022). In the regime of extremely low accretion rates (e.g. λ Edd < 10 −4 ), Γ has been found to be saturated (e.g. Corbel et al. 2006;Sobolewska et al. 2011;Gültekin et al. 2012;Plotkin et al. 2013). This switch between correlation behaviors may suggest transition between different accretion modes (e.g. Cao et al. 2014;Yang et al. 2015). We test the correlation with our sample in the right panel of Fig. 8. Our data are all in the range of λ Edd > 1% and divided into two branches on the Γ-λ Edd diagram. In the lower branch, the data of V404 Cyg are not from canonical outbursts but strong repeated flaring events that last less than 1 ks (see Sec. B.4). Moreover, the flaring events might be related to the transient jet activities instead of changes of the mass accretion rate (Walton et al. 2017). These may explain why Γ stays stable even though the X-ray luminosity changes by more than one order of magnitude. Besides V404 Cyg, the data of GS 1354-64 are also deviating from the trend of other sources. If considering only the upper branch, there is a strong positive correlation between the two parameters (r = 0.83, 1 − p = 99.96%). Moreover, we plot in Fig. 8 the statistical relation found in AGNs by Shemmer et al. (2008), which is in good agreement with our upper branch. This indicates that the mass accretion rate changes the physical properties of the hot corona in a similar manner for XRBs and AGNs. 4.7. The corona properties in the compactness-temperature plane The coronal properties can be studied on the compactness-temperature ( − Θ) plane. The compactness is defined as: where L is the luminosity, R is the source radius assuming a spherical geometry, σ T the Thomson cross section and m e the electron mass. By studying a sample of NuSTAR measurements AGNs and black hole XRBs on the, Fabian et al. (2015) find that most of the measurements are clustered close to the limit set by electronpositron pair production (see also Ricci et al. 2018). This supports the idea that pair production is an important ingredient in the corona. This process could lead to runaway pair production to share the energy and thus limits the coronal temperature below the pair line (e.g. Svensson 1984;Zdziarski 1985;Stern et al. 1995;Coppi 1999). Moreover, the fact that many sources are above the electron-electron coupling line means that energetic electrons do not have enough time to thermalize,    which suggests the possibility of a magnetized and hybrid plasma that contains both thermal and non-thermal particles (e.g. Zdziarski et al. 1993;Grove et al. 1998). The inclusion of a non-thermal component would tend to lower the temperature of the thermal component and produce a high-energy excess compared to the spectra of purely thermal Comptonization. This behavior matches well with the data of GRS 1739-278 which has the lowest observed coronal temperature in our sample. Fabian et al. (2017) investigated the hybrid plasma scenario with a sample of AGNs and found that a non-thermal fraction of 10%-30% could account for most of the objects.
We plot the measurements of our sample on the compactness-temperature diagram along with the theoretical predictions of hybrid plasma to understand the physical properties of the corona in the hard state ( Fig. 9). The calculations of the theoretical lines are done with the eqpair 13 model (Coppi 1999). The nonthermal fraction is defined as the ratio between the compactness parameter of non-thermal and the total heating power ( nth / h ). For a fixed heating power, a higher non-thermal fraction results in a lower equilibrium temperature. We assume a corona size of 10 R g , which is a reasonable value given existing measurements on AGNs with X-ray reflection modelling or reverberation analysis (e.g. Fabian et al. 2009;Wilkins & Fabian 2011;Emmanoulopoulos et al. 2014). Fig. 9 shows that all of our measurements are be- 13 The description of the model can be found here: http://www. astro.yale.edu/coppi/eqpair/eqpap4.ps. The model is available in XSPEC: https://heasarc.gsfc.nasa.gov/xanadu/xspec/ manual/XSmodelEqpair.html. low the theoretical line of purely thermal plasma. A large fraction (15/21) of our data are below the line for nth / h = 30% which suggests an even higher nonthermal fraction. We note that the hybrid plasma model has been applied to a number of black hole XRBs (e.g. Gierliński et al. 1999;Zdziarski et al. 2001;Droulans et al. 2010;Parker et al. 2015;Zdziarski et al. 2021) and in some cases a significant non-thermal fraction is indeed required in the hard state (Wardziński et al. 2002;Nowak et al. 2011;Del Santo et al. 2016).
For individual sources, we can see that GS 1354-64 and IGR J17091-3624 follow the expected trend for a pairdominated plasma that the coronal temperature drops when it gets radiatively more compact. This is even true for the flux-resolved data of V404 Cyg although the timescale of the flaring activities we analyze here is as short as ∼ ks (see Sec. B.4). For H 1743-322, there are two outliers with higher temperatures compared to the other six observations with a similar level of radiation compactness. This may suggest a lower non-thermal fraction but it is also possible that the assumption of 10 R g for the corona size is too small for the two observations. The latter may not be the main reason since the measurements of R in for the two observations are consistent with 10 R g . We also note that the two outliers are from the so-called "failed" outbursts (e.g. Stiele & Kong 2021) during which the source stays in the low hard state only.
We further investigate these H 1743-322 outliers by considering the correlation between X-ray and radio luminosity. This source is known to switch between two different radio-X-ray correlations: the 'radio loud' and 'radio quiet' branches (Coriat et al. 2011). We therefore explore whether our inferred high/low non-thermal fraction corresponds to the source being on the radio loud/quiet branch, which may be expected if, for example, the non-thermal fraction is related to the jet properties. We find no such correspondance. Only one of our outliers has a simultaneous radio observation, and it belongs to the radio quiet branch (1.28 GHz flux density of 0.88 ± 0.05 mJy; Williams et al. 2020), as do several of the non-outliers.

CONCLUSION
In this work, we analyzed the broadband X-ray spectra of six stellar-mass black hole XRBs in the hard state with data from NuSTAR and Swift. The purpose is to test the effect of X-ray reflection by a high-density accretion disk. The main results are as follow: • The model with the disk electron density fixed at 10 15 cm −3 systematically overestimates the ionization degree of the disk atmosphere. The measurement of the inner disk radius is not affected by the assumption of disk density.
• The X-ray spectra of black hole XRBs require a higher disk density than that of AGNs. For the selected observations in this work, the accretion disk should be dominated by gas pressure. However, the measured densities are lower than either the prediction of gas or radiation pressure-dominated disk. The discrepancy might represent the vertical density structure of the disk. • From our analysis and previous measurements, we find that the inner disk radius can be close to the ISCO in the hard state.
• We find that the reflection strength is not correlated to the Eddington ratio or photon index. There is a strong correlation between the Eddington ratio and photon index which is in good agreement with the statistical relation found in AGNs.
• The coronal temperature is lower than the prediction of a purely thermal plasma and can be explained by a hybrid plasma with a non-thermal fraction around 30%.
(a * = 0.8 ± 0.2), an intermediate inclination (20 • -40 • ) and an inner disk radius close to the ISCO. We obtain similar measurements with our reflection modelling (see Tab 4 and Tab 5). The very low temperature of the corona has already been revealed by Miller et al. (2015). There is another NuSTAR observation in the archive (obsID: 80101050002) from September 2015 but the reflection features are too weak (Fürst et al. 2016).
GS 1354-64 GS 1354-64 is a Galactic black hole candidate discovered in 1987 (Makino 1987). Only a lower limit (∼ 8 M ) has been determined for its black hole mass and its distance is also not well-constrained (25-61 kpc, Casares et al. 2009). The first two observations in June and July in our analysis have been studied by El-Batal et al. (2016) without including the Swift data. The authors found a large truncation of the accretion disk for the June observation with a low density relxill model, which is consistent with our results in Tab. 5. The high disk inclination and small inner disk radius measured from the July observation are also in agreement with previous studies (e.g. El-Batal et al. 2016;Xu et al. 2018a).
IGR J17091-3624 IGR J17091-3624 is a Galactic black hole candidate at a distance of 11 kpc to 17 kpc (Rodriguez et al. 2011, but see Altamirano et al. 2011 for debates on its distance.). The black hole spin parameter is still quite uncertain. Previous studies have reported both negative spin (Rao & Vadawale 2012;) and high spin (a * > 0.9, Reis et al. 2012a) measurements. The three NuSTAR and Swift observations in this analysis have been studied by Xu et al. (2017) with relxill-based low density reflection models. By fixing a * = 0.998, the authors found a truncated disk (R in ∼ 20R g ) and an intermediate disk inclination angle (i ∼ 30 • − 40 • ), which are confirmed with our analysis (see Tab. 5).

V404 Cyg
V404 Cyg is a dynamically confirmed black hole XRB located at 2.39 ± 0.14 kpc (Miller-Jones et al. 2009). The black hole mass should be 9.0 +0.2 −0.6 M ( Khargharia et al. 2010). The two observations of V404 Cyg analyzed in this work show strong flaring variability (see Fig. 2 of Walton et al. 2017) with the count rate (per FPM) changing from 100 to more than 10 4 ct s −1 . It was shown by Walton et al. (2017) that there are strong relativistic reflection features in the X-ray spectra and it is possible to constrain the black hole spin parameter if dividing the data according to flux levels and neglecting the data with deep absorption edge. We process the data in the same way as Walton et al. (2017) and obtain the spectra for five flux levels (see Tab. 1 of Walton et al. (2017)).
To fit the spectra of V404 Cyg, it requires to include additional components: a neutral local absorption layer, an ionized absorption layer and a distant reflection component. The ionized absorption is modeled with a grid calculated with XSTAR (Kallman & Bautista 2001;Kallman et al. 2004). The full model of our Model 2 for this source reads as: constant * tbabs gal * (tbabs local * xstar * (cflux * nthcomp + cflux * relxillcp) + cflux * xillvercp). The column density for the Galactic absorption is fixed at 10 22 cm −2 and the absorption column locally to the source is a free parameter. The ionization state of the distant reflection component (xillvercp) is fixed at log(ξ) = 0 and the other parameters are linked to the relativistic reflection component. H 1743-322 H 1743-322 was discovered in 1977(Kaluzienski & Holt 1977. It is located at 8.5 ± 0.8 kpc and the jet inclination is 75 • ±3 • (Steiner et al. 2012a). The continuum fitting method indicates a black hole spin around a * ∼ 0.2 (Steiner et al. 2012a). Molla et al. (2017) estimated the black hole mass to be 11.21 +1.65 −1.96 M using the two-component advective flow (TCAF) solution and the correlation between photon index and frequency of the quasi-periodic oscillation (QPO).
There are eleven NuSTAR observation in the archive, three of which do not show reflection features. The other eight observations analyzed in this work are all in the hard state. The two observations from 2016 have been analyzed by Chand et al. (2020) and a truncated disk was found. However, our analysis with relxillcp model requires a disk extends to the ISCO. This difference is possibly due to the fact that we treat the inclination angle as a free parameter while in Chand et al. (2020) it is fixed at 75 • . The two 2018 observations have also been analyzed with reflection models assuming a fixed inclination angle (Stiele & Kong 2021) and again their measurements on R in are larger than ours.
MAXI J1535-571 MAXI J1535-571 is a black hole candidate discovered in 2017 (Negoro et al. 2017). It is located at a distance of 4.1 +0.6 −0.5 kpc (Chauhan et al. 2019). Previous reflection analysis indicates a high black hole spin and high inclination angle (e.g. Xu et al. 2018b;Miller et al. 2018). Following Xu et al. (2018b), we include a disk component and a distant reflection component to the full model to fit the NuSTAR data. The high inclination angle is well recovered with our fit (see Tab 5).

2222.3/2157
Note: The flux of the power-law and reflection component are estimated in the 0.1-100 keV band. The symbol P denotes the lower or higher limits. Note that for V404 Cyg, the column density (nH) of the Galactic absorption is fixed at 1 × 10 22