Reflection Modeling of the Black Hole Binary 4U 1630$-$47: the Disk Density and Returning Radiation

The models, (a) and (b), are defined explicitly in Table 3. Model (a) represents the standard low-density reflection model, and model (b) is the high-density reflection model. The component TBabs accounts for interstellar absorption, implementing the elemental abundance tables of Wilms et al. (2000). We adopt the Verner et al. (1996) atomic cross sections. The simplcut component is a model for inverse-Compton (IC) scattering in a coronal plasma (Steiner et al., 2017), and is an extension of the model simpl (Steiner et al., 2009b). simplcut acts as a convolution model, scattering the disk photons (as described by the multitemperature disk blackbody component diskbb; Mitsuda et al. 1984) that traverse hot coronal plasma. It uses the IC spectral calculations given by the Xspec model nthComp (Zdziarski et al., 1996; Życki et al., 1999). The temperature of the photon pool for scattering is tied to the diskbb value. The parameter $f_{\rm sc}$ sets the proportion of photons which are scattered in the corona. Parameterizing the fraction of photons scattered into the power-law component in this way simplifies the two-part problem: in reality some fraction of photons enter the coronal plasma (a covering fraction), and then some portion of those photons undergo scattering (the optical depth, $\tau$, of the corona sets this quantity). The $f_{\rm sc}$ parameter combines these two quantities, by assuming all photons enter the corona. Thus $f_{\rm sc}$ is a lower bound on the optical depth of the corona, where $\tau\geq-\ln(1-f_{\rm sc})$. From a modeling perspective, implementing simplcut as a model for IC scattering in a corona is preferable to treating the disk and corona components independently (with their own normalizations) because the former enforces self-consistency in the fitting procedure (the coronal IC component can only exist with some total number of disk photons, thus placing limits on the disk flux in the modeling procedure).

Models (a) and (b) are identical in regard to the disk and coronal spectral components, but they differ in the nature of the reflection spectrum. Model (a) assumes a disk density fixed at $10^{15}~{}{\rm cm^{-3}}$. The reflection spectrum is calculated assuming illumination atop the disk by an IC coronal spectrum, as given by the nthComp model (Zdziarski et al., 1996; Życki et al., 1999), and thus coincides with the illumination spectrum given by simplcut$\otimes$diskbb.

The model reflionxHD is the latest version of the long-standing reflection model reflionx (Ross & Fabian, 2005, 2007), which operates as a grid of reflection spectra which can be implemented in Xspec to perform statistical modeling of X-ray spectra. reflionxHD differs from the original reflionx model in two key ways: (i) the irradiating continuum ionizing the disk is based on the nthComp Comptonization model, as opposed to the approximation of a cutoff power law; (ii) the disk number density is a model parameter extending from $10^{15}~{}{\rm cm^{-3}}$ up to $10^{22}~{}{\rm cm^{-3}}$. The model component relconvlp is part of the relxill distribution of models (García et al., 2014; Dauser et al., 2014) for relativistic reflection, and acts as a convolution model to relativistically smear any input spectrum, assuming a lamppost geometry for the irradiating source. Thus the convolution procedure relconvlp$\otimes$reflionxHD applies relativistic smearing to the non-relativistic reflionxHD model. We adopt relconvlp, the lamppost version of the coronal geometry, based upon previous constraints on the dimensionless BH spin, $a_{\star}=0.985$, and a steep emissivity, $q>9$ (King et al., 2014).

In all fits to the PCA data shown in Figure 1 and Table 1 we treat the model parameters as follows. Based on explorative preliminary fits to all the data, we fix the hydrogen column density ($N_{\rm H}$, a parameter of the TBabs spectral component) to $1.4\times 10^{23}~{}{\rm cm^{-2}}$. This decision was based on small yet impactful differences between the constrained value of $N_{\rm H}$ as a function of spectral state and data counts, and an approximate range of values between $1.3$ and $1.5\times 10^{23}~{}{\rm cm^{-2}}$. We further justify this choice due to the uncertainties on $N_{\rm H}$ in the literature, as measured in the X-ray band (Tomsick et al., 1998; Tomsick & Kaaret, 2000; Gatuzz et al., 2019), as well as the presence of dust scattering along the line of sight (Kalemci et al., 2018). The simplcut model has three variable parameters: $\Gamma$, the power law photon index, $f_{\rm sc}$, the scattering fraction, and $kT_{\rm e}$, the coronal thermal electron temperature. The diskbb component is characterized by the disk normalization $N_{\rm disk}$, and the inner disk temperature $T_{\rm in}$. The reflionxHD model has parameters for the inner disk temperature $T_{\rm in}$ (this sets the input spectrum for IC scattering, see Jiang et al. 2020, submitted to ApJ), photon index $\Gamma$, and electron temperature $kT_{\rm e}$—these are all tied to their respective values as given by simplcut and diskbb. Further reflionxHD parameters are disk ionization $\xi$, iron abundance $A_{\rm Fe}$, density $n_{\rm e}$, and normalization $N_{\rm refl}$. We fix $A_{\rm Fe}=5$ in our fits, again based upon preliminary fits in which typical constrained values were found to be close to $5$, with fairly large uncertainties.

The model relconvlp has parameters which pertain to the source geometry: the black hole spin $a_{\star}$, inner disk radius $R_{\rm in}$, outer disk radius $R_{\rm out}$, disk inclination $i$, the lamppost height $h$, and $\Gamma$, the powerlaw index of the IC spectrum. We fix the black hole spin to maximal, $a=0.998$, and inclination to $i=64^{\circ}$, based on the reflection modeling results of King et al. (2014), which also made use of NuSTAR observations of the soft state of 4U 1630$-$47 (observation 40014009001). The inclination constraint also agrees with the lower and upper bounds provided by Tomsick et al. (1998) ($>60^{\circ}$) and Seifina et al. (2014) ($<70^{\circ}$) respectively. However, we note that there is some evidence that the inner disks of BHBs in general, as probed via reflection modeling, may be misaligned with the outer disk, and thus the orbit in general (see, e.g., Connors et al. 2019). In addition, we fix the inner disk radius to the ISCO ($R_{\rm in}=R_{\rm ISCO}$), again based upon preliminary fits which revealed its typical value to be very close to $R_{\rm ISCO}$, and poorly constrained. Thus, the reflection properties in models (a) and (b) vary only in the disk ionization $\xi$, the normalization $N_{\rm refl}$, lamppost height $h$, and the disk density $n_{\rm e}$, which varies freely in model (b)—up to $10^{22}~{}{\rm cm^{-3}}$—but is fixed at low density in model (a). With this carefully curated modeling approach we are able to capture the general differences between low and high density reflection modeling by inspecting the key parameter trends with maximal degrees of freedom in the spectral fits.

Finally, during preliminary fits with both Models (a and b) we also noticed that the PCA detects the well-known Fe XXVI absorption feature at $\sim 6.9$ keV, associated with an ionized wind, and previously detected in 4U 1630$-$47 NuSTAR spectra during the soft state (see, e.g., King et al. 2014; Hori et al. 2014). The feature is visible in spectra from the $2002\mbox{--}2004$ outburst at hardness ratios below $\sim 0.4$ (see Figure 1 for reference). As a result, we apply the gabs model to those datasets in which the feature is present, in order to factor out those residuals. We fix the centroid energy of the line to $6.9$ keV, the Gaussian width to $0.01$, and allow only the strength of the line to vary freely.

The results of our modeling of the PCA data are summarized visually in Figure 2, with numerical parameter values and their confidence limits presented in Tables 5 and 6 in Appendix A. The first and most fundamental result to highlight is that the disk density, $n_{\rm e}$, typically trends towards values $>10^{20}~{}{\rm cm^{-3}}$ across all spectral states. Secondly, there are some clear and significant differences in the properties of the corona when comparing low and high density reflection models. As the model reaches higher densities, lower values of $\Gamma$ and $kT_{\rm e}$ appear to be preferred, particularly during the more intermediate states—this is far clearer in our results for the $1998$ outburst data, whereas in our modeling of $2002\mbox{--}2004$ the systematic differences are much less obvious, which is likely because the model does not fit as well to those observations. Thus it appears that when higher BHB disk densities are assumed (which is arguably more appropriate based on simple BH mass-scaling arguments, see, e.g., García et al. 2016), one should systematically expect to derive a cooler corona with a harder IC continuum. We discuss the significance of this result in more detail in Section 4.

During the $1998$ outburst, in which the source was caught in a harder state, there are clear differences in the disk properties ($T_{\rm in}$ and $N_{\rm disk}$) between Model (a) and Model (b), as well as a weaker constraint on the coronal height, $h$, allowing a lamppost at almost $100~{}R_{\rm g}$. The caveat of such statements is that observation 30178-01-01-00 has a low number of counts, and as shown by the reduced $\chi^{2}$ panel at the bottom of Figure 2, the data are over-fit. Nonetheless, it is curious that the introduction of softer reflected emission due to high-density effects appears to be more consistent with a hotter disk with a smaller inner radius.

We note that Figure 2 does not show results for the low-density reflection modeling of observations 80117-01-03-00G and 80117-01-07-01 (at hardness ratios of 0.60 and 0.53 respectively). This is because the low-density model (Model a) struggles to fit the data within our parameter setup (with the inclination and inner radius frozen). Rather than show results of the improvements we can make to the modeling with additional free parameters, we choose to show the stricter comparisons across hardness.

Figure 2 also shows that generally the disk ionization, presented in log units ($\log[\xi/{\rm(erg~{}cm~{}s^{-1})}$), is $\sim 4$, but in the much softer states reached during the $2002\mbox{--}2004$ outburst it drops two orders of magnitude to $<2$. This significant decline may not be meaningful because the reflection model dominates the fit over the IC component—thus the constrained value of $\log\xi$ is less reliable, despite the trend coinciding with the drop in coronal flux as expected. In addition, we should expect the bright disk emission to be contributing to the disk ionization, which makes the drop in flux even more suspect. This issue was highlighted in full by Connors et al. (2020), and is possibly evidence that a IC continuum is an inappropriate description of the irradiating continuum in BHB soft states (due to the superior disk flux with respect to the IC emission). We address this in modeling of the NuSTAR and Swift-XRT data in Section 3.2.

There are not many clear differences in the parameter trends between fits to the $1998$ and $2002\mbox{--}2004$ outbursts. However, one can notice higher disk temperatures during the $2002\mbox{--}2004$ outburst, particularly as the source reached a softer spectral state, where we would expect a higher temperature disk. Constraints on quantities such as the coronal electron temperature, $kT_{\rm e}$, scattering fraction, $f_{\rm sc}$, and coronal height, $h$, are too weak to properly distinguish the differences between outbursts. The electron temperature $kT_{\rm e}$ is best constrained in the intermediate states, where we find good agreement between outbursts, despite the doubling in count rate during the latter outburst in $2002\mbox{--}2004$ (though as Tables 5 and 6 in Appendix A show, there are slight numerical differences, with the higher luminosity observations appearing to show a cooler coronal gas).

The two reflection models we use are relconv$\otimes$reflionxHD and relxillNS. The model relxillNS was implemented in Connors et al. (2020), and will be described in detail in a forthcoming paper (Garcia et al., in preparation). It is a variant of the relxill consortium of relativistic reflection models (García et al., 2014; Dauser et al., 2014), and behaves similarly to relxill, the difference being that the irradiating continuum is a single-temperature blackbody spectrum characterized by $kT_{\rm refl}$, as opposed to a cuttoff power law or IC spectrum.

We use relxillNS purely as a consequence of empirical results from modeling, showing that power-law-like reflection proves inadequate during the soft spectral state (Connors et al., 2020). The fundamental driver of this disagreement is the strength of the Compton hump in the reflection component, which is prominent in reflection models which assume an IC irradiating continuum. Reflection from an irradiating blackbody spectrum, on the other hand, naturally results in a weaker Compton hump due to the softer irradiating continuum. The relxillNS component acts as a proxy for the reflection spectrum resulting from returning disk radiation illuminating the disk. This is not a model of returning radiation, and the blackbody spectrum we use is a single-temperature one, as opposed to a disk spectrum, so there are strong caveats to using this model. However, it provides us with a way to compare the expected outcome when the disk is reflecting a blackbody-like spectrum as opposed to a power-law-like spectrum.

This model comparison is suitable given how disk-dominated the source was during the NuSTAR/Swift observations, as shown by simulated data presented in Figure 3, alongside the superior energy resolution of NuSTAR (a feature which allows for stricter reflection model comparisons), and the soft X-ray coverage provided by Swift-XRT. The simulated PCA observations are based on model fits to the NuSTAR spectra where the assumed model is an absorbed Comptonized disk with an absorption line from an ionized wind: TBabs(simplcut*diskbb)*gabs. One can see that all three NuSTAR/Swift observations occurred during the soft state, and basic fits reveal typical disk temperatures of $\sim 1.4$ keV.

We described the data reduction and grouping in Section 2.2. Here we describe the details of the model treatment and fitting procedure. The first model is similar to Model (b) as presented in Section 3.1, instead now we require additional calibration corrections to fit the NuSTAR and Swift-XRT spectra simultaneously, as well as the inclusion of dust scattering effects in the model (Kalemci et al., 2018).

In addition, all the NuSTAR spectra clearly display a strong $\sim 6.9$ keV absorption line, commensurate with a previously confirmed ionized disk wind (Díaz Trigo et al., 2014; King et al., 2014). Given the increased spectral resolution of NuSTAR with respect to RXTE, we replaced our simplistic Gaussian model for the absorption line with the full, self-consistent xstar photoionization model (Kallman & Bautista, 2001). We generate a grid of models assuming solar metallicity, and an input blackbody spectrum at a temperature of $1.4$ keV, based on initial absorbed disk$+$powerlaw fits to the simultaneous NuSTAR and Swift-XRT observations, showing a disk inner temperature of $1.4$ keV. We fix the gas density to $n=10^{15}~{}{\rm cm^{-3}}$, based on typical X-ray wind densities (Miller et al., 2008). We chose an ionizing luminosity of $10^{38}~{}{\rm erg~{}s^{-1}}$, and a turbulent velocity of $1000~{}{\rm km~{}s^{-1}}$, based upon the variation in its best fit value identified by Díaz Trigo et al. (2014). Typical wind velocities can fall in the range of $300\mbox{--}3000~{}{\rm km~{}s^{-1}}$, based upon rough estimates of the Keplerian velocity range corresponding to the radial launching window allowed for thermally driven winds (Ponti et al., 2012). Our grid then has just three variable parameters, the column density, $N_{\rm H,wind}$, and the ionization, $\log\xi_{\rm wind}$, and the outflow velocity of the wind, $v_{\rm wind}$.

The models, now (b.2) and (c), are outlined in Table 3. We replaced the lamppost version of the relativistic convolution model, relconvlp, with relconv, which instead parameterizes the emissivity profile of irradiation from $R_{\rm in}$ to $R_{\rm out}$ via the index $q$, where $\epsilon(r)\propto r^{-q}$. This is based on the lack of a lamppost version of the model relxillNS, and it therefore allows a more direct comparison of models (b.2) and (c).

The model crabcorr (Steiner et al., 2010) applies a correction to the detector response of a given instrument to match the instrument-calibrated normalization and power-law slope of the Crab spectrum (Toor & Seward, 1974). Thus its model parameters are $N_{\rm CC}$ and $\Delta\Gamma_{\rm CC}$, which renormalize and apply a shift in power law slope respectively. We fix $N_{\rm CC}=1$ and $\Delta\Gamma_{\rm CC}=0$ in FPM A, and allow only $N_{\rm CC}$ to vary for FPM B. In accordance with the results of Steiner et al. (2010), we apply a $\Delta\Gamma_{\rm CC}=-0.04$ to the Swift-XRT spectra (which are contemperaneous with NuSTAR observations 40014009001 and 90002004004). $N_{\rm CC}$ is then left free to vary as a calibration constant for the Swift-XRT data relative to NuSTAR.

The model xscat (Smith et al., 2016) accounts for dust scattering of source photons. For sources with a large hydrogen column density ($N_{\rm H}$) along the line of sight it is common for a dust scattering halo to form—this is because photons at soft X-ray energies mostly undergo small-angle scattering with dust grains along the line of sight. The effects of dust scattering were explored fully by Kalemci et al. (2018) in the case of 4U 1630$-$47. They found that the bulk of the dust scattering is likely due to a molecular cloud located close to the source ($\sim 0.94D$, where $D$ is the distance to the source). The xscat model implementation in Xspec has a few parameters, the most important being $X_{\rm pos}$, the position of the dust scattering halo along the line of sight to the source, $N_{\rm H}$, the hydrogen column density of the scattering source (units of $10^{22}~{}{\rm cm^{-2}}$), and $R_{\rm ext}$, the radius of the circular extraction region in arcseconds. In all our fits we fix $X_{\rm pos}=0.9$, $R_{\rm ext}=80$” and $R_{\rm ext}=50"$ for the NuSTAR and Swift-XRT spectra respectively. The value of $0.9$ for the position of the scatterer is based on the fact that it is likely that multiple scattering regions exist along the line of sight, despite there being a dominant one located at $X_{\rm pos}=0.94$ (Kalemci et al., 2018). We then tie $N_{\rm H}$ to the line of sight hydrogen column given by the TBabs model. This modeling approach limits the different fits to the same level of dust scattering, as it relates to the interstellar gas absorption ($N_{\rm H}$).

We perform the spectral fits to all three NuSTAR observations (along with the simultaneous Swift-XRT observations) jointly, tying model parameters we have a priori determined should remain constant. We describe the full parameter setup in the following subsection.

Given the superior spectral energy resolution of NuSTAR over the PCA, we allow more freedom in the models. We let the iron abundance, $A_{\rm Fe}$, disk inclination, $i$, BH spin, $a_{\star}$, and absorption column density, $N_{\rm H}$, vary freely, tying them across the three observations such they are assumed to be constant over time. In addition, the disk density, $n_{\rm e}$, is allowed to vary in model (b.2), again with its value tied across the three observations, given all three are at similar count rates, and disk density is dependent on accretion rate.

As shown by Figure 3, all three observations are during a soft state, and so we assume the thin accretion disk extends to the ISCO ($R_{\rm in}=R_{\rm ISCO}$). We let $q$ vary freely in fits with Model b.2, based upon the low-$h$ constraints found in our fits to the PCA data (see Section 3.1). We fix $q=3$ in our Model c fits, based on preliminary predictions regarding the shallow emissivity profile of returning disk radiation (Wilkins et al., 2020). Whilst relxillNS is not itself a physical model of returning radiation—which should be noted as a significant caveat in our spectral fitting approach—the work of Wilkins et al. (2020) shows that the emissivity profile should be flatter than that of a compact corona irradiating the disk. In Model c we also separate the temperature of the irradiating blackbody component from the intrinsic inner disk temperature in diskbb: $kT_{\rm refl}\neq kT_{\rm in}$. This is based on the inherent uncertainty regarding the relativistic energy shifts experienced by the returning disk photons. We also fix the disk density in relxillNS (Model c) to the highest available value the reflection table model affords us, $10^{19}~{}{\rm cm^{-3}}$.

We let the xstar parameters vary freely for each individual dataset. These are the wind column density, $N_{\rm H,wind}$, gas ionization, $\log\xi_{\rm wind}$, and wind velocity, $v_{\rm wind}$. The column density parameter has limits in our generated xstar table given by $5\times 10^{21}\mbox{--}5\times 10^{23}$, and the ionization has limits given by $2\mbox{--}5$ in log units. The wind velocity is determined by a freely varying redshift parameter on the xstar table model.

The resulting spectral fits are shown in Figure 4, and corresponding parameters with 90% confidence limits are shown in Table 4. The key result to highlight is that when modeling all three observations Model (c) performs slightly better than Model (b.2). The superiority of reflected returning radiation to fit the broadband X-ray spectra is clear both from the lower total $\chi^{2}$ for one additional degree of freedom, and the patterns in residuals shown in Figure 4—Model b.2 underfits the high energies and introduces curvature around $10\mbox{--}20$ keV. Thus our results make a strong case for returning radiation as a contributor to the observed reflection spectrum in BHB soft states (Connors et al., 2020).

The hydrogen column density is consistently within the range of $9.5\mbox{--}9.7\times 10^{22}~{}{\rm cm^{-2}}$, in agreement with recent high-resolution spectroscopic observations with Chandra (Gatuzz et al., 2019). The higher value ($14\times 10^{22}~{}{\rm cm^{-3}}$) assumed in our fits to the PCA data (see Section 3.1) may reflect the increased effects of dust scattering due to the large aperture size of the PCA detectors.

We find that the disk density (when applying Model b.2) exceeds $9\times 10^{21}~{}{\rm cm^{-3}}$, which is consistent with our findings from modeling of PCA data shown in Section 3.1. However, inspection of Figure 4 also shows that given the source is in the soft state, reflection from illumination of the disk by the corona (relconv$\otimes$reflionxHD) tends to over-fit the power law component in the data, introducing an undesirable curvature in the residuals. Indeed the parameter constraints in Table 4 shows this quantitatively in the systematic difference in coronal parameters between Models (b.2) and (c): Model c retrieves similar values of $f_{\rm sc}$, but significantly lower values of $\Gamma$, or in other words a much harder power law component.

The BH spin is high in both Model b.2 and Model c fits, with $a_{\star}=0.989^{+0.001}_{-0.002}$ and $a_{\star}=0.85\pm 0.07$ respectively. The Model b.2 constraint on BH spin thus agrees fairly well with those found by King et al. (2014) in their analysis of the 2013 observations of 4U 1630$-$47, but the Model c constraint is considerably lower, albeit remaining a high BH spin. This is encouraging given that the models differ slightly, with King et al. (2014) adopting the model refbhb (Ross & Fabian, 2007) which includes the transmitted emission from internal disk dissipation.

The Model b.2 inclination constraint agrees well with binary orbit inclination limits ($60^{\circ}<i<70^{\circ}$; Tomsick et al. 1998; Seifina et al. 2014), as well as with modeling by King et al. (2014). However, we find a systematic difference in the inclination obtained from Model c (returning radiation) fits, with $i=37^{+1}_{-2}~{}{\rm deg}$. We note that relxillNS is a very preliminary model lacking physical consistency, and so we cannot derive much meaning from this difference in inclination. One can speculate that a contrast in inclination with respect to power-law reflection modeling may be expected since there is additional curvature around the Fe line region, as well as a weaker Compton hump—once the reflection line and continuum adjust significantly, one may expect the inclination to re-adjust to fit the data. We refrain from going beyond this basic explanation, in an effort to avoid speculation on the physical interpretation.

Our ionizing wind parameter constraints reveal a similar dichotomy outlined by King et al. (2014), between ionization state and blueshift of Fe absorption lines. We find that some of our Model/observation combinations (see Table 4) prefer wind velocities on the order of $\sim 2000~{}{\rm km~{}s^{-1}}$ with higher wind ionization levels, $\log\xi_{\rm wind}/[{\rm erg~{}cm~{}s^{-1}}]\sim 4$. Others prefer an ultra-fast outflow with $v_{\rm wind}>10,000~{}{\rm km~{}s^{-1}}$, and ionizations closer to $\log\xi_{\rm wind}/[{\rm erg~{}cm~{}s^{-1}}]\sim 3$. Whilst the absolute values differ from the results of King et al. (2014), the trends are similar. However, high-resolution spectroscopy provided the best constraints on the disk wind, revealing a relatively low velocity, high ionization, and high column density (Kubota et al., 2007; Díaz Trigo et al., 2014). This is indeed what we find for four of the six model fits. There is no clear and obvious distinction between Models b.2 and c in terms of this dichotomy in the wind properties, so we suggest that the dichotomy is predominantly due to the limiting resolution of the NuSTAR spectra compared to, for example, Chandra high energy transmission grating (HETG) spectra and XMM-Newton spectra. This explains why we King et al. (2014) derived a similar dichotomy.