A NuSTAR view of GRS 1716−249 in the hard and intermediate states

ABSTRACT

We present a detailed analysis of the spectral properties of the black hole transient GRS 1716−249, based on the archival Swift and NuSTAR observations taken during the outburst of this source in 2016–2017. The first six NuSTAR observations show that the source is in a canonical hard state, where the spectrum is dominated by a power-law continuum. The seventh NuSTAR observation is taken during the intermediate state where both a disc thermal component and a power-law continuum are shown. All of our observations show a broad emission-line feature in the iron band and a Compton hump above 10 keV. We model the broad-band spectra using a high-density disc reflection model, where the soft X-ray emission in the hard state is interpreted as part of the disc reflection component. This model enables us to constrain the disc density parameter of GRS 1716−249 in the range of 10¹⁹–10²⁰ cm⁻³. We only obtain an upper limit of the inner disc radius using high-density disc reflection spectroscopy and the results indicate either a non-truncated disc or a slightly truncated disc with R_in ≲ 20 r_g.

1 INTRODUCTION

X-ray binaries (XRBs) are accretion-powered binary systems that are luminous in X-rays, and typically made up of a compact object accreting from a ‘normal’ stellar companion. Many black hole (BH) XRBs show transient events that are characterized by outbursts in the X-ray band. Some outbursts last for only few weeks (e.g. MAXI J1659−152, Muñoz-Darias et al. 2011) while some last for years before they return to the quiescent state (e.g. Swift J1753.5−0127, Soleri et al. 2013). Two main different spectral states are commonly seen during the outbursts of BH XRBs: (1) the hard state where the spectrum is dominated by a power-law-shaped continuum; (2) the soft state where the spectrum is dominated by a disc-blackbody emission (e.g. Oda et al. 1971). In the hard state, the spectrum often shows broad Fe K emission line and a strong Compton hump above 10 keV in addition to the power-law emission that is associated with the corona (e.g. Kitamoto et al. 1990). These features are very similar with the ones found in Seyfert galaxies despite different BH mass scales (e.g. Walton et al. 2012) and together are referred to as the disc reflection component. In some cases, BH XRBs show intermediate states between the soft and hard states, where a hard power-law continuum and a disc thermal emission component make approximately equal contribution to the total spectrum (e.g. Steiner et al. 2012; Tomsick et al. 2018).

Previous reflection-based spectral analysis of BH XRB observations commonly found a supersolar disc iron abundance (e.g. García et al. 2018). A possible solution to this problem is a high-density disc reflection model (e.g. Tomsick et al. 2018; Jiang et al. 2018, 2019b). However, a fixed disc density (n_e = 10¹⁵ cm⁻³) was typically assumed in previous models. This assumption is appropriate for the disc of a massive supermassive BH (e.g. BH mass m_BH = M_BH/M_⊙ > 10⁸ and mass accretion rate |$\dot{M}\gt 20{{\ \rm per\ cent}}$||$\dot{M}_{\rm Edd}$|⁠, Shakura & Sunyaev 1973). At higher disc densities, the temperature of the disc surfaces increases due to stronger free–free absorption. The soft X-ray band of the disc reflection spectrum turns into a blackbody-shaped emission (Ross & Fabian 2007; García et al. 2016). When modelling the broad-band spectra with a high-density reflection model, the continuum flux level of the disc reflection component increases and the inferred iron abundance parameter drops for the spectral fitting purpose (Jiang et al. 2019a). The blackbody-shaped emission feature in the high-density disc reflection model may be able to explain the weak blackbody components commonly seen in the hard-state spectra of BH XRBs (Jiang et al. 2019b) and the blackbody-shaped soft excess emission in Seyfert galaxies (e.g. Mallick et al. 2018; García et al. 2019).

GRS 1716−249 is a low-mass XRB first discovered in 1993 (Ballet et al. 1993; Harmon et al. 1993). It reached a peak flux of 1.4 Crab in the 20–100 keV band within five days. By studying the Na D absorption lines, della Valle, Mirabel & Rodriguez (1994) obtained an upper limit of d ≈ 2.8 kpc for the distance of GRS 1716−249. In our work, we adopt d = 2.4 kpc, the same value applied for luminosity calculations in della Valle et al. (1994). Masetti et al. (1996) estimated a lower limit of m_BH > 4.9 for the mass of the central BH by studying the period of the light modulations during the previous outburst. GRS 1716−249 was then detected in another outburst in 1995 with a simultaneous radio flare (Hjellming et al. 1996).

GRS 1716−249 was found in an outburst 21 yr later by MAXI on 2016 December 18 (Masumitsu et al. 2016; Negoro et al. 2016). By analysing the spectra from the Swift monitoring programme of this outburst, Bassi et al. (2019) found that the source approached the soft state three times but never reached the canonical soft state. Simultaneous radio observations also show GRS 1716−249 is located on the radio-quiet branch of the X-ray–radio luminosity plane (Bassi et al. 2019).

In this paper, we conduct detailed spectral analysis of the seven NuSTAR (Harrison et al. 2013) observations of GRS 1716−249 triggered during the outburst in 2016–2017. Simultaneous Swift snapshot observations in the soft X-ray band are considered as well if available. The first three NuSTAR observations did not have simultaneous Swift observations because of the sun violation for Swift. In Section 2, we introduce the data reduction processes; in Section 3, we present the broad-band spectral modelling by using high-density disc reflection model, study the spectral variability, and conduct multi-epoch spectral analysis to obtain the disc viewing angle and the disc iron abundance; in Section 4, we discuss our spectral analysis results; in Appendix, we study the parameter degeneracy when both the disc thermal and the disc reflection components are considered for the hard-state observations of GRS 1716−249.

2 DATA REDUCTION

2.1 NuSTAR data reduction

We use the standard pipeline nupipeline V0.4.6 in heasoft V6.25 to reduce the NuSTAR data. The calibration file version is v20171002. In addition to the science mode (mode 1) data, we also extracted the products from the spacecraft science (mode 6) event files by following Walton et al. (2016). The mode 6 events are considered because the star tracker was blinded by the sun during the early observations of GRS 1716−249 and was unable to provide an aspect solution. Therefore, the source coordinate in the sky was calculated using the spacecraft bus attitude solution for mode 6 observations. The inclusion of the mode 6 data maximizes the usable exposure of our observations. See Table 1 for the exposure time in two different modes. The mode 1 source spectra are extracted from a circular region with radii of 120 arcsec, and the background spectra are extracted to be circular regions near the source region on the same chip. The mode 6 source spectra are extracted from a circular region with radii of 110–140 arcsec. We generate spectra using the NUPRODUCTS command. The FPMA and FPMB spectra are grouped to a minimum signal-to-noise ratio (S/N) of 6 and to oversample the energy resolution by a factor of 3. The NuSTAR data are modelled over the full 3–78 keV band in this work.

2.2 Swift data reduction

We use the standard pipeline xrtpipeline V0.13.5 to reduce the Swift XRT data (Burrows et al. 2005). The calibration file version is V20171113. The source spectra are extracted from an annulus region to avoid pile-up effects. The inner radius of the annulus is 5 arcsec and the outer radius of the annulus is 40 arcsec. The XRT spectra are grouped to a minimum S/N of 6 and to oversample the energy resolution by a factor of 3. The Swift data are modelled over the 1–8 keV band in this work. We ignore the spectra below 1 keV due to the known calibration issue with the Window Timing mode.

3 DATA ANALYSIS

In this section, we first present an overview of the spectral transition through the outburst of GRS 1716−249 in 2016–2017. Second, we conduct detailed broad-band spectral modelling using a high-density disc reflection model. In the end, we discuss the changes of the disc-coronal geometry and the disc properties during the outburst by conducting a multi-epoch spectral analysis.

3.1 Spectral state transition

The observed NuSTAR and Swift flux of GRS 1716−249 is shown in the top panel of Fig. 1. The MAXI light curves of the same period in different energy bands are shown in the middle panel. The outburst of this source lasted for 8 months until it returned to the quiescent state. The X-ray flux rose for the first 2 months of the outburst before starting a slow decay. The bottom panel of Fig. 1 shows the hardness ratio between the 4–10 and 2–4 keV bands. The hardness ratio curve from MAXI indicates that the source spent majority of the outburst in the hard state before approaching the soft state near the end of the outburst, although the total X-ray flux continued to decay across this period. The decrease of the hardness ratio is also confirmed by the NuSTAR and Swift observations. GRS 1716−249 returned to the hard state before fading back to quiescence, indicating an outburst with failed transition to the soft state. We refer interested readers to Bassi et al. (2019) for a complete hardness–intensity diagram (HID) of this outburst obtained by Swift, which monitored GRS 1716−249 throughout this outburst.

The black arrows in Fig. 1 mark the dates when the NuSTAR observations were taken during the outburst. The spectra of all the NuSTAR and coordinated Swift observations unfolded through a constant model are shown in Fig. 2. The left-hand panel of Fig. 2 shows the spectra of the first four observations during the rising phase of the outburst. These four observations show very similar spectral shape but different flux levels. The right-hand panel shows the spectra extracted from the rest of the observations. Obs 5 and 6, which were taken 2 months after the fourth observation, show a softer broad-band spectral shape but a similar X-ray flux level compared with the observations in the rising phase of the outburst. Such a change in the spectral shape matches the indication from the MAXI light curves in Fig. 1. Obs 7 was taken 7 months after the beginning of the outburst. The spectra of obs 7 are shown in the right-hand panel of Fig. 2 for comparison with previous observations. The last observation shows the lowest total X-ray flux level but a very strong blackbody-shaped soft emission below 3 keV.

3.2 Spectral modelling

We use xspec V12.10.1b (Arnaud 1996) for spectral analysis and C-stat is considered in this work. The upper limit of the Galactic column density in the line of sight towards GRS 1716−249 is 3.97 × 10²¹ cm⁻² (Willingale et al. 2013). We fix the column density of the Galactic absorption at this value when fitting the spectra with absorbed power-law models. Then we model with the column density as a free parameter for the final results.

We first fit the spectra with absorbed power-law models. Data/model ratio plots for these models are shown for each epoch in Fig. 3. All the spectra show a broad emission-line feature in the iron band and a Compton hump above 10 keV. In addition to these relativstic reflection features, the soft energy band of obs 7 is dominated by a strong excess emission. Based on the evidence for broad Fe K emission and the Compton hump in each of the spectra, we model the broad-band spectra with a combination of a coronal emission component and a relativistic disc reflection component. The former component is modelled using the simple power-law model cutoffpl with a high-energy exponential cut-off. The latter component is calculated using relconvlp * reflionx. A more developed version of reflionx (Ross & Fabian 2005, 2007) is used to model the disc reflection model in the rest frame with following parameters: disc iron abundance (Z_Fe), disc ionization (log (ξ)), disc density (n_e), high energy cut-off (E_cut), and photon index (Γ). The solar abundance in reflionx is calculated by Morrison & McCammon (1983). The ionization parameter is defined as ξ = 4πF/n in units of erg cm s⁻¹, where n is the hydrogen number density and F is the illuminating flux. The high-energy cut-off parameter and the photon index in the reflection model are linked to the corresponding parameters of the cutoffpl component. The convolution model relconvlp (Dauser et al. 2013) is used to apply relativistic effects to the rest-frame reflection model. This calculates the emissivity profile according to a simple lamp-post model parametrized by the height of the corona above the disc. Other free parameters in this model are the disc inner radius r_in/r_ISCO and the disc inclination angle i. The spin parameter is fixed at the maximum value 0.998 to fully study the inner disc radius parameter. The outer radius of the disc is assumed to be at 400 r_g as the fits are insensitive to this parameter. The convolution model cflux is used to calculate the flux of each model component between 1 and 78 keV band. We define the reflection fraction as f_refl = F_refl/F_pl, where F_refl and F_pl are the best-fitting flux values given by cflux. Note that this empirical reflection fraction is different from the physical reflection fraction discussed in Dauser et al. (2016). A constant model is used to account for the cross-calibration uncertainity between different instruments. The total model is constant * tbabs * (cflux * (relconvlp * reflionx) + cflux * cutoffpl) (MODEL1) in the xspec format.

The combination of a relativistic disc reflection and a power-law-shaped coronal emission models offers a very good fit to the spectra of first six observations with C-stat/ν ≈ 1.0, where ν is the number of degrees of freedom. For example, we obtain a good fit with no evidence for structured residuals and C-stat/ν = 2665.32/2613 by modelling the spectra of obs1 with MODEL1. In order to test for any narrow emission line in the iron band, we also tried adding one Gaussian line model gauss with the line energy fixed at 6.4 keV. This additional line model only improves the fit by ΔC-stat = 2 with two more free parameters for obs 1. Only an upper limit of the equivalent width of the line model is obtained (EW < 3 eV). Similar conclusion can be found when analysing the spectra of obs 2–6. Therefore, we conclude that there is no evidence for a distant reflector in our observations.

As shown in Figs 2 and 3, the low-energy spectrum of obs 7 shows clear evidence for an additional blackbody-shaped component. By modelling the Swift and NuSTAR spectra simultaneously with MODEL1, we fail to obtain a good fit for obs 7. The high-density disc reflection model is unable to model the strong ‘soft excess’ emission with C-stat/ν = 5250.68/1512. By adding an additional diskbb component, we improve the fit by ΔC-stat = 3486 with 2 more parameters. The full model is constant * tbabs * (cflux * (relconvlp * reflionx) + cflux * cutoffpl + diskbb) (MODEL2) in the xspec format. The requirement for an additional disc-blackbody component agrees with previous analysis (e.g. Armas Padilla & Munoz-Darias 2017; Bassi et al. 2019), marking a transition to an intermediate state where both the disc thermal emission and the coronal emission contribute significantly to the total observed X-ray flux. Similarly, a combination of high-density disc reflection and thermal emission from the disc has been also found in the intermediate state of Cyg X-1 (Tomsick et al. 2018).

In order to test whether the spectra of the hard-state spectra (obs 1–6) require an additional diskbb component, we re-analysed the spectra with MODEL2. MODEL2 provides an equivalent fit as MODEL1 with C-stat decreasing by 2 and two more free parameters for obs 4 as an example. See Appendix for more details about our further Markov chain Monte Carlo (MCMC) analysis for the hard-state observations. We conclude that an additional diskbb component is statistically not required for the spectral fitting for the hard-state observations as it does not improve the spectral fitting significantly. However, we do note that there is a degeneracy between the inner disc temperature parameter (kT_e) of the diskbb model, the disc density parameter in the reflection model, and the column density of the line-of-sight absorption. Similar conclusion has been found in the previous analysis of the hard-state observations of GX 339−4 with a variable disc density parameter (Jiang et al. 2019b).

We also notice that the FPMA and FPMB data of obs 2 and 7 show difference of 6 per cent in the spectra below 4 keV. Similar disagreement between FPMA and FPMB spectra has also been found in other NuSTAR observations (e.g. Madsen et al. 2017) due to possible calibration uncertainty.

3.3 Multi-epoch spectral analysis

In the previous section, we successfully model the spectra of obs 1–6 with MODEL1 and obs 7 with MODEL2 individually. The best-fitting disc inclination angles i for all the epochs are consistent with i ≈ 30° and the disc iron abundances Z_Fe are consistent with the solar iron abundance. Therefore, we conduct a multi-epoch spectral analysis by linking the disc iron abundance and the disc inclination angle parameters. MODEL1 is used for obs 1–6 and MODEL2 is used for obs 7. The best-fitting model for each observation is shown in Fig. 4 and the best-fitting parameters are shown in Table 2. We obtain a very good fit with C-stat/ν = 10299.69/9299. A small disc viewing angle of |$i = {31^{+3}_{-2}}^{\circ }$| and a near-solar iron abundance (⁠|$Z_{\rm Fe} = 1.1\pm 0.5\, \mathrm{Z}_{\odot }$|⁠) are found. All the observations require a disc density significantly higher than n_e = 10¹⁵ cm⁻³, which was assumed in most of previous reflection-based analysis for other XRB.

Fig. 5 presents the changes of the best-fitting parameters with the Eddington ratio. |$L_{\rm 0.1-100\, keV}$| is the 0.1–100 keV band absorption-corrected luminosity calculated using the best-fitting model. Only a lower limit of m_BH is currently available (e.g. m_BH > 4.9, Masetti et al. 1996). A BH mass of 10|$\, \mathrm{M}_{\odot }$| is assumed for the calculation of the Eddington luminosity for simplicity. The relative luminosity change between epochs does not change even if this mass assumption is incorrect. Because the same mass value is used for the calculation of the Eddington luminosity for all the epochs.

The first panel of Fig. 5 shows the best-fitting column density of the X-ray absorption. The upper limit of the Galactic column density along the line of sight is estimated being 3.97 × 10²¹ cm⁻² (Willingale et al. 2013). Most of the best-fitting values obtained in this work are consistent with this value. The column densities obtained for obs 3 and 7 are significantly higher than 4 × 10²¹ cm⁻², possibly indicating a variable line-of-sight neutral absorption. More work concerning possible degeneracy between the column density and other parameters in this analysis can be found in Appendix when the thermal emission from the disc is considered in the model for the hard-state observations.

The next four panels of Fig. 5 present the evolution of the coronal properties through the outburst. The absolute flux of the coronal emission (F_pl) increases with the total X-ray flux. The photon index of the coronal emission remains consistent with Γ = 1.4 in obs 1–4, although F_pl increases by a factor of 2. As shown in Fig. 2, the coronal emission becomes softer in obs 5 and 6 with Γ = 1.5. The spectra show the softest continuum in obs 7 when the source is in an intermediate state. The high-energy cut-off parameter (E_cut) shows an anticorrelation with the Eddington ratio, except for obs 4 when the source is in the highest flux state. Obs 6, where a softer continuum is found compared to obs 1–4, has the lowest E_cut, corresponding to a cooler corona. In comparison, obs 7, where the source is in the lowest flux state, has the highest E_cut, corresponding to a hotter corona. This suggests the temperature of the corona may respond to the accretion rate of the disc. At a higher luminosity/accretion rate, more disc photons are upscattered in the coronal region, which increases the radiative cooling and lowers the coronal temperature. The fifth panel shows the height of the corona h for each epoch. Only upper limits of h are found. No obvious correlation between h and the X-ray luminosity is found. Future observations with a higher spectral resolution in the iron band may be able to better constrain the geometry of the corona in BH transients by measuring the emissivity profile of the broad Fe K emission line in more detail.

The variability of the disc reflection is shown in the last three panels of Fig. 5. The absolute flux of the disc reflection component remains consistent when |$L_{\rm 0.1-100\, keV}/L_{\rm Edd}\lt 1{{\ \rm per\ cent}}$| and increases when |$L_{\rm 0.1-100\, keV}/L_{\rm Edd}\gt 1{{\ \rm per\ cent}}$|⁠. Due to the decrease in the absolute flux of the coronal emission, obs 7 shows the highest reflection fraction in all the epochs. Only upper limits of the inner radius have been obtained in all the observations. However, a larger upper limit is found at low luminosities, indicating a possible disc truncation (e.g. r_in ≈ 20r_g) at |$L_{\rm 0.1-100\, keV}/L_{\rm Edd} = 0.5{{\ \rm per\ cent}}$|⁠.

4 DISCUSSION

We have undertaken an analysis of the broad-band (1–78 keV) X-ray spectra of the BH XRB GRS 1716−249 taken by NuSTAR and Swift during its recent ‘failed’ outburst. Following our recent work on GX 339−4 (Jiang et al. 2019b), we focus on modelling the data with reflection models that allow the density of the disc to be fit as a free parameter (see also Tomsick et al. 2018).

We first summarize the spectral fitting results: the first four observations were taken at the rising phase of the outburst. The spectra show a consistent spectral shape but an increasing X-ray luminosity (from 0.8 per cent to 1.2 per cent of L_Edd). The inner radius of the disc is found to be <8r_ISCO in obs 4 when the source reaches the highest X-ray luminosity, indicating a disc extending to the ISCO or a slightly truncated disc. The spectra of obs 5 and 6 become softer than obs 1–4 although the source remains at a similar X-ray luminosity as the previous observations. During this transition, obs 6 shows the lowest high-energy cut-off in the coronal emission, indicating a cool coronal temperature (e.g. T = 30–50 keV¹). At the same time, the density in the surface of the disc increases by a factor of 5 compared to the density in obs 1–4. The last observation, obs 7, was taken 7 months after the beginning of the outburst. The spectra of obs 7 show a combination of a disc thermal emission (kT = 0.4keV), a softer coronal emission (Γ = 1.76), and a strong disc reflection component, indicating the source is in a very low flux intermediate state. GRS 1716−249 shows the lowest X-ray luminosity in obs 7 compared to obs 1–6, but has a similar absolute flux from the disc reflection component to these higher flux observations, indicating a possible change in the geometry of the corona. However, only upper limits of the coronal height are obtained.

Second, we discuss the disc density measurements during the outburst of GRS 1716−249. The disc densities measured in this work are shown in brown in Fig. 6. The green circles represent the observations of GX 339−4 in the hard state, and the green circle represents the high flux soft-state observations of the same source in 2015 (Jiang et al. 2019b). The blue triangle shows the disc density of Cyg X-1 in the intermediate state (Tomsick et al. 2018). The mass accretion rate |$\dot{m}$| is estimated using |$\dot{m} = L_{\rm Bol}/\epsilon L_{\rm Edd}\approx L_{\rm 0.1-100\, keV}/\epsilon L_{\rm Edd}$|⁠, where ϵ is the accretion efficiency. A BH mass of |$10\, \mathrm{M}_{\odot }$| and an accretion efficiency of |$\epsilon = 20{{\ \rm per\ cent}}$| are assumed in our calculation. Assuming the distance measurements of GX 339−4 and GRS 1716−249 are accurate and the BH masses used for these calculations are correct, the hard states of two sources show a similar Eddington ratio. However, the disc density in GRS 1716−249 is 10 times lower than the density in GX 339−4 (Jiang et al. 2019b). The difference of the disc density could be due to the different BH mass, the accretion efficiency, and the true intrinsic bolometric luminosity. The intermediate state of GRS 1716−249 shows a similar disc density as the intermediate state of Cyg X-1 but a much lower luminosity in the X-ray band. Nevertheless, we conclude that the disc density required for the broad-band spectral analysis of GRS 1716−249 is significantly larger than n_e = 10¹⁵ cm⁻³, which was assumed in previous reflection-based spectral modelling.

The blue solid line in Fig. 6 shows the disc density at r = 2R_S of a radiation pressure-dominated disc varying with the parameter |$\dot{m}$|⁠. r_in = 1R_S is assumed for the disc inner radius. We refer interested readers to Svensson & Zdziarski (1994) for more details and Jiang et al. (2019b) for the complete equations. In this thin disc model, the density of the disc has an anticorrelation with the BH mass accretion rate (⁠|$n_{\rm e} \propto \dot{m}^{-2}$|⁠). However, a significantly lower disc density has been found in various BH transients (e.g. Tomsick et al. 2018; Jiang et al. 2019b) compared to the radiation pressure-dominated disc model. The reasons of the offset might be (1) the uncertainties of the BH mass and distance measurements. This may also be able to explain a lower disc density in GRS 1716−249 than in GX 339−4 due to the underestimation of |$\dot{m}$|⁠. We note, however, that the distance and the BH mass of Cyg X-1 is well constrained but the density given by reflection spectroscopy is still lower than the classical disc model; (2) the low accretion rate (⁠|$\dot{m}\lessapprox 1{{\ \rm per\ cent}}\, \dot{m}_{\rm Edd}$|⁠) in the hard state suggests that the disc might be dominated by gas pressure in the inner region instead of radiation pressure (Shakura & Sunyaev 1973; Svensson & Zdziarski 1994). The vertical structure of a gas pressure-dominated disc can be estimated by assuming an isothermal model (e.g. Frank, King & Raine 2002). For a stellar-mass BH, such as GX 339−4 and GRS 1716−249, the disc density within one Thomson optical depth in the surface of the disc is approximately 100 times lower than the density in the mid-plane (Jiang et al. 2019a). However, the disc density in the standard disc density model shown in Fig. 6 is calculated assuming a uniform density across the height of a thin disc (Shakura & Sunyaev 1973; Svensson & Zdziarski 1994). This assumption has been proven to be invalid in later studies of discs (e.g. Ionson & Kuperus 1984; Turner 2004).

ACKNOWLEDGEMENTS

JJ acknowledges support by the Cambridge Trust and the Chinese Scholarship Council Joint Scholarship Programme (201604100032), the Tsinghua Shuimu Scholar Programme and the TAO Fellowship. DJW acknowledges support from an STFC Ernest Rutherford Fellowship. ACF acknowledges support by the ERC Advanced grant 340442. MLP and FF are supported by European Space Agency (ESA) Research Fellowships. This work was based on the reflionx code written by Randy R. Ross. We acknowledge his big contribution to this project. This work made use of data from the NuSTAR mission, a project led by the California Institute of Technology, managed by the Jet Propulsion Laboratory, and funded by NASA. This research has made use of the NuSTAR Data Analysis Software (NuSTARDAS) jointly developed by the ASI Science Data Center and the California Institute of Technology. This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester.

Footnotes

REFERENCES

APPENDIX: CONSTRAINING THE DISC THERMAL COMPONENT

Previous spectral modelling of BH transients in the hard state commonly requires an additional weak disc-blackbody component when assuming a fixed disc density parameter (n_e = 10¹⁵ cm⁻³) for the reflection model. The origin of such a disc-blackbody component remains unclear. For example, the hard state of GX 339–4 shows |$L_{\rm X} \lessapprox 1{{\ \rm per\ cent}} L_{\rm Edd}$| and the spectral modelling suggests the presence of a disc-blackbody component with a temperature of kT > 0.8 keV (e.g. Wang-Ji et al. 2018). Such a high temperature is comparable to the one in the soft state of the same source, where the X-ray band is dominated by the disc thermal component and has a luminosity of |$L_{\rm X} \approx 20{{\ \rm per\ cent}} L_{\rm Edd}$| (e.g. Parker et al. 2016; Jiang et al. 2019b).

In this work, the blackbody-shaped emission shown in the hard state of GRS 1716−249 is instead interpreted as part of a high-density disc reflection component. A disc density of 10¹⁹–10²⁰ cm⁻³ is found using high-density disc reflection spectroscopy. See Section 3.3 for more details. Similarly, Jiang et al. (2019b) analysed the NuSTAR and Swift spectra of GX 339−4 in the hard state and found that the disc thermal component is statistically degenerate with the column density when both disc thermal component and high-density disc reflection are considered. In this section, we apply the same method as in Jiang et al. (2019b) to the hard-state spectra of GRS 1716−249. Both relativistic disc reflection and disc thermal components are considered (MODEL2). We do not consider obs 1–3 as they do not have soft X-ray coverage. Only obs 4 is shown here as an example. Similar conclusions can be found for the rest of the hard-state observations. Note that the Swift and NuSTAR spectra of obs 4 cover 1–78 keV. The Swift spectrum below 1 keV is ignored due to calibration issue.

The best-fitting parameters given by MODEL2 are shown in Table A1, in comparison with the parameters given by MODEL1 where only reflection is considered. MODEL2 gives a slightly worse fit compared to MODEL1 with ΔC-stat = 5 and two more free parameters. In conclusion, both two models offer similar good fits to the spectra. However, only upper limits of the disc density and the disc thermal component are obtained when an extra diskbb model is added. The column density shows a larger measurement uncertainty calculated by using the ERROR command in xspec when MODEL2 is used, although they are consistent with the results given by MODEL1. Different uncertainty ranges given by MODEL2 likely indicate parameter degeneracy when both reflection and disc thermal component are considered. Therefore, we checked the constraints of all the parameters by using the MCMC algorithm. The xspec/emcee code by Jeremy Sanders based on the python implementation (Foreman-Mackey et al. 2013) and the MCMC ensemble sampler (Goodman & Weare 2010) was used. We use 100 walkers with a length of 400 000, burning the first 5000. A convergence test has been conducted and the Gelman–Rubin scale-reduction factor R < 1.3 for every parameter. Fig. A1 shows the output distributions of some of the parameters that affect the soft X-ray spectral modelling.

From our MCMC analysis, we find that:

• The column density N_H of the line-of-sight absorption towards GRS 1716−249 is well constrained between (0.5–6) × 10²¹ cm⁻² (3σ). This result is consistent with the upper limit of the Galactic absorption towards GRS 1716−249 (e.g. N_H = 4.7 × 10²¹ cm⁻², Willingale et al. 2013). However, N_H shows degeneracy with the temperature of the disc blackbody component kT and the disc density parameter n_e. When kT and n_e are higher, the continuum model produces a stronger emission in the soft X-ray band and thus a higher N_H is required to fit the spectrum.

• When kT is lower, a higher value of n_e is required to model the blackbody emission in the soft X-ray band. For example, our MCMC analysis shows that (1) n_e is in the range of 10¹⁶–10²⁰ cm⁻³ (3σ) when kT < 0.4 keV; (2) in comparison, n_e is in the range of 10¹⁵–10^19.5 cm⁻³ (3σ) when kT ≈ 0.5 keV. The former fit is slightly preferred than the latter fit according to our MCMC analysis, although they have similar likelihood within the 2σ uncertainty range.

• Only an upper limit of kT is obtained due to the lack of data below 1 keV. Future observations, such as from NICER, will be able to help constraining the soft X-ray spectral shape in BH XRBs with more details.

To sum up, we can only obtain an upper limit of the contribution of the disc thermal component in the hard state of GRS 1716−249 due to the modest column density along the line of sight and the lack of data below 1 keV. The disc thermal component is also degenerate with the density parameter in the high-density disc reflection model used in this work. We note that when n_e = 10¹⁵ cm⁻³ is used, a disc blackbody component with kT ≈ 0.5 keV is required to model the soft X-ray emission. This result is consistent with the analyses in Bassi et al. (2019).