© The Authors 2024

1. Introduction

Reflection spectra (RS) modeling in X-ray binaries (XRBs) or active galactic nuclei (AGNs) offers a promising approach to studying the space-time geometry around black holes (BHs) (Miller 2007; Fabian & Ross 2010; Reynolds 2014; Bambi et al. 2021). To model the RS, it is typical to simulate the RS theoretically for a specific incident radiation illuminating a slab. Changes in incident radiation and/or the properties of the slab result in different final spectra, which are stored as tables known as reflection tables, such as reflionx (Ross & Fabian 2005; Tomsick et al. 2018) and XILLVER (García & Kallman 2010; García et al. 2013, 2016). Very recently, similar reflection tables have been built including polarization as well (Podgorný et al. 2022). Theoretically simulated RS used to be convolved relativistically and integrated over the disk to fit the observed RS from an accretion disk (Fabian et al. 1989; Laor 1991; Dabrowski et al. 1997; Dovčiak et al. 2004; Dauser et al. 2010; García et al. 2014). During the early development of reflection tables, the density (n) of the slab was fixed at 10¹⁵ cm⁻³, which is the typical density for the accretion disk of AGNs. However, for XRBs, we expect that the density of the accretion disk, especially in the inner region, can be much higher than 10¹⁵ cm⁻³ (Frank et al. 2002; Kato et al. 2008). Nonetheless, observations of XRBs are also typically fitted using reflection tables of 10¹⁵ cm⁻³ with some arbitrary normalization factor to match the observed flux. The reasoning behind doing this is twofold. First, making reflection tables for high-density slabs is much more computationally and physically demanding (Kallman et al. 2021). With respect to RS, we are primarily concerned with the shape of the iron line, which is very much dependent on how much the slab is ionized, namely, the ionization (ξ) parameter of the reflection tables (Ross & Fabian 2005, 2007; García & Kallman 2010; García et al. 2013; Tripathi et al. 2020). Arbitrary normalization also hides the uncertainties in the distance to the source and its inclination. Thus, we can uncover the physical picture of the source to a great extent by focusing only on the shape of the observed spectra. However, it was recently reported that density can play a crucial role in altering the shape of the iron line, the Compton hump, and can provide thermalized blackbody emission in soft X-ray that is due to the reprocessing of the illumination (García et al. 2016; Tomsick et al. 2018). The higher density disk in XRBs (compared to AGNs) makes the fitting of XRB data with high-density tables more reasonable.

As an essential consequence of the reflection, a fraction of incident radiation is always absorbed and emitted in lower energies. The albedo gives the fraction of incident radiation that is reflected back. From early studies, we know that the albedo of an ionized disk can vary in the range ∼0.3–0.7 (Basko et al. 1974; Nayakshin & Kallman 2001; Psaltis 2002) depending on ionization, which means a significant fraction of the incident radiation will always be absorbed and emitted as blackbody and fluorescent lines. With the increment of the density of the illuminated slab, if we need the same ξ to fit the RS, more illumination is required; this subsequently increases the thermalized blackbody radiation (i.e., the temperature at which the blackbody peak will appear). This is the reason for the appearance of thermalized blackbody radiation in the soft X-ray at higher slab densities, which is shown clearly in the left panel of Fig. 1. Due to the high dependence of the absorbed flux on illumination, blackbody radiation will also appear in soft X-rays for very high ionization cases, even at low density, as shown in the right panel of Fig. 1. If the observed RS demands sufficiently high ionization and density of the disk, then consequently the reprocessed blackbody will appear in soft X-ray. The effect of galactic absorption also sharply decreases from ultraviolet to X-ray. This makes it possible to observe the reprocessed blackbody component in soft X-ray if it is produced in the local frame of the source. However, in low/hard states of the XRBs, we do see RS but do not see any such blackbody radiation (Remillard & McClintock 2006; Dunn et al. 2010). This is the motivation behind our investigation of the consistency of the disk reflection of highly irradiating flux observed as a power law specifically in low/hard states of XRBs. One general shortcoming of using the XILLVER or reflionx tables is that none of them include intrinsic viscous dissipation. While this generates radiation in ultraviolet for AGNs and can be ignored for X-ray spectra, in the case of XRBs, this will certainly be significant (especially in high/soft state). The absence of any blackbody component in the low/hard state of the XRBs supports the usage of reflionx or XILLVER, which we present here. Although this issue has been addressed in refhiden reflection table (Ross & Fabian 2007), the lower limit of the temperature of the blackbody due to viscous dissipation is 0.4 keV, making it suitable only in bright high/soft state.

Fig. 1. Variation in EFE with density (n) and ionization (ξ) from reflionx_hd table with appropriate normalization. The corresponding n and ξ are written in the legend. The reflionx_hd table has three input parameters: n, ξ, and Γ. For a variation of either n or ξ, the fixed values of other input parameters are noted on the plot. More on the normalization of these spectra in footnote 9.

A significant amount of work has been undertaken very recently about fitting XRB RS with high-density tables (Tomsick et al. 2018; Jiang et al. 2019b, 2022; Chakraborty et al. 2021, even in the hard state, Liu et al. 2023, for AGNs see Jiang et al. 2019a; Mallick et al. 2022). However, even with high density tables, the arbitrary normalization factor is always used to scale the flux as per the requirement. For reflection tables, ξ and n defines completely the illuminating flux (F_inc). Again, RS is the radiation emitted by the slab when it is in radiative equilibrium with the incident radiation; this makes the reflected flux (F_refl) equal to F_inc (García & Kallman 2010) and makes it possible to define the reflected flux in the following way:

$$\begin{array}{c}\hfill F_{\mathrm{refl}}=F_{\mathrm{inc}}=\frac{\xi n}{4\pi }.\end{array}$$(1)

Here, F_refl indicates the radiation emitted by the disk per unit area per unit time in its local frame. Unless stated otherwise, we always consider the full energy range 1 eV to 1 MeV to compute F_refl, which also includes the reprocessed blackbody component. Therefore, once ξ and n are fixed, we know F_refl. From the fitting of RS, not only do we get ξ and n, but also the inner (r_in) and outer (r_out) radii of the disk (in units of gravitational radius), giving the area emitting F_refl. For some XRBs, we know the distance, inclination, and mass of the BH quite well or at least the narrow range in which they lie (Tetarenko et al. 2016). Using the mass of the BH, we can find the gravitational radius and consequently the area of the disk in physical units. Using F_refl and total emitting area of the disk we can calculate the total luminosity emitted by the disk in its local frame as reflection. As the disk emits as a cosine source, from the distance and inclination, we can then estimate observed reflected flux (F_refl, o) from the computed luminosity (Equation (A.5)). These simple steps provide exact value of observed reflected flux from the fitted RS and takes away the requirement of arbitrary normalization. Therefore, if we know the mass of the BH, along with the distance and inclination of the source, it should be possible to fit the observed RS, considering both the shape and the flux. This may put a strong constraint on reflection in the low/hard state of XRBs due to reprocessed blackbody.

Zdziarski & De Marco (2020) (hereafter ZDM20) illustrated the concern over the absence of thermalized blackbody radiation in the hard state of XRBs by linking the primary radiation (the power-law) with the disk area. If the corona is on the axis and not very high from the disk, the illuminating flux at the innermost radius of the disk inversely depends (approximately) on the square of the radius (Equation (3) in ZDM20). If the disk is extended to the innermost stable circular orbit (ISCO), the illuminating flux at the innermost radius will be larger than in truncated disk geometry, creating more significant thermalized blackbody radiation due to reprocessing. For typical values of the albedo, the peak for thermalized blackbody should appear in the range 0.1–1.0 keV for the observed primary power law in the hard state of XRBs, which we do not observe (Fig. 3 of ZDM20, Zdziarski et al. 2021a). In summary, the disk in the hard state of XRBs cannot approach too closely to the BH due to the restriction from thermalized blackbody as a consequence of illumination of the disk by primary hard X-ray. However, the above study requires the reflected flux to be connected with the primary power-law, which consequently includes the uncertainty in albedo as well as the geometry of the source. The amount of illumination (and thus also reflection) on the disk depends on the illuminating solid angle, namely, the geometry of the power-law emission region and its position with respect to the reflecting disk. Usually, the normalization of the reflection component is not linked with the normalization of primary power-law while fitting the observation due to geometrical uncertainty. The ratio between the two is typically denoted as reflection fraction (Dauser et al. 2016). Therefore, in this study, we also kept the normalization of the primary power-law free, independent of reflection component. However, as mentioned in the earlier paragraph, for RS, we made sure the fitting was consistent in terms of the shape and flux simultaneously, depending on the area of the disk, distance, and inclination to the source; with ξ and n define the amount of illumination and, thus, the reflection as well. With this setup, we can verify whether disk extending to ISCO or truncated geometry is better for high-density disk reflection incorporating self-consistent normalization of the RS simultaneously from its shape and flux.

Another pivotal issue is the emissivity profile, which requires attention in the context of this investigation. It is natural that due to a change in the illumination of the disk with radius, the reflection by the disk also changes. The change in emission of the disk with radii is represented through the emissivity profile:

$$\begin{array}{c}\hfill I_{\nu }(r)\propto r^{-q},\end{array}$$(2)

where r is in units of R_g ($R_{\mathrm{g}}=\frac{\mathit{GM}}{c^2}$, gravitational radius) and q is the emissivity index (Wilkins & Fabian 2012; Dauser et al. 2013). Also, q plays a crucial role in fitting the broad Fe emission line present in the observed RS (Svoboda et al. 2012; Fabian et al. 2014). In lamp-post geometry, for large values of r for the disk, compared to the height of the corona, the q value will be 3. For the Shakura-Sunyaev disk, the q value is also 3 (applicable for a radius much larger than zero torque inner boundary of the disk, Equation (3.52) in Kato et al. 2008). This is relevant if the corona lies just above the disk, does not extend very close to the BH and is fed by the accretion process. We fix q = 3.0 as done in ZDM20, Zdziarski et al. (2021b), Marino et al. (2021) etc. Later in this work, we also assessed the robustness of our final conclusion on the assumed q value.

A discussion of the source and observation is given in Section 2. We describe how we modeled the RS in Section 3. First, we reproduce the expected inconsistency of reprocessed blackbody with observation in Section 3.1. In Section 3.2 we build fully self-consistent model and report its results by simultaneously fitting the shape and flux of the RS including reprocessed blackbody. We discuss a few crucial points including possible geometry of the source in Section 4 and present our conclusions in Section 5. Another alternative model for fitting the shape and flux simultaneously is presented in Appendix A. For completeness, we present the fitting results of NuSTAR and NICER observations in Appendix B.

2. Observations and source details

We used NuSTAR observation of MAXI J1820+070 with ID 90401309002 during its low/hard state (Epoch 1 in Zdziarski et al. 2021b; Marino et al. 2021). We reduced the observation using NuSTARDAS following the standard procedure¹ with heasoft version 6.33.1. The reasoning behind the choice of this observation is that our study is relevant in the low/hard state of XRBs when the thermal component is the weakest. In addition, following the detailed studies of this source (Torres et al. 2019; Atri et al. 2020; Torres et al. 2020; Mikołajewska et al. 2022), we fixed the BH mass to 8 M_⊙, distance to 3 kpc, inclination to 66°, and hydrogen column density for absorption in inter stellar medium (N_H of the TBabs model) to 1.4 × 10²¹ cm⁻² (Kajava et al. 2019; Zdziarski et al. 2021b, 1.5 × 10²¹ cm⁻² is used in Chakraborty et al. 2020).

Although we used NuSTAR FPMA data within the energy range 3.0–78.0 keV for the fitting procedure, the effect of thermalization can be more prominent in the soft X-ray region below 3 keV. Therefore, we considered the quasi simultaneous observation by NICER with ID 1200120104 for energy range 0.5–10.0 keV. We reduced the NICER data using NICERDAS software version 2024-02-09_V012 in heasoft 6.33.1, following the standard procedure. We used the latest NICER calibration files of CALDB release “xti20240206” to extract level 2² data and, finally, level 3 spectral data, rmf, and arf files³ (which include a 1.5% spectral systematic error in the energy range 0.5–10.0 keV). However, the low-energy data from NICER typically reveal some calibration features and require some phenomenological models (e.g., gabs, edge) for the fitting procedure (Wang et al. 2020, 2021; Ludlam et al. 2021, 2022; Kumar et al. 2022; Feng et al. 2022; Moutard et al. 2023; Svoboda et al. 2024b,a). In addition, we have to take into account the constant factors for cross-normalization across different instruments to fit simultaneous data from NuSTAR and NICER. Therefore, to keep the fitting procedure simple and to concentrate mainly on the reprocessed blackbody due to reflection, we fit the NuSTAR FPMA data only in the energy range of 3.0–78.0 keV; we show the NICER data in the plot for energy range of 0.5–10.0 keV. This efficiently demonstrates whether the fitted model for NuSTAR data is overestimating the reprocessed blackbody (or not) with respect to the simultaneous low-energy data from NICER. If the model flux remains below the NICER data, then we can safely say that reprocessed blackbody due to illumination of the disk does not contradict the non-observation of any blackbody component in the hard state. In Appendix B, we present the complete fitting of NICER, NuSTAR FPMA and NuSTAR FPMB data with our fully self-consistent model (details in Section 3.2.1) which fits the shape and flux of the RS simultaneously including reprocessed blackbody. We do not focus on the spin measurement of the BH here, instead freezing it to 0.998 for the entire study.

3. Modeling the RS

From the existing literature, it is rather clear that one reflection component is insufficient to fit the observed RS for most of the sources (Wang-Ji et al. 2018; Sridhar et al. 2020; Liu et al. 2023). For the observations of MAXI J1820+070 as well, the inclusion of two reflection components (a broad component for relativistic reflection and a narrow component for distant reflection) improves the spectral fits significantly (Buisson et al. 2019; Chakraborty et al. 2020; Zdziarski et al. 2021b; Marino et al. 2021). Therefore, we also include two reflection components to fit the RS. To model the interstellar absorption and primary power law, we used TBabs and cutoffpl, respectively, in XSPEC⁴ (Arnaud 1996). Therefore, our complete model for fitting becomes TBabs×(cutoffpl + reflection + reflection). To model the reflection, we have two choices available, which are mostly used by the community: reflionx and XILLVER. However, publicly available version of XILLVER tables⁵ only cover the energy range from 0.1 keV–1 MeV; whereas our work (specifically model C in Appendix A) incorporates the reprocessed blackbody demands estimation of flux in broader energy range of 1 eV–1 MeV. The flux estimation in this broad energy range also plays a crucial role in improving our understanding of the geometry of the source (see Section 4.2). For this reason, we used reflionx_hd table⁶ to model the high-density disk reflection. This table is prepared by illuminating the slab with a hard X-ray modeled by cutoffpl and fixing the cutoff energy at 300 keV (Tomsick et al. 2018; Ross & Fabian 2005). Here we fit the primary hard X-ray with cutoffpl setting the cutoff energy at 300 keV, and the photon index of the primary is linked to the photon index of the reflionx_hd table, which ensures the consistency between primary and reflection. Also, the reflionx_hd table assumes the iron abundance to be solar. The final models differ, depending on how the reflection components are modeled.

At first, to reproduce the expected discrepancy of observation with reprocessed blackbody, we fit the shape of the reflection component, assuming the disk density to be 10¹⁵ cm⁻³, and estimated the reprocessed blackbody from the reflection (see Section 3.1). Next, we built a model with high-density disk reflection that includes reprocessed blackbody in the fitting and self-consistent methodology to fit the observation simultaneously by its shape and flux, as described in Section 3.2.

3.1. Inconsistency with the reprocessed blackbody

To fit the reflection component, it is typical to convolve (using kyconv: Dovčiak et al. 2004, or relconv: Dauser et al. 2010; García et al. 2014) the reflection spectrum of the reflionx_hd table to include the general relativistic effects and the contribution from the whole disk. To understand the problem of reprocessed blackbody better and to reproduce it, we first followed ZDM20 to a great extent (but not fully) to fit the observed spectra by treating the reprocessed blackbody separately from the reflection component. This is achieved by fixing the density of the slab for reflection to be 10¹⁵ cm⁻³, making our first model, model A = TBabs×(cutoffpl + relconv⊗reflionx_hd (10¹⁵ cm⁻³)+ relconv⊗reflionx_hd (10¹⁵ cm⁻³)).

Here we kept the normalization parameter of the reflionx_hd table in XSPEC free while fitting. The best fit with this model (χ²/d.o.f. = 2157/1867) gives a reasonable estimate of the reflection. To estimate the unabsorbed reflected flux, we computed the flux corresponding to reflection component within the energy range of 0.1–1000.0 keV (different from 1 eV–1 MeV range used in rest of the paper) to ignore any possible contribution from thermalization (Dovčiak et al. 2022). We assumed that there is no reprocessed blackbody radiation above 0.1 keV for density of 10¹⁵ cm⁻³ (same as in ZDM20). Although the boundary of energy between reflection and reprocessed blackbody is vague, this is the best we can do if we treat those separately. We note that the RS in reflionx_hd table includes the reprocessed blackbody radiation within full energy range, but by fixing the density to 10¹⁵ cm⁻³ and lower boundary of energy to 0.1 keV, we may be able to ignore the reprocessed blackbody from the reflection. Assuming an average value of albedo of 0.5, we can consider that the flux emitted as reflection should give rise to similar blackbody flux, which we add as a separate component using the model diskbb (Mitsuda et al. 1984). From the fitted inner radius of the disk (r_in), and using distance (D) and inclination (θ) of the source, we can easily fix the normalization (A_dbb) of diskbb component, given by:

$$\begin{array}{c}\hfill A_{\mathrm{dbb}}={\left(\frac{r_{\mathrm{in}}/\mathrm{km}}{D/(10\mathrm{kpc})}\right)}^{2}cos\theta .\end{array}$$(3)

We can then find the temperature that equates the diskbb flux to reflected flux (Equation (10) in Zdziarski et al. 2021a). We note that here r_in/km indicates the inner radius of the disk in units of km. Flux in 0.1 keV–1000 keV energy range for the inner and outer reflection component becomes 1.0 × 10⁻⁸ and 8.1 × 10⁻⁹ ergs cm⁻² s⁻¹, respectively. Also, the fitting gives r_in to be 7.5R_g and 167.5R_g which set the T_in of the added diskbb to 0.35 keV and 0.08 keV, respectively, for the inner and outer reflection component (using D = 3 kpc, θ = 66°). The crucial difference between the above procedure and what followed in ZDM20 is that here we equate the blackbody flux with reflected flux instead of connecting it with the primary⁷.

Figure 2 shows that the reprocessed blackbody (diskbb component) due to inner reflection overshoots the data in the soft X-ray region significantly. This raises concerns because the reprocessed blackbody component must arise as a consequence of reflection, although we never observe any blackbody component along with the reflection component in the hard state of XRBs (as pointed out in ZDM20). For the outer reflection component, a high value for r_in makes the corresponding diskbb normalization large. Therefore, to achieve similar flux, T_in for the outer reflection becomes much lower compared to the inner reflection. The above exercise shows that despite the inclusion of two reflection components instead of one, and even estimating the reprocessed blackbody from reflected flux instead of connecting it with primary radiation, the inner reflection still leads to strong reprocessed soft X-ray emission, which is not present in the observed data. However, two crucial points are overlooked here. We added the reprocessed blackbody components on top of the already fitted model, which makes the result not robust. Ideally, the fitting should be done by including the reprocessed blackbody components because multiple components can eventually be adjusted among themselves and ultimately result in a consistent fit. Treating the reprocessed blackbody separately from reflection is another tricky undertaking, due to the ambiguous boundary between the two, which is very sensitive to the density and ionization of the slab (Fig. 1). Therefore, to draw any inference about the thermalization, we chose more accurate models in the next section (Section 3.2) to fit the observations.

Fig. 2. Best-fit model A for NuSTAR data. Expected reprocessed blackbody as a consequence of reflection are modeled as diskbb components, which are added separately over the best-fitted model A. diskbb(inner) and diskbb(outer) denote the diskbb components corresponding to thermalization due to inner and outer reflection, respectively. Simultaneous NICER data were presented to compare the fitted model in the lower energy. We present the unfolded spectrum to show the data of different instruments on the same scale.

3.2. Fitting the shape and flux of RS simultaneously including reprocessed blackbody

As explained in Section 1 (and later in more detail), from the fitting of shape of the RS emitted by the accretion disk, we can get all the required information to find the observed reflection flux. So, we should be able to connect the shape of the RS with its flux (i.e., photon counts) and avoid the arbitrary normalization in fitting. On the other hand, in the full energy range of 1 eV–1 MeV, RS from reflionx_hd includes both the reflection (i.e., Compton hump, emission lines) and the reprocessed blackbody. Therefore, the best model to infer the reprocessed blackbody would fit the shape and flux of the RS together, including the reprocessed blackbody. This is achieved completely by incorporating reflionx_hd within the kynrefionx model, including all the general relativistic effects (as detailed in Section 3.2.1).

3.2.1. Modeling RS with kynrefionx

Dovčiak et al. (2004) developed a fitting scheme within XSPEC to study the X-ray spectra in strong gravity regime. It gives us flexibility to compute the spectra emitted by an accretion disk because of the opportunity to include the local emission of the disk separately within the code. Once the local emission (emission by the disk in its co-rotating local frame) at different radii of the disk are defined, the KY code integrates it over the disk incorporating the general relativistic effects on radiation in Kerr space-time. There are different versions⁸ of the code to model different features of the accretion disk as outlined in Table 2 of Dovčiak et al. (2004). One such version is kynrefionx, which specifically can read reflionx tables and fit the observed RS, while including or excluding the primary. We refer to Dovciak (2004) for more details.

We modeled the RS independent of primary power-law and fix the emissivity index to 3.0, as done for model A in the previous section. We incorporated the reflionx_hd table as the local emissivity of the disk within kynrefionx. We assumed that the density of the disk remains same with radius and only ξ (and consequently flux) changes with radius depending on the value of the emissivity index. The code basically takes spectra corresponding to different ξ from the reflection table for different radii of the disk and integrates them over the disk area, including all general relativistic effects to give the final observed spectra. Due to the inclusion of RS of different ionization at different radii of the disk, this model is more accurate compared to model A, where the convolution model is applied externally; thus, it takes only one spectrum to convolve corresponding to one ξ value for all radii. Another crucial step forward is the exact normalization of the spectra. In kynrefionx, local emission of disk for reflection at a specific radius is taken from the reflionx_hd table, depending on ξ, n at that radius and normalized such that the total reflected flux in the energy range of 1 eV–1 MeV becomes equal to (ξn/4π), following Equation (1)⁹. As the local emission from each radius is normalized appropriately, the kynrefionx code integrates these local reflected spectra throughout the whole disk from r_in to r_out, which basically gives the reflected luminosity as a function of energy. To find the RS (and flux) in observer’s frame, the code then renormalizes the spectra by the factor (cosθ/πD²), as done in Equation (A.5), and also modifies the spectra according to all general relativistic effects. The normalization parameter (=(Mpc/D)² = 1.11 × 10⁵ for D = 3 kpc) in the kynrefionx model incorporates distance to the source and inclination is provided in the model as another input parameter, which we fixed to 66°. The model spectra produced by the above procedure are then used by XSPEC to fit the observed RS. The above modeling procedure enables us to fit the flux of the observed RS along with its shape because the reflection spectrum itself (from kynrefionx) is normalized appropriately and removes all arbitrariness in the normalization. All of these advantages lead to our final model: model B = TBabs×(cutoffpl + kynrefionx(reflionx_hd) + kynrefionx(reflionx_hd)).

The best-fit model B gives the density and ionization of reflection modeling, which are consistent in terms of both the flux and the shape, including all the general relativistic effects and radial variation of ionization.

It is also possible to fit both the shape and flux of the spectra keeping the model the same as in model A (i.e. reflection ≡ relconv⊗reflionx_hd), but extending the energy range from 1 eV to 1 MeV, along with the iteration of the fitted density and ionization parameter. Although this procedure is not as accurate as that of model B, it provides more flexibility to extend the same procedure to different reflection tables (if available) in the full energy range from 1 eV to 1 MeV. The method of iteration to find a consistent fit with this model and its fitting results, which are very similar to the results with model B are detailed in Appendix A.

3.2.2. Results

Figure 3 and Table 1 give the best-fit results of NuSTAR FPMA data with model B. The presence of a reprocessed blackbody hump in the soft X-ray region (especially for the inner reflection) can be noticed easily due to the higher density value. To find the consistency of reprocessed blackbody with observed data, a simultaneous NICER observation in 0.5–10.0 keV is also very helpful; its unfolded spectrum is shown in Fig. 3. Even then, the reprocessed blackbody radiation does not lead to any problems; that is, the model prediction is below the observed data, as opposed to what is shown in Fig. 2 of current work or in Fig. 3 of ZDM20. This result implies that using the high-density tables, we can fit the low/hard state of MAXI J1820+070 self-consistently. Although the model fits the energy range below 40 keV very well, some residual remains in the higher energy. However, the significance of the data above 40 keV is low and possibly other adjustments might improve the high energy fitting (Zdziarski et al. 2021b; Marino et al. 2021); namely, using nthcomp to model the primary hard X-ray or making the cutoff energy free in the existing cutoffpl model. The most crucial point to notice is the inner radius of the disk (r_in), which turns out to be 2.1 R_g. This exercise concludes that even within the curved space-time regime, the high-density disk can reach very near to ISCO, remaining consistent with the observed flux and shape simultaneously in low/hard state of XRBs. As the observed spectra is very similar to the typical spectra in the low/hard state of XRBs, we can generally say that it is most likely the reprocessed blackbody due to disk reflection is consistent, even if the disk extends close to the BH. We note that the compatibility of simultaneous NICER observation with the model spectra (from fitting only NuSTAR FPMA) in the lower energy in Fig. 3 is sufficient to draw a conclusion on the consistency of the reprocessed blackbody with the reflection. The primary reason for not fitting NICER data is the requirement of more model components (e.g., gabs, edge) to fit the different instrument’s calibration features at low energy, which would make the modeling unnecessarily complicated. For completeness, in Appendix B, we fit the NICER data as well along with NuSTAR FPMA and FPMB.

Fig. 3. Best-fit model B for the NuSTAR observation with reflection ≡ kynrefionx(reflionx_hd). The fitted model parameters are presented in Table 1. Simultaneous NICER observation is also overplotted in the energy range of 0.5–10.0 keV which is sufficient to understand the consistency or inconsistency of reprocessed blackbody. The primary reason for not fitting NICER data is the requirement of more model components (e.g., gabs, edge) to fit the different calibration features of the instrument at low energy which will make the model unnecessarily complicated. For completeness, in Appendix B we fit the NICER data as well with NuSTAR. Please note that we are showing unfolded spectra instead of actual data to match the Y-axis for two different instruments. Data were rebinned using the command “setplot rebin 100 10” to present in the plot.

As the disk approaches closer to the BH, fixing emissivity index to q = 3 may not be a good assumption for fitting the inner reflection component (Svoboda et al. 2012; Dauser et al. 2013). If we keep q for the inner reflection free, then the fit is improved slightly and the best fit gives q = 6.6, which is more acceptable (χ²/d.o.f. = 2145/1865). We also tried more realistic configuration by assuming radial variation of density of the disk same as Keplerian disk. We assumed the density of the disk for inner and outer reflection to vary as r^−1.65 and r^−1.875, respectively, following the middle and outer region of the Keplerian thin disk (Equations 3.66 and 3.68 in Kato et al. 2008). This imposition also improves the fitting slightly (χ²/d.o.f. = 2134/1866). However, with all these sophisticated adjustments, the fitting gives r_in ∼ 2 − 2.5R_g, keeping the flux distribution in different components same as for model B. Even if we assume the BH to be slowly spinning (a = 0.5), the fitted r_in becomes the same as the ISCO (4.23R_g for a = 0.5), making the fit slightly worse (χ²/d.o.f. = 2167/1866). Considering the flux contribution from different components, model B.1 (i.e., model B with a modification in the normalization of the inner reflection) is introduced (see Section 4.2.1) to make the geometry physically acceptable. Also with respect to model B.1, the disk extends very close to the BH, similarly to the results we see from model B.

4. Discussion

4.1. Reasons behind achieving a consistent good fit

A higher value for r_in typically (unless it is very close to r_out) requires the density and ionization to be lower to reach the same flux, if the other parameters remain the same (Equation (A.6)). This is why a greater flux in the outer reflection component still requires lower density and ionization than the inner reflection (Table 1). This is exactly in line with the logic put forward by ZDM20, and higher r_in values put the thermal blackbody radiation at lower energy. It must be emphasized that if we treat the reflection and reprocessed blackbody separately, and then add the resulting blackbody as a consequence of reflection on top of an already fitted model (Section 3.1), then it is very probable that model will overshoot the observed spectra (Fig. 2 in this paper, Fig. 3 in ZDM20), making the reprocessed blackbody inconsistent with the observed reflection. On the contrary, fitting the shape and flux simultaneously, while including the reprocessed blackbody (Section 3.2), results in very good fit (Fig. 3, Table 1) even if the disk extends very close to ISCO. Another reason for the appearance of thermal hump in ZDM20 is to the connection between primary and reflection, which we typically ignore when we fit the spectra, as done in this study due to uncertainty in the geometry of the corona and its orientation with respect to the disk. The best-fit result with model B suggests that the inner reflection component is heavily ionized (Fig. 3, Table 1) and practically does not give any Fe line feature. Nevertheless, it plays a key role in the fit by providing some curvature near the Fe line and in the overall continuum. It is evident that if the shape of the spectra demands high ξ and n, then naturally r_in has to be low enough to match the flux requirement (Equation (A.6)). The perfect balance between ξ, n, and r_in gives a fitting result that consistent in terms of both the shape and the flux.

4.2. Geometry of the source

Although our best fit suggests that the disk is extended very close to the BH, making the spectral fits consistent from both the shape and the flux, it does not tell us anything about the geometry of the source. The flux present in different spectral components can guide us in understanding this aspect. Due to the large energy band (1 eV–1 MeV), the flux represents the total radiation available in respective components. From the fitting results of model B (Table 1), we find that the total flux in inner and outer reflection is 3.2 × 10⁻⁸ ergs.cm⁻².sec⁻¹ which is 0.64 times of the flux in power-law component. This fact indicates that the flux reaching to the observer from the corona is greater than that reaching to the disk. Another very puzzling fact is that the flux in outer reflection component is four times the flux in inner reflection component. Stronger Compton hump for the outer reflection compared to inner reflection in many other fitting results also indicates larger flux in outer reflection component compared to inner reflection in energy range, 1 eV–1 MeV (for example Fig. 5 in Zdziarski et al. 2021b for epochs 1–4; Fig. 4 in Marino et al. 2021 for epoch 2; Fig. 10 in Liu et al. 2023 for source V404 Cyg). As all the released gravitational energy will be concentrated in the inner region due to proximity to the BH, it is natural to think that the power-law continuum that is the most dominant component of the spectra will also be concentrated in the inner region. However, with the centrally concentrated coronal radiation, it is difficult to imagine that the inner region of the disk gets only ∼1/10th of the power-law flux¹⁰. The simplest possibility for resolving the issue of very low inner reflection is to assume that the inner disk is being covered by the corona, which reduces the observed flux in the inner reflection component.

4.2.1. Fitting that includes the covering of the inner disk

If the inner region of the disk is covered by the optically thick corona, we would expect that most of the inner reflection component is up-scattered by the corona and contributes to the power-law component. Thus, we introduced our model B.1 which is based on model B, but with the normalization of inner reflection reduced. For model B, the normalization is fixed according to the distance to the source; whereas for model B.1, the normalization is further reduced to 10% of this value to mimic the covering of the inner disk by the corona (10% is chosen as it is approximately the ratio of the observed power-law flux and inner reflection flux in model B). The best-fit parameters are given in the last column of Table 1. We were able to find an equally good fit with model B.1, but with ten times higher luminosity and ionization of inner reflection compared to model B. It indicates the crucial fact that the fitting of inner reflection is not driven by a relativistic iron line, but by its contribution to the continuum (as seen in Fig. 3). This is due to already very high inner disk ionization in model B. By increasing the ionization further, the inner reflection in model B.1 displays a similar spectral shape, but higher normalization. The goodness of the fit with model B.1 indicates that the inner disk covered with an optically thick corona can offer a plausible geometry for the observation analyzed in this work. Although, this geometry indicates that the disk is extended close to the BH, the obscuration of a major part of the inner disk by corona resembles a truncated disk geometry.

4.2.2. Comptonization of inner reflection using MONK

To check the possibility of reduction of inner reflection by its Comptonization in the slab corona above the disk, we used the relativistic Monte Carlo radiative transfer code MONK (Zhang et al. 2019). In the hard state, we did not see any disk blackbody emission due to intrinsic dissipation. Therefore, we assumed only the Comptonization of the inner reflection component here, which also includes the reprocessed blackbody. To keep things simple, we restricted our approach to the Newtonian regime and provided the unabsorbed inner reflection component of model B.1 (Section 4.2.1) as the seed spectrum for Comptonization. We assumed that the disk is emitting the seed spectrum as a cosine source and the thin slab of hot electrons as corona is present just above the disk. The electron temperature and optical depth of the corona are set according to the requirement of two factors: (i) the observed power law and (ii) the ratio of observed flux in inner reflection and power-law component in model B.1, which is 0.1 at 66° inclination of the source. The unscattered flux that leaks through the corona becomes the observed inner reflection component. The result of this experiment is presented in Fig. 4(A), as broadband spectral output from MONK. To achieve the observed constraints, we require the optical depth of the corona to be 0.282 and temperature of the electrons to be 160 keV. The integrated term is used to denote the total spectrum integrated over 2π solid angle (upper hemisphere above the disk or lower hemisphere below the slab corona). Otherwise, the spectrum at an inclination of 66° is indicated specifically. The components of unabsorbed model B.1 (by setting N_H = 0) are overplotted with corresponding MONK outputs in Fig. 4(B). It shows that the above-mentioned coronal properties reproduce the observed power law quite accurately. Although there is some change in the shape of the unscattered spectrum compared to the incident spectrum, it appears that the covering of the disk by such corona can credibly lead to the required reduction in the observed inner reflection flux.

Fig. 4. Comptonization of inner reflection photons of model B.1 through slab corona of optical depth 0.282 and electron temperature 160 keV using MONK. Panel (A) shows the broadband spectral output from MONK. Integrated term is used to denote total spectrum integrated over 2π solid angle, whereas 66° inclination is mentioned to indicate observation at that specific inclination. Seed spectrum is provided as the seed photons for Comptonization in MONK. Comptonized spectrum is made of photons which are scattered one or more times by the hot electrons in the corona. Unscattered spectrum is made of photons which are able to leak through corona without scattering. Flux is always computed in energy range 1 eV–1 MeV. Panel (B) shows the unscattered spectrum and Comptonized spectrum of MONK output at inclination 66° along with model components of the best fitted model B.1 presented in Table 1.

However, although we would be able to achieve the observed power law and inner reflection with the slab corona lying above the disk, it is crucial to check whether this configuration can stand as a possible steady-state solution or not. In the MONK simulation, we started with the seed spectrum emitted by the disk for Comptonization (pink dashed line in Fig. 4(A)). As a result, not only we get the observed power law, observed unscattered inner reflection, but also a Comptonized spectrum illuminating the disk (yellow dashed line in Fig. 4(A)). This illumination will be reprocessed by the disk and emitted as the seed spectrum for Comptonization in the next step. A steady-state solution demands that the total integrated flux (in energy band 1 eV–1 MeV) of the Comptonized power law illuminating the disk should be the same as the seed spectrum for Comptonization; however, this is not the case here. The Comptonized power law illuminating the disk integrated over 2π solid angle (3.69 × 10⁻⁷ ergs cm⁻² s⁻¹) has 2.75 times more energy flux compared to the integrated seed spectrum (1.34 × 10⁻⁷ ergs cm⁻² s⁻¹), also denoted in the top left corner of Fig. 4(A). Thus our model predicts more illumination, thus, more seed photons than we originally used to produce the power law. This indicates that static slab corona lying above the disk can not be a steady-state equilibrium solution to explain the required reduction of inner reflection component and some kind of anisotropy in the coronal emission may solve the problem as discussed in next section. To find the effect of Comptonization by different coronal geometries for spectra in hard state, we also refer to Poutanen et al. (2018).

4.2.3. Possibilities for required anisotropy

The discrepancy between the Comptonized flux illuminating the disk and the seed spectrum required for Comptonization (described in Section 4.2.2) could, in principle, be overcome by introducing some anisotropy in the emission of the corona. One such possibility would be outflowing corona which is in line with the observed jet for this source (Bright et al. 2020). Following the procedure from Section 4.8 of Rybicki & Lightman (1979), we find that moving corona with velocity ∼0.2 − 0.3c (c is the speed of light) perpendicular to the disk will lead to the integrated seed flux received by the corona from the disk reduced by ∼2 − 3 relative to the flux emitted by the corona toward the disk. From MONK simulation with outflowing slab corona covering the disk, we also find 0.3c velocity is sufficient to make the model as a steady-state solution for this observation. Beloborodov (1999) carried out a similar analysis of bringing anisotropy in the coronal emission through its bulk motion to explain the observed power law and low reflection fraction in the hard state of Cyg X-1, which resulted in a velocity of ∼0.3c. Recent observations of polarization from Cyg X-1 even requires the slab corona to be moving at ∼0.6c (Poutanen et al. 2023). If we did not consider any obscuration of the disk by the corona like in a lamp-post geometry, the velocity of the corona would have to be ∼0.7 − 0.8c to give the observed fraction of inner reflection flux relative to power law at 66° inclination (0.13 for model B, 0.1 for model B.1).

4.2.4. Past analyses and comparison with current work

As our primary objective here is reflection, we fit the RS considering its shape and flux, and we kept its normalization disconnected from the power law to avoid the uncertainty in disk illumination by the corona due to its uncertain geometry and/or velocity. It is only after we find the best fit, that we can discuss whether it can be interpreted with some corona geometry. In a similar line, there are many works that have fitted the observed data from the same source with arbitrary normalization of the power law and reflection components; even without worrying about the consistency of reflection from its shape and flux. Previous analyses of MAXI J1820+070 through its different states of the outburst were studied by Buisson et al. (2019), Chakraborty et al. (2020), with two-component corona in lamp-post geometry illuminating the inner and outer disk. They have also found that an accretion disk may reach very close to the BH. Although their outer reflection remain lower compared to inner reflection as expected with centrally concentrated corona, their inner reflection is still too low with respect to the power-law component (first panel, Fig. 6 in Buisson et al. 2019¹¹), similar to our result. This is impossible to achieve with static lamp-post corona above the disk.

On the other hand, Zdziarski et al. (2021b), who also considered two power-law corona components, reached the conclusion that the disk had to be truncated. In their work, they fit the same data with more acceptable values of inclination (the same as used in our current work) and iron abundance. The recent broadband study of Banerjee et al. (2024) also offered the same conclusion. The physical model of a jet emitting disk (JED) acting as the corona (Marino et al. 2021) fit the data equally well with truncated disk geometry. In this scenario, hot advective accretion flow or JED serves as the corona, which lies in the same plane of the disk. The solid angle subtended by the disk to the corona is smaller compared to the 2π solid angle of the upper hemisphere, resulting in a lower observed flux in the reflection component than in the power-law component. However, to achieve more flux in outer than inner reflection component, some obscuration of the inner disk region and/or flaring of the outer disk is still needed as suggested in Fig. 4 of Zdziarski et al. (2021b) or in Fig. 10 of Marino et al. (2021).

The observed timing properties of this source in the low/hard state show large lags between the coronal emission and the disk reflection, which Kara et al. (2019) interpreted with an elongated corona (modeled as two component lamp-post corona) illuminating a disk reaching close to the BH. The same lags can be explained also by the illumination of the outer disk by the central corona in the truncated disk geometry (De Marco et al. 2021; Kawamura et al. 2022). Although, we did not try to model the timing properties, from the spectral modeling, we propose a possible geometry of outflowing slab corona covering the disk. The inner region of the disk remains covered by the slab corona and, from the observed flux, we can see that the outer reflection is more dominant than the inner reflection. Therefore, the proposed geometry becomes very similar to the truncated disk geometry; as a consequence, this could, in principle, explain the observed lag between the power law and the reflection.

Although all the works mentioned above are able to fit the data component wise, none of them have fitted the reflection component by simultaneously considering its shape and flux, as done in this work. Additionally, the proper modeling of Comptonization and obscuration by the corona was only carried out only in few works to describe the observations, for instance, Steiner et al. (2017), Poutanen et al. (2018) (for AGNs, see Petrucci et al. 2001; Malzac et al. 2005; Wilkins & Gallo 2015). To investigate further about the geometry with appropriate modeling of Comptonization, we need to include more coronal geometries in the existing MONK code; for instance, outflowing slab corona or truncated disk geometry. The Comptonization of the reflection spectrum in the general relativistic framework also has to be further developed.

4.3. Limitations of our study and possibilities for its extension

Naturally, there are a few limitations of our study that need to take into account to understand the robustness of the results. The atomic database of reflionx_hd table is not as rich as in XILLVER tables. The fixed high energy cutoff value of reflionx_hd table used in this work also restricts us to a narrow parameter regime. In addition, there can be two different Comptonizing components giving rise to two different power law, as investigated in Zdziarski et al. (2021b). The self-consistent normalization of the RS from both the flux and the shape is possible only because of the well known distance and inclination of the source. We do not have this freedom for most of the observed XRBs. To study the geometry, we only considered the Comptonization of the inner reflection emitted by the disk. A self-consistent modeling of coupled disk-corona system with possible energy exchange between them is beyond the scope of this work. The incorporation of intrinsic viscous dissipation in the fitting also demands development of new reflection table. Above all, we have analyzed only one spectra in low/hard state of one source in this study. In addition to addressing the difficulties presented in this work to explain the spectral flux distribution in fitted spectral components, we also need to study the observed timing properties to understand the true geometry of the source.

Despite of the above limitations, it will be interesting to extend our study to other XRBs in the low/hard state and find the inner radius of the disk. If the source has well-known distance, inclination, and mass values for the BH (e.g. V404 Cyg, Miller-Jones et al. 2009; Khargharia et al. 2010), then the simultaneous fitting of shape and flux of RS should give very accurate value of inner radius of the disk. Even though the distance and inclination of GX 339-4 is uncertain (Basak & Zdziarski 2016), constraining them to some acceptable range and simultaneous fitting the flux and shape of the RS for multiple observations together can shed some light on the inner radius of the disk in low/hard state for that source. Due to repeated outbursts of this source it is observed and studied many times. There is a large scatter about the inferred inner radius of the disk (Fig. 11 in García et al. 2015).

5. Conclusions

In this work, we have fit both the shape and the flux of the observed RS in the low/hard state of the XRB MAXI J1820+070. We have examined whether we can exclude the possibility of disk extending close to the BH from the perspective of reprocessed black body due to disk illumination. From the flux consideration we have attempted to constrain the geometry of the source. The key findings are summarized below.

• Consistently carrying out the fitting, with both the shape and the flux of the observed RS taken into consideration, gives the exact estimate of the density and the ionization of the disk. The consistent best fit with one power-law component suggests that the disk is extended very close to the BH in the low/hard state of XRB. The fitting of inner reflection component gives quite high disk ionization; thus, it is primarily driven by reflection continuum shape and its contribution to the total flux rather than by the shape of the relativistic iron line.

• The density and ionization of the disk are high enough to give rise to the observed reflected flux and low enough to make the reprocessed blackbody negligible. It is not possible to exclude the possibility for disk extending very close to the BH in low/hard state from the perspective of reprocessed blackbody.

• The total flux in reflection is smaller than the flux in the power-law component. At the same time, the inner reflection is much lower than the outer one, which is impossible for a centrally located corona. We show that covering of the inner disk by static slab corona may lead to the necessary level of the power-law component as a result of the Comptonization of re-processed emission and, at the same time, to the apt level of transmitted unscattered inner reflection. On the other hand, the flux in generated power-law component illuminating the disk is higher than the assumed seed photon flux for Comptonization.

• An outflowing corona with some bulk velocity (∼0.2 − 0.3c) could, in principle, reduce the illumination of the disk and lead to the same flux in the seed spectrum seen for Comptonization and in a Comptonized power law emitted by the corona toward the disk. This has to be investigated thoroughly with the appropriate Comptonization in the outflowing corona.

Acknowledgments

SRD and MD thank Andrzej A. Zdziarski and Chris Done for many stimulating scientific discussions. SRD also thanks Peter Boorman, Jiří Svoboda and Craig A Gordon for their help in learning XSPEC. The authors also thank the reviewer for many critical constructive comments which make the conclusion of the manuscript more robust. SRD, MD, MB, JH and VK thank GACR project 21-06825X for the support, as well as the institutional support from RVO:67985815. WZ acknowledges NSFC grant 12333004, and support by the Strategic Priority Research Program of the Chinese Academy of Sciences, grant no. XDB0550200.

References

Appendix A: Fitting shape and flux simultaneously through iteration of density

A.1. The model

Besides model B (presented in section 3.2), here we build another alternative model with which also we can fit simultaneously the shape and flux of the RS (reflection spectra). Here reflection is modeled by (relconv⊗reflionx_hd) as done in model A (Section 3.1). Fitting the RS using the model (relconv⊗reflionx_hd) gives values of ξ, n and r_in of the disk. reflionx_hd model in XSPEC has normalization parameter which we keep free while fitting. First we fit the spectra keeping the density parameter free which essentially fits only the shape of the spectra because using arbitrary density and normalization XSPEC rescale the model spectrum to the desired flux value. For fixed values of parameters (e.g. photon index, ξ, n) of reflionx_hd table, it gives a specific reflection spectrum. Now when we apply the relconv model on reflionx_hd table, it does the following: (i) takes that specific spectrum depending on the input parameters of the reflionx_hd table, (ii) adds contribution from each radius according to r^−q, i.e., emissivity index parameter of relconv model, (iii) integrates these radial flux contributions over the disk area from r_in to r_out after incorporating all the relativistic effects. Therefore, in XSPEC, this final model spectrum is used to fit the observed reflection spectrum whenever we use (relconv⊗reflionx_hd) to model the reflection. Note that although the relconv model takes radial flux contribution according to emissivity profile, it takes only one single spectrum from reflionx_hd table corresponding to one ξ value for all radii. So, the fitted value of ξ gives an average value of ξ(r) which best fits the shape and we term it as ξ_shape.

From the point of view of the flux, the assumed radial emissivity profile (Equation 2) demands the same radial profile of ξ(r) as we are not considering any radial variation of density,

$$\begin{array}{c}\hfill \xi (r)={\xi }_1{r}^{-q},\end{array}$$(A.1)

where ξ₁ is the ionization at r = 1 and r is always in units of R_g. Plugging the radial profile of ξ in Equation 1:

$$\begin{array}{c}\hfill F_{\mathrm{refl}}(r)=\frac{n\xi (r)}{4\pi }=\frac{n{\xi }_1}{4\pi }{r}^{-q}.\end{array}$$(A.2)

It is important to remember that the above relation holds true because the slab is in radiative equilibrium with the incident radiation which makes the total reflected flux in the energy range 1 eV - 1 MeV equal to the incident flux. As ξ_shape represents some average value of ξ(r) fitting the shape of the spectra, to make a comparison with it, we must find some average value of ξ(r) corresponding to flux, namely, ξ_flux. To find that, we take flux weighted average of ξ over the disk area, which gives

$$\begin{array}{cc}\hfill {\xi }_{\mathrm{flux}}& =\frac{{\displaystyle {\int }_{r_{\mathrm{in}}}^{r_{\mathrm{out}}}}\xi (r)F_{\mathrm{refl}}(r)2\pi rdr}{{\displaystyle {\int }_{r_{\mathrm{in}}}^{r_{\mathrm{out}}}}{F}_{\mathrm{refl}}(r)2\pi rdr}=\frac{{\displaystyle {\int }_{r_{\mathrm{in}}}^{r_{\mathrm{out}}}}{\xi }_1{r}^{-q}\left(\frac{n{\xi }_1}{4\pi }\right)r^{-q}2\pi rdr}{{\displaystyle {\int }_{r_{\mathrm{in}}}^{r_{\mathrm{out}}}}\left(\frac{n{\xi }_1}{4\pi }\right)r^{-q}2\pi rdr}\hfill \\ \hfill & ={\xi }_1\left(\frac{2-q}{2-2q}\right)\left(\frac{r_{\mathrm{out}}^{2-2q}-r_{\mathrm{in}}^{2-2q}}{r_{\mathrm{out}}^{2-q}-r_{\mathrm{in}}^{2-q}}\right).\hfill \end{array}$$(A.3)

The total radiation emitted by one side of the disk is

$$\begin{array}{cc}\hfill L_{\mathrm{refl}}& ={\displaystyle {\int }_{r_{\mathrm{in}}}^{r_{\mathrm{out}}}}{F}_{\mathrm{refl}}(r)2\pi rdr\times R_{\mathrm{g}}^2=\frac{{\xi }_{1}n}{2(2-q)}(r_{\mathrm{out}}^{2-q}-r_{\mathrm{in}}^{2-q})\times R_{\mathrm{g}}^2\hfill \\ \hfill & ={\xi }_{\mathrm{flux}}n\frac{(1-q)}{{(2-q)}^2}\frac{{(r_{\mathrm{out}}^{2-q}-r_{\mathrm{in}}^{2-q})}^2}{(r_{\mathrm{out}}^{2-2q}-r_{\mathrm{in}}^{2-2q})}\times R_{\mathrm{g}}^2.\hfill \end{array}$$(A.4)

At inclination θ and distance D from the source, the observed reflected flux will be

$$\begin{array}{c}\hfill F_{\mathrm{refl},\mathrm{o}}=L_{\mathrm{refl}}\times \left(\frac{cos\theta }{\pi D^2}\right).\end{array}$$(A.5)

For the limiting case of r_out ≫ r_in and q = 3, ${\xi }_{\mathrm{flux}}={\xi }_1{r}_{\mathrm{in}}^{-3}/4$, L_refl = 2ξ_fluxn(r_inR_g)², leading to

$$\begin{array}{c}\hfill F_{\mathrm{refl},\mathrm{o}}=2{\xi }_{\mathrm{flux}}n{(r_{\mathrm{in}}{R}_{\mathrm{g}})}^2\times (cos\theta /\pi D^2).\end{array}$$(A.6)

From the fitting of observed spectra, we find the unabsorbed reflected flux within the energy range 1 eV - 1 MeV using the flux command in XSPEC, which gives true F_refl, o. This large energy range covers essentially all possible radiation emitted by the slab. Note that this becomes possible only because the energy grid of the reflionx_hd table covers this large range. Using Equation (A.4) & (A.5), from F_refl, o we find ξ_flux which should match with ξ_shape from fitting. For n = 10¹⁵ cm⁻³, always ξ_flux turns out to be larger than ξ_shape. From Equation (A.4), it is clear that by increasing n, we can reduce ξ_flux and meet the demand of high flux while ξ_flux = ξ_shape. After fixing the n corresponding to ξ_flux = ξ_shape, we again fit the data since n changes the spectral shape. ξ_shape and r_in are modified due to fitting and again new n is required to converge. We iterate until a fit is achieved when ξ_shape become within the 5% range of ξ_flux. Therefore, model C for our fitting becomes model C = TBabs×(cutoffpl + relconv⊗reflionx_hd (ξ_flux = ξ_shape)+ relconv⊗reflionx_hd (ξ_flux = ξ_shape)). The final fit that follows ξ_flux = ξ_shape will possibly not be as good as when there is no restriction in density. Through iteration of density we are able to match ξ_shape and ξ_flux, and thus fit the observation simultaneously by its shape and flux.

Fig. A.1. Reasonably good-fit model C for NuSTAR FPMA observation along with its components when reflection ≡ relconv⊗reflionx_hd. Densities for reflection components are iterated to make the fitting consistent from shape as well as flux. The fitted model parameters are presented in Table A.1. Simultaneous NICER observation is also overplotted in the energy range of 0.5-10.0 keV which is sufficient to understand the consistency or inconsistency of reprocessed blackbody in soft X-ray. The best fitted model components turn out to be very similar with the results of model B which is shown in Fig. 3.

A.2. Results

For the fitting with model C, densities for both the reflection components were first kept free and the best fit was achieved. Progressing toward consistency, we varied the density and fitted iteratively, as described at the end of section A.1. We fixed the ionization to the most neutral case provided by the tables, ξ = 100 for the outer reflection component. The reasonable good fit model is shown in the Fig. A.1, and the relevant parameters are presented in Table A.1. The statistics of the consistent fit is χ²/d.o.f = 2141/1865, slightly worse than the best fit when densities were kept free (χ²/d.o.f = 2123/1865). The density to make ξ_flux = ξ_shape for inner and outer reflection becomes 1.5 × 10¹⁹ cm⁻³ and 1.1 × 10¹⁸ cm⁻³, respectively. The iteration in model C is delicate due to the strong dependence of flux on r_in, n, and ξ. Therefore, once we achieve a reasonably good fit, we avoid finding the best fit. Consequently, we cannot estimate the errors in the model parameters using XSPEC because that is possible only after achieving the best fit. That is the reason why no error for parameters of model C is presented in Table A.1. With the consistent parameters, the theoretically calculated flux remains within the 5% range of unabsorbed flux estimated through XSPEC for the energy range 1 eV - 1 MeV (based on the criterion we follow for convergence after iterations as mentioned in Section A.1). This confirms that the fit is almost consistent simultaneously with the shape and flux of the observed RS. The fitting results are presented in Table A.1 and in Fig. A.1 which seem very similar to the results of more accurate model B presented in Section 3.2.

Appendix B: Fitting NuSTAR and NICER together for MAXI J1820+070

Here we present the fitting of quasi-simultaneous data of NuSTAR FPMA, NuSTAR FPMB and NICER with self-consistent normalization considering both the shape and the flux of the RS. The energy range for fitting NuSTAR data is 3.0-78.0 keV whereas for NICER we fit 0.5-10.0 keV. One such model with self-consistent normalization of RS is model B (Section 3.2.1). Fitting of only NuSTAR FPMA data with model B is presented in Fig. 3 and Table 1 (Section 3.2.2). To fit the data from different instruments, we need to incorporate a constant factor for cross-normalization. In addition, following previous works (Wang et al. 2020, 2021; Ludlam et al. 2021, 2022; Kumar et al. 2022; Feng et al. 2022; Moutard et al. 2023; Svoboda et al. 2024b,a), to handle the different calibration features present in the low energy data from NICER, we incorporate a Gaussian absorption line (gabs) and absorption edge (edge) to fit the line feature below 1 keV and edge around 2.4 keV. These models (gabs and edge) are applied only to fit the NICER data. This makes model B.2 for fitting NuSTAR + NICER data together. model B.2 = constant×gabs×edge×model B = constant × gabs × edge × TBabs × (cutoffpl + kynrefionx(reflionx_hd) + kynrefionx(reflionx_hd)).

The relevant parameters of the best fitting are mentioned in Table B.1. Corresponding unfolded spectra and the residuals are plotted in top and bottom panel of Fig. B.1. The best fitting leads to χ²/ν = 4519/3904 ∼ 1.16 (for NuSTAR FPMA χ² = 2246 for 1874 bins, for NuSTAR FPMB χ² = 2165 for 1874 bins, and for NICERχ² = 108 for 171 bins respectively). Note that the χ² value for NICER becomes much less than the number of bins which is basically due to added 1.5% systematic error while reducing the level 3 data following complete product pipeline of NICERDAS. On the contrary, if we do not add gabs or edge component then the χ² value becomes reasonable but residuals clearly indicate the presence of such components.

Residuals in Fig. B.1 (bottom panel) indicate some feature still present around iron line emission. This is possibly due to some abrupt change in the iron line shape in reflionx_hd table with change in density or ionization. As an example, notice the change in iron line shape around 6.4 keV in the right panel of Fig. 1 for ξ = 278.3 with respect to ξ = 100.0 and ξ = 774.3. This kind of diminished iron line can be noticeable for some other values of parameters also. We verify this argument by including another Gaussian line with existing two reflection components in fitting which diminishes such feature present in residuals significantly. Coarser energy grid of the reflionx_hd table (energy spacing between two grid points is around 0.1 keV at 6.4 keV) may have some effect in fitting iron line. Additionally, the data we are fitting from NICER and NuSTAR are not exactly simultaneous. This non-simultaneity also can contribute to the residuals as the fitting is done assuming same model parameters for the observations from two instruments. Despite the fact that the iron line is not fitted perfectly well for simultaneous data from NICER and NuSTAR, reprocessed blackbody from reflection still become sufficient to fit the low energy data and remain consistent with the non observation of any other blackbody component in the hard state. Compared to fitting only the NuSTAR FPMA data (Table 1, model B), primarily the density (n) and ionization (ξ_in) of the inner reflection component change to fit the spectra of NICER in low energy along with NuSTAR. There is also slight change in ionization and density of outer reflection. The flux in different components remain almost same. This result denotes that density can also play a crucial role in fitting the broadband RS, along with the ionization parameter.

Fig. B.1. Spectral fit for NuSTAR FPMA (magenta), NuSTAR FPMB (red) and NICER (blue) with model B.2. Top panel shows the unfolded spectra and bottom panel shows the data residuals. Although residuals indicate that this is not a perfectly fitted iron line, the reprocessed blackbody remain sufficient and consistent to fit the data from NICER along with NuSTAR. The fitted model parameters are presented in Table B.1. Data is rebinned using the command “setplot rebin 100 10” to present in the plot.

All Tables

All Figures

	Fig. 1. Variation in EFE with density (n) and ionization (ξ) from reflionx_hd table with appropriate normalization. The corresponding n and ξ are written in the legend. The reflionx_hd table has three input parameters: n, ξ, and Γ. For a variation of either n or ξ, the fixed values of other input parameters are noted on the plot. More on the normalization of these spectra in footnote 9.
In the text

Fig. 2. Best-fit model A for NuSTAR data. Expected reprocessed blackbody as a consequence of reflection are modeled as diskbb components, which are added separately over the best-fitted model A. diskbb(inner) and diskbb(outer) denote the diskbb components corresponding to thermalization due to inner and outer reflection, respectively. Simultaneous NICER data were presented to compare the fitted model in the lower energy. We present the unfolded spectrum to show the data of different instruments on the same scale.

In the text

Fig. 3. Best-fit model B for the NuSTAR observation with reflection ≡ kynrefionx(reflionx_hd). The fitted model parameters are presented in Table 1. Simultaneous NICER observation is also overplotted in the energy range of 0.5–10.0 keV which is sufficient to understand the consistency or inconsistency of reprocessed blackbody. The primary reason for not fitting NICER data is the requirement of more model components (e.g., gabs, edge) to fit the different calibration features of the instrument at low energy which will make the model unnecessarily complicated. For completeness, in Appendix B we fit the NICER data as well with NuSTAR. Please note that we are showing unfolded spectra instead of actual data to match the Y-axis for two different instruments. Data were rebinned using the command “setplot rebin 100 10” to present in the plot.

In the text

Fig. 4. Comptonization of inner reflection photons of model B.1 through slab corona of optical depth 0.282 and electron temperature 160 keV using MONK. Panel (A) shows the broadband spectral output from MONK. Integrated term is used to denote total spectrum integrated over 2π solid angle, whereas 66° inclination is mentioned to indicate observation at that specific inclination. Seed spectrum is provided as the seed photons for Comptonization in MONK. Comptonized spectrum is made of photons which are scattered one or more times by the hot electrons in the corona. Unscattered spectrum is made of photons which are able to leak through corona without scattering. Flux is always computed in energy range 1 eV–1 MeV. Panel (B) shows the unscattered spectrum and Comptonized spectrum of MONK output at inclination 66° along with model components of the best fitted model B.1 presented in Table 1.

In the text

Fig. A.1. Reasonably good-fit model C for NuSTAR FPMA observation along with its components when reflection ≡ relconv⊗reflionx_hd. Densities for reflection components are iterated to make the fitting consistent from shape as well as flux. The fitted model parameters are presented in Table A.1. Simultaneous NICER observation is also overplotted in the energy range of 0.5-10.0 keV which is sufficient to understand the consistency or inconsistency of reprocessed blackbody in soft X-ray. The best fitted model components turn out to be very similar with the results of model B which is shown in Fig. 3.

In the text

Fig. B.1. Spectral fit for NuSTAR FPMA (magenta), NuSTAR FPMB (red) and NICER (blue) with model B.2. Top panel shows the unfolded spectra and bottom panel shows the data residuals. Although residuals indicate that this is not a perfectly fitted iron line, the reprocessed blackbody remain sufficient and consistent to fit the data from NICER along with NuSTAR. The fitted model parameters are presented in Table B.1. Data is rebinned using the command “setplot rebin 100 10” to present in the plot.

In the text