© The Authors 2026

1. Introduction

The first successful model describing angular momentum transport and the formation of an accretion disk, and produced multi-temperature blackbody radiation from low-mass X-ray binaries (XRBs), was formulated by Pringle & Rees (1972) and Shakura & Sunyaev (1973). This model was soon extended into a general relativistic framework (Novikov & Thorne 1973; Page & Thorne 1974), which remains the most widely used model for spectral modeling of accretion disks to this day (see Frank et al. 2002; Kato et al. 2008 for further details). However, this “cold disk model” cannot account for the higher-energy spectral data (≳10 keV) observed in XRBs (Lightman & Shapiro 1975). To explain the origin of hard X-rays, incomplete thermalization in an optically thin hot flow must be invoked (Shapiro et al. 1976). To render this hot flow thermally stable, some degree of advection is required (Abramowicz et al. 1995; Chakrabarti & Titarchuk 1995; Narayan & Yi 1995). In the unified picture of accretion flows, a truncated disk geometry is therefore proposed, where an optically thick, geometrically thin disk is truncated at a radius inside of which a hot, advective, optically thin, and geometrically thick flow exists. The truncation radius moves inward with increasing accretion rate (Esin et al. 1997). The hot electron medium responsible for producing a hard power law spectrum via inverse Comptonization is referred to as the corona. This truncated geometry is supported by numerous observational studies, both spectral and timing-based (Basak & Zdziarski 2016; Marcel et al. 2019; Zdziarski & De Marco 2020; Marino et al. 2021; Banerjee et al. 2024; Chand et al. 2024; De Marco et al. 2021; Kawamura et al. 2022).

Although truncated disk geometry remains a plausible configuration for accretion flows, the observation of broad fluorescent iron lines from XRBs, and from active galactic nuclei (AGNs), suggests that the optically thick disk often extends close to the black hole (BH) and is illuminated by a hard X-ray–emitting corona even in the hard state (Fabian et al. 1989; Ross et al. 1996; Miller et al. 2006) (for more recent observations, see Wang-Ji et al. 2018; Buisson et al. 2019; Chakraborty et al. 2020; Sridhar et al. 2020; for AGNs, see Tanaka et al. 1995; Nandra et al. 1997; Reynolds & Begelman 1997; Fabian et al. 2000). If the disk extends close to the BH, then the corona must lie above the disk. Buoyant magnetic flares and an optically thin disk atmosphere produced via disk evaporation could also serve as a hot corona (Meyer et al. 2000). However, whether the inner region of the optically thick disk can evaporate completely, thus giving rise to the truncated disk scenario, remains an open area of research (Różańska & Czerny 2000; Liu et al. 2002; Cho & Narayan 2022; Bambic et al. 2024). Other geometric signatures, such as steep emissivity profiles or the reprocessed blackbody component from reflection, have not proven sufficient to distinguish among competing geometries (Svoboda et al. 2012; Datta et al. 2024). Recent polarization measurements in the hard state of XRBs suggest that the corona is extended horizontally, perpendicular to the disk axis (Krawczynski et al. 2022; Ingram et al. 2024; Podgorný et al. 2024; for a review, see Dovčiak et al. 2024). In this work, we also adopt the assumption of a two-phase accretion flow, consisting of a hot, extended slab corona situated above the disk.

For typical parameters of a standard Shakura–Sunyaev disk in X-ray binaries (XRBs), the thermal and viscous timescales are on the order of minutes and days, respectively (Pringle 1981). These are much longer than the light travel time between the disk and the corona, which is approximately GM/c³ ∼ 10⁻⁵ seconds (where G is the gravitational constant, M is the black hole mass, and c is the speed of light). Moreover, the typical timescale for XRB observations is ≳10 ks. Therefore, in a two-phase accretion flow, we expect the disk and corona to achieve equilibrium through mutual interaction.

For this work we considered the disk–corona connection solely through two processes: the supply of seed photons from the disk for Comptonization, and illumination of the disk by the corona, an approach primarily motivated by observational evidence. The first modeling of such a two-phase accretion flow in mutual equilibrium was presented by Haardt & Maraschi (1991, 1993) (hereafter HM93), and was later developed further by Stern et al. (1995), Poutanen & Svensson (1996), and Malzac et al. (2001) (with a dynamic corona). Malzac et al. (2005) extended this by treating the disk as a slab with either constant density or constant pressure. The most recent advancement in this direction is by Poutanen et al. (2018), who incorporated reflection spectra using the XILLVER tables (García & Kallman 2010; García et al. 2013) to more accurately model reprocessing from cold media.

A key result emerging from these models is that, for a given coronal temperature and geometry, there exists an upper limit to the Comptonized power the corona can produce in equilibrium. Beyond this limit, no energy-balanced solution is possible, even if the entire accretion power is assumed to be dissipated in the corona. This implies that one cannot arbitrarily increase the optical depth at fixed coronal temperature to obtain a harder Comptonized spectrum. As the optical depth increases, the coronal output increases, which enhances disk illumination and, in turn, raises the seed photon flux for Comptonization. For a sufficiently hot and optically thick corona, this feedback creates a runaway cascade, preventing equilibrium from being established. This limiting behavior offers an important means of constraining the disk–corona geometry. The general conclusion is that the hardest observed spectra in XRBs cannot be explained by a static slab corona located above the disk. Viable alternatives include an outflowing corona (Malzac et al. 2001) or a truncated disk (Poutanen et al. 2018). In the truncated disk scenario, a sufficient supply of seed photons for inverse Comptonization must still be maintained, either via overlap between the disk and corona or through non-thermal processes such as bremsstrahlung or synchrotron emission originating within the corona itself.

All of the aforementioned studies assume a local energy balance, i.e., they do not consider radial variations in disk emission or feedback from the corona. Moreover, general relativistic effects are typically neglected, even though they are expected to be significant in the vicinity of the black hole. On a different front, recent developments such as the Monte Carlo radiative transfer code MONK (Zhang et al. 2019) allow the computation of spectra from disk–corona systems around black holes within a fully general relativistic framework. This code can accommodate a range of extended coronal geometries, including slab, spherical, wedge, and conical structures, both static and outflowing. However, it does not incorporate feedback from the corona to the disk, making it inconsistent with the expected light travel timescale between the two regions and limiting its self-consistency.

For this work we further developed the MONK code to incorporate the mutual interaction between the disk and the corona. This extension allowed us to find equilibrium solutions for a given black hole, disk, and coronal configuration, thereby establishing a robust framework for constraining system parameters from observations. A key distinction from previous studies lies in our treatment of the fraction of accretion energy dissipated in the corona. Earlier works generally assumed that all of the accretion energy is dissipated in the corona, and that the seed photons for Comptonization arise solely from the reprocessing of coronal illumination by the cold disk (see also Section 2.4 of Poutanen et al. 2018, where a small intrinsic disk dissipation component was introduced as a free parameter to achieve equilibrium). This approach was primarily aimed at modeling the hardest spectra achievable and demonstrating the limitations of static slab corona in explaining the hard state of accreting systems. Under that assumption, the equilibrium coronal temperature was derived for a given optical depth. In contrast, we compute the equilibrium solution for given disk and coronal properties by adjusting the fraction of accretion power dissipated intrinsically in the disk through viscous processes, parameterized by α, for a fixed coronal electron temperature and optical depth. This is done under the constraint that the total luminosity from the disk–corona system equals the total gravitational power available through accretion. This approach allows us to determine the physically consistent coronal luminosity fraction, not only for the hard state, but also for intermediate states, where a significant contribution from the intrinsic disk blackbody emission may be present.

We describe the complete structure of MONK step by step to find the equilibrium disk-corona solution in Sect. 2. In Sect. 3 we show the effects of imposing equilibrium on different spectral components and the limiting behavior on the hardest spectra possible to achieve. The effect of BH spin, coronal temperature, and disk albedo are studied. Finally after discussing a few crucial points and our limitations in Sect. 4, we conclude in Sect. 5.

2. The complete model

The Monte Carlo radiative transfer code MONK (Zhang et al. 2019) computes Comptonized spectra in Kerr space-time, providing the emergent spectrum from a disk–corona system consisting of a disk blackbody and a Comptonized power law from an extended corona. To generate the final spectrum, the disk and coronal properties must be specified. For the disk, MONK assumes a thin Novikov–Thorne (NT) profile (Novikov & Thorne 1973; Page & Thorne 1974), requiring input parameters such as the black hole mass and spin, the mass accretion rate, and the radial extent of the disk. These determine the blackbody photon distribution emitted from different annuli. For the corona, one must provide the optical depth, electron temperature, and geometric configuration. Although MONK accurately models disk–corona spectra in a fully relativistic framework, it does not currently incorporate any feedback from the corona to the disk. However, based on earlier studies of two-phase accretion flows, a self-consistent treatment of disk and corona emission is expected to emerge naturally when such feedback is included. This is the gap our work aims to address.

2.1. Steps in MONK

Here we briefly describe the procedure by which MONK computes the final spectra, as well as the iteration scheme used to achieve equilibrium between the disk and corona. The steps are as follows:

(i) For each radial annulus of the disk, null geodesics are calculated in Kerr geometry to track photon propagation. A photon emitted from the disk can follow one of four possible trajectories: reaching to the observer at infinity, reaching to the corona, hitting back to the disk (self-irradiation) and lost within the horizon of the BH.

(ii) Photons reaching the corona undergo Compton scattering with hot electrons, governed by the Klein–Nishina cross section. The number of scatterings and the associated energy exchange are determined by the specified coronal properties. For further details on steps (i) and (ii), see Zhang et al. (2019). The subsequent steps describe the upgrades to MONK implemented in this work to incorporate disk–corona interactions and to achieve an equilibrium solution via iterative procedures.

(iii) Comptonized photons from the corona can either escape to infinity, illuminate the disk, or fall into the black hole. An albedo parameter is introduced to specify the fraction of the illuminating flux that is reflected by the disk; the remaining fraction is absorbed and added to the intrinsic disk dissipation. Reflection is modeled using Chandrasekhar’s multiple scattering formalism, assuming the disk acts as a semi-infinite electron atmosphere.

(iv) As are coronal photons, the self-irradiated photons are reflected and/or absorbed by the disk depending on the albedo.

(v) The reflected photons have the same four possibilities as the original blackbody photons. If they reach to the corona, they are further Comptonized. We note that here we ignore the second and higher order reflection or absorption of the photons after first reflection.

(vi) All the above steps lead to the total luminosity from the disk-corona system L_tot considering the photons available at infinity and the photons lost in the BH,

$$\begin{array}{cc}\hfill L_{\mathrm{tot}}=L_{\mathrm{disk},\inf }& +L_{\mathrm{corona},\inf }+L_{\mathrm{refl}-\mathrm{selfirr},\inf }\hfill \\ \hfill & +L_{\mathrm{refl}-\mathrm{corona},\inf }+L_{\mathrm{disk},\mathrm{BH}}+L_{\mathrm{corona},\mathrm{BH}}.\hfill \end{array}$$(1)

Different luminosity components are categorized depending on the different possible sources and destinations of the photons. The first word in the subscript (i.e. disk, corona, refl-selfirr, refl-corona) represents the source of the photon, whereas the last word (i.e., inf, BH) represents the final destination of the photon. For example, L_{refl − selfirr, inf} is the luminosity reaching to infinity from the reflection of self-irradiated photons, and L_corona, BH is the luminosity lost due to coronal photons ending in the BH.

(vii) Absorption of the illuminating photons by the disk increases its effective temperature, thereby modifying the disk blackbody photon distribution that serves as the primary source of seed photons for Comptonization. This, in turn, alters the coronal luminosity. Additionally, Comptonization of reflected photons further modifies the total luminosity, L_tot. Therefore, we iterate the above steps until L_tot converges, i.e., changes by less than 1% between two consecutive iterations. This converged value represents the total luminosity of the disk–corona system in equilibrium. Hereafter, whenever we refer to L_tot, we mean this final, converged luminosity.

(viii) If the coronal electrons are hotter compared to the seed blackbody photons on average, the corona contributes to the luminosity significantly and L_tot becomes larger than the total accretion power available,

$$\begin{array}{c}\hfill L_{\mathrm{acc}}=\eta \dot{M}{c}^2,\end{array}$$(2)

where η is the accretion efficiency governed by the spin of the BH, Ṁ is the mass accretion rate. If this situation lasts a while, the accretion energy should be dissipated in the corona to keep the electrons hot, otherwise the corona will cool down almost instantly. Therefore, to achieve the balance, we introduce a fraction α, which splits the viscous power generated by accretion into part that is radiated away and part (1 − α)L_acc that is used to energize the corona (L_intr and L_sc in Appendix A). With this modification, we note that the total available accretion power remains the same (L_acc in Equation 2); the factor α only scales the intrinsic disk dissipation and consequently the disk blackbody photons distribution. Therefore the fundamental condition for the iteration of α becomes L_tot = L_acc. We make no assumption about the mechanism that transfers (1 − α) fraction of the disk viscous energy to the corona to heat up electrons.

(ix) Typically, we start the computation with α = 1. α remains fixed until convergence is achieved in L_tot through iterating the MONK computations for the given input parameters. After reaching to the converged L_tot we modify the α so that

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{new}}=\frac{L_{\mathrm{acc}}}{L_{\mathrm{tot}}}{\alpha }_{\mathrm{old}},\end{array}$$(3)

where α_new is the value for the next run and α_old is the value from the last run. Once α is modified, it changes the seed photon distribution for Comptonization and as a consequence will change the L_tot. Therefore, we again iterate to reach the convergence of L_tot keeping α the same. And the cycle repeats. Basically iteration of α is done to reach to the equilibrium solution where each iteration step consists of iterations on L_tot. The nonsimultaneity of iterations for α and L_tot keeps the iteration scheme stable. Equation (3) ensures that as the L_tot approaches to L_acc, the modification factor for α approaches to unity.

(x) Finally, we terminate the computation when convergence is achieved in both L_tot and α simultaneously–that is, when the changes in both quantities lie within 1% between two consecutive iterations. This ensures that the total luminosity from the disk–corona system equals the available accretion power, L_acc, with a fraction α (uniform across all radii) dissipated in the disk and the remaining fraction dissipated in the corona. The value of α is determined self-consistently and depends on the specified disk and coronal properties.

(xi) It is important to note that the above method does not always guarantee a steady-state equilibrium solution. If the corona is too hot and/or too optically thick, the feedback loop between the disk and corona becomes unstable, and no self-consistent equilibrium can be reached. In such cases, the computation is terminated once the iterated value of α drops below 10⁻¹⁰.

The luminosity computation is relatively insensitive to the grid size and the number of photons injected per geodesic. Therefore, we perform the iterative steps using a low-resolution setup, which significantly reduces computational time and resource usage. Once the equilibrium is achieved, we increase the resolution in the final step to obtain the detailed spectra.

2.2. Modeling the reflection

MONK uses a superphoton scheme to do the whole computation (see Sect. 2 in Zhang et al. 2019 for more details). Each superphoton represents a bundle of identical photons which is parameterized by: energy at infinity (E_∞), weight (w), polarization degree (δ), four-position (x^μ), wavevector (k^μ), and polarization vector (f^μ). The weight (w) signifies the photon generation rate per unit time in a distant observer’s frame. Therefore, in MONK, when we say photon, we basically mean a superphoton. The reflection from the semi-infinite electron atmosphere of the disk is diffuse in nature which is spanned in all directions above the disk. Therefore, for each incident photon one reflected photon is produced in a direction chosen randomly from uniform distribution of cos θ (within 0 and 1) and ϕ (within 0 and 2π) where θ and ϕ are the polar and azimuthal angle in spherical polar coordinates in the local rest frame above the disk. The energy of the reflected photon becomes

$$\begin{array}{c}\hfill E_{\mathrm{refl}}=\frac{p_{\mu }^{\mathrm{inc}}{u}^{\mu }}{p_{\mu }^{\mathrm{refl}}{u}^{\mu }}{E}_{\mathrm{inc}}.\end{array}$$(4)

Here, E_refl and E_inc denote the energies of the reflected and incident photons, respectively, as measured by an observer at infinity. The four-momenta of the incident and reflected photons in Boyer-Lindquist coordinates are represented by p_μ^inc and p_μ^refl, and u^μ is the four-velocity of the disk particle at the location where the incident photon strikes the disk. The weight of the reflected photon is calculated using Chandrasekhar’s multiple scattering formula for a semi-infinite electron atmosphere (Equation 163, Section 70.3, Chapter X in Chandrasekhar 1950). This scattering matrix depends on the polar angles of the incident and reflected directions, as well as the difference in their azimuthal angles. Chandrasekhar’s formulation ensures conservation of total energy flux between the incident and reflected radiation. To incorporate partial reflection, we multiply the computed weight by the disk albedo (ranging from 0 to 1), which specifies the fraction of incident flux that is reflected. The remaining fraction (1 − albedo) is absorbed and added to the intrinsic viscous dissipation of the disk. In addition to determining the weight, Chandrasekhar’s scattering matrix also provides the polarization of the reflected photon. We validate our implementation of self-irradiation reflection by benchmarking it against the results of Taverna et al. (2020), which in turn was benchmarked against the earlier work of Schnittman & Krolik (2009). A comparison for a spin parameter a = 0.998 and an inclination of 75° is presented in Fig. 1.

Fig. 1. Benchmarking the self-irradiated reflection computed using MONK, assuming the disk to be a semi-infinite electron atmosphere, against the results of Taverna et al. (2020) for a Kerr black hole with spin parameter a = 0.998 and an inclination angle of 75°. The data points from Taverna et al. (2020) were extracted using the WebPlotDigitizer..

To improve photon statistics at high energies, photons are split into multiple photons with appropriately scaled weights during scattering in the optically thin corona (Dolence et al. 2009; Zhang et al. 2019). These higher-energy photons then illuminate the disk and contribute to the coronal reflection component. Through iterative coupling between the disk and the corona, we also account for the Comptonization of reflected photons that re-enter the corona. However, to prevent an exponential increase in the total number of photons during successive iterations, photon splitting during Comptonization is applied only to the seed blackbody photons emitted from the disk.

2.3. Absorption of illumination

The value of the albedo depends on the ionization state of the disk. Even for highly ionized disks, a significant portion of the incident radiation is absorbed and reprocessed into lower-energy emission (Ross & Fabian 2005; García et al. 2013). Typically, the albedo is expected to be 0.5 or lower. Consequently, absorption of the illuminating flux plays a crucial role in heating the disk and maintaining the overall energy balance between the disk and the corona. To compute the illuminated flux per unit area in the local corotating frame of the disk, we transform the area element and proper time from the distant observer’s frame into the local frame, following the same prescription as in Section 2.1 of Zhang et al. (2024). This procedure is effectively the inverse of that used for calculating blackbody photon emission from the disk and transforming it to the observer’s frame, as described in Zhang et al. (2019). The illumination for each radial annulus is computed by including contributions from both the corona and self-irradiation. As noted in Section 2.1, we scale the effective mass accretion rate, and hence the energy flux from intrinsic disk dissipation, by the parameter α. Therefore, the effective temperature of the NT disk at a radial distance r from the BH is

$$\begin{array}{c}\hfill \sigma T_{\mathrm{eff}}^4(r)=\alpha F_{\mathrm{NT}}(r)+(1-\mathrm{albedo})F_{\mathrm{ill}}(r).\end{array}$$(5)

Here, F_NT(r) denotes the flux due to intrinsic viscous dissipation at radius r, determined by the black hole mass and mass accretion rate following Novikov & Thorne (1973), while F_ill(r) is the flux incident on that annulus from illumination. We assume that both the energy dissipation fraction α and the albedo are uniform across the disk. After each iteration of α, seed blackbody photons are generated using the updated effective temperature, following the standard prescription in MONK. The disk atmosphere is modeled as a semi-infinite planar atmosphere (Chandrasekhar 1950), which MONK uses to compute the polarization properties and angular distribution of blackbody emission (see Section 2.5 of Zhang et al. 2019 for details). From Equation (5), it is evident that even if the intrinsic dissipation component decreases due to a reduction in α (as required for equilibrium in a given disk-corona system), the total effective temperature–and thus the disk’s blackbody emission–can still remain significant, depending on the amount of illumination and the assumed albedo.

2.4. Photon index of the Comptonized radiation

In MONK, the corona is characterized by its geometry, optical depth, and electron temperature. However, from an observational standpoint, the photon index (Γ) of the Comptonized power law component is a more directly relevant quantity. To extract this, we fit the Comptonized photon flux reaching infinity using a simple power law model, N(E) = N₀, E^−Γ, where N(E) is the photon flux per unit energy. The fitting is performed using the curve_fit function from SciPy. A nontrivial aspect of this procedure is the determination of the appropriate energy range for fitting, which can vary with the input parameters. Since the power law shape emerges from the superposition of contributions from different orders of Compton scatterings (Rybicki & Lightman 1979), the fitting range must lie between the average energy of singly scattered photons and that of the highest significant scattering order. This ensures that we capture the regime where multiple orders of scattering overlap to form a power law. MONK provides information about the scattering order and energy distribution, which we use to define this fitting range. Given the low optical depths considered (τ < 1), the number of photons decreases rapidly with increasing scattering order, making the statistics poor at higher orders. To avoid this, we limit the analysis to those scattering orders that account for at least 0.1% of the total number of scattered photons. This criterion ensures reliable statistical sampling without losing physically important contributions. With increasing optical depth and keeping rest of the parameters the same, the energy range for fitting becomes broader and broader. After finding the energy range, we interpolate the scattered radiation of the MONK output in linear energy grid with bin size of 0.25 keV to do the fitting. We note that the estimation of the energy range for fitting does not depend on the inclination. After finding the energy range, we fit the scattered photons at different inclinations to infer about the photon index.

3. Results

Although MONK is capable of performing computations for a wide range of extended coronal geometries – such as spherical, slab, conical, wedge, outflowing slab, and truncated configurations – here we focus solely on the classical case of a two-phase accretion flow (Haardt & Maraschi 1991): a static, thin slab corona that fully covers the accretion disk and corotates with it. The primary motivation is to benchmark the current developments and facilitate direct comparison with results from existing literature. In a follow-up study, we will extend this analysis to explore the effects of alternative coronal geometries.

3.1. Reference values of input parameters in MONK

As reference values of the input parameters, we assume the mass of the BH, M = 10 M_⊙, spin, a = 0, color correction factor, f_col = 1.8, albedo = 0.5, Ṁ = 0.1 Ṁ_Edd (Ṁ_Edd = 2.45 × 10^{18} g/s) where Ṁ_Edd = L_Edd/ηc², L_Edd is the Eddington luminosity. The equatorial razor thin NT disk is extended from the innermost stable circular orbit (ISCO) to 500 R_g. The slab corona, corotating with the Keplerian disk, has vertical extension from 0 to 0.1 R_g and lie above the disk with the same radial extension as the disk. We fix the electron temperature of the corona to be 120 keV (higher boundary of the electron temperature typically observed in hard state, Poutanen et al. 2018). In addition to the above input parameters, we need to fix the optical depth of the corona to find the equilibrium solution. In Sect. 3.2.1 we fix it to a specific value, whereas in later sections we find the equilibrium solutions with increasing optical depths till we reach α ∼ 0. In Sect. 3.2.3 when we study the effect of different input parameters, except the investigated parameter, we fix all the other input parameters to their reference values.

3.2. Equilibrium solution with global energy balance

Here we present the results of the equilibrium solutions for a given disk–corona system. Equilibrium is achieved by iteratively adjusting the disk power fraction α that splits the viscous power between the disk and corona as described in Section 2. The fundamental condition for reaching equilibrium is that the total radiative luminosity from the system matches the total available accretion power, i.e., L_tot = L_acc. Unless stated otherwise, all results presented below correspond to the final equilibrium solutions obtained after iterating over α. The only exception is the α = 1 case discussed in Section 3.2.1, shown as the dashed line in Fig. 2.

Fig. 2. Change in the observed spectrum at an inclination of 60° due to the imposition of global energy equilibrium (Ltot = Lacc). The coronal optical depth is fixed at 0.15, and all other parameters are set to their reference values as described in Section 3.1. Different spectral components are shown using distinct colors. The solid lines correspond to the equilibrium solution obtained through iteration of the disk dissipation fraction α, which converges to a value of 0.07. The dashed lines represent the case with α = 1, where no iteration is performed and the equilibrium condition Ltot = Lacc is not enforced. The dotted black line shows the total spectrum for the α = 1 case, but renormalized such that Ltot = Lacc. The solid magenta line represents a power law fit to the scattered spectrum of the equilibrium solution, yielding a photon index of 1.82. The two vertical lines mark the energy range of 2.5−43.7 keV over which the power law fit is performed.

3.2.1. Change in spectra due to equilibrium

Along with the reference set of parameters mentioned in Sect. 3.1, here we also fix the optical depth of the corona τ = 0.15 and study the effect of imposing equilibrium on the disk-corona system. First we do not iterate α and fix its value to unity. For α = 1 case disk luminosity only due to intrinsic viscous dissipation becomes L_acc (Equation 2). In addition corona adds luminosity through inverse Comptonization and L_tot (Equation 1 after reaching convergence in luminosity as described in step (vii) of Sect. 2.1) becomes much larger than L_acc. This is shown by the dashed lines in Fig. 2. Different colors represent different spectral components. To impose consistency of L_tot = L_acc we iterate α (Equation 3) and reach to the final equilibrium solution between disk and corona, spectral components of which are shown as the solid lines in Fig. 2. The final α for equilibrium decreases to 0.07 after iterations. The reduction in α to reach an equilibrium solution indicates that the luminosity contributed by the inverse Comptonization of seed photons by the hot coronal electrons is compensated by the reduction in luminosity radiated by the disk due to intrinsic dissipation. Basically (1 − α) fraction of L_acc is energizing the corona to keep it hot by balancing its cooling through inverse Comptonization. By dotted black line we also present the total spectrum for α = 1 case (no iteration on α) but renormalized to L_tot = L_acc. The difference between the solid and dotted black lines show that the imposition of L_tot = L_acc not only changes the luminosity but also the spectra. For the disk-corona in equilibrium with α = 0.07 (α = 1) the fraction of disk in the total luminosity at 60° inclination becomes 0.22 (0.24). If we consider the energy range 2−20 keV then the disk fraction becomes only 0.06 (0.23).

We fit the scattered radiation with a power law following the method as described in Sect. 2.4 which is shown as the solid magenta line in Fig. 2. The energy range for fitting becomes 2.5−43.7 keV (4.4−58.1 keV for α = 1) which shifted toward lower energy for equilibrium solution because the seed blackbody from the disk itself shifted to lower energy due to low value of α. The two vertical dashed lines indicate the energy range for the fitting. The fitted photon index at 60° inclination becomes 1.82. However, it is noticeable in Fig. 2 that a single power law does not fit the higher-energy side very well although the energy range for fitting is computed appropriately where multiple scattering orders overlap with each other as described in Sect. 2.4. We find this fact typical for other values of input parameters as well. The deviation of Comptonized radiation from a simple power law is well known as the anisotropy break when seed photons for Comptonization are coming from one side of the corona instead of injected in the mid-plane (Fig. 1 in Petrucci et al. 2000, Sect. 3.3 in this draft). In addition to that, general relativistic effects may also contribute to the deviation.

3.2.2. Limit on Comptonization

To study the effect of optical depth – and to determine the hardest spectra achievable for a given disk–corona system – we fix all input parameters of the disk and corona to their reference values (as described in Section 3.1), varying only the coronal optical depth. For each value of optical depth, we compute the corresponding equilibrium solution. Each point in Fig. 3 represents one such equilibrium, with the different panels illustrating how optical depth influences various output parameters. The accretion efficiency (η) and the total mass accretion rate onto the disk–corona system (Ṁ) are held constant throughout, ensuring the total available accretion power (L_acc) remains unchanged. As optical depth increases, the hot corona contributes more significantly to the total luminosity (L_tot) due to enhanced Compton scattering. To maintain the energy balance required by the global equilibrium condition L_tot = L_acc (see Section 2.1), this increase must be compensated by a decrease in other radiative components. This is achieved via a reduction in the disk dissipation fraction α, found through our iterative scheme. A lower α leads to a reduced disk thermal radiation, thereby decreasing the seed blackbody photon input for Comptonization and counterbalancing the increased coronal output at higher optical depths. This behavior is shown in Fig. 3a, which demonstrates a gradual decrease in α with increasing τ. However, if the corona becomes too optically thick, no physically meaningful value of α (i.e., 0 ≥ α ≥ 1) can restore equilibrium. In such cases, the system would require α < 0 to satisfy L_tot = L_acc, which is unphysical. As a result, we find equilibrium solutions only for optical depths τ ≲ 0.16, corresponding to α ≳ 0. This constraint imposes an upper limit on the allowed optical depth (τ_max), and consequently a lower limit on the photon index (Γ_min), for any given disk–corona configuration.

Fig. 3. Equilibrium solutions for increasing coronal optical depth, computed using the reference set of parameters. Panel (a): Variation of the disk dissipation fraction α with optical depth, illustrating a reduction in intrinsic disk dissipation as optical depth increases. Panel (b): Fitted photon index as a function of optical depth for different viewing inclinations. Panel (c): Fractional contributions of the disk, corona, and reflection components to the total luminosity (integrated over all energies) at an inclination of 60°.

The variation of the fitted photon index with increasing optical depth at different inclinations is shown in Fig. 3b. As optical depth increases, the coronal spectrum becomes progressively harder due to a greater number of scatterings. Additionally, the global equilibrium condition causes the seed blackbody spectrum for Comptonization to shift to lower energies, owing to the reduction in α with increasing optical depth. At very low optical depths, the number of scattered photons is small, leading to poor statistics and unreliable fitting. Moreover, with only one or two scatterings contributing significantly, the resulting spectrum exhibits a bumpy structure, making the fitted photon index for τ = 0.001 an outlier. For an electron temperature of 120 keV and an albedo of 0.5, a static slab corona above the disk produces a Comptonized power law with photon index Γ ≳ 1.7. We find that at higher optical depths (τ > 0.11), the Comptonized spectrum becomes harder at lower inclinations, consistent with previous studies (Table 1 of Haardt & Maraschi 1993, Section 4 of Zhang et al. 2019). In contrast, for lower optical depths (0.01 ≤ τ ≤ 0.11), the trend reverses: the spectrum becomes softer at lower inclinations. This change in behavior with inclination likely arises from the interplay between photons traveling vertically and those moving sideways through the corona. In our setup, the slab corona lies just above the disk, and each radial point of the disk emits blackbody photons isotropically in its local frame (Equation 26 in Zhang et al. 2019). The specified optical depth refers to the vertical direction. Consequently, photons emitted at higher inclinations (i.e., at more oblique angles to the disk normal) traverse a longer path through the corona and thus experience a greater effective optical depth than those moving vertically. At low optical depths, sideways-traveling photons dominate the scattered radiation, leading to harder spectra at higher inclinations. As optical depth increases, vertically emitted photons begin to scatter more significantly, reversing the trend. This shift in the dominant direction of seed photon scattering is responsible for the inclination-dependent transition in the spectral hardness behavior.

For a given disk–corona configuration, the global equilibrium condition not only determines the specific values of α and Γ, but also uniquely sets the fractional luminosity contributions from the disk, corona, and reflection components. Figure 3c illustrates how these relative contributions vary with optical depth at an inclination of 60°, integrated over the full energy band. As the optical depth increases, the disk’s contribution to the total luminosity decreases, compensating for the growing emission from the corona. Interestingly, while the coronal luminosity continues to rise with increasing optical depth, the reflection component shows an early steep increase followed by saturation. This behavior arises because at high optical depths a significant fraction of the reflected photons are unable to escape the corona and instead undergo additional scatterings along with the disk blackbody photons. In this computation, we adopt an albedo of 0.5, which leads to nearly equal contributions from the disk and reflection components at τ = τ_max = 0.16. At this point, α ∼ 0, implying that intrinsic dissipation within the disk is nearly negligible. As a result, all the disk emission originates from the reprocessing of coronal illumination rather than from internal viscous heating.

3.2.3. Effect of BH spin, coronal temperature and albedo

In addition to the a = 0 case presented before as reference set in Sect. 3.2.2, we find the equilibrium solution with a = 0.9 and 0.998 keeping all the other parameters same as in Sect. 3.1. For different spins, accretion efficiency becomes different which makes the mass accretion rate Ṁ = 0.1 Ṁ_Edd different in physical units but keeps the available accretion power same (L_acc = 6.2 × 10³⁷ erg/s). Fig. 4a represents the variation of photon index at 60° inclination with the BH spin. At the same optical depth of the corona, we obtain harder spectra with higher spin (ΔΓ ∼ −0.15 for Δa = 0.998). Although for all the spins total accretion power remains the same, due to the larger area of the disk (as the inner boundary of the NT disk, ISCO, moves closer to the BH for higher spin) and the increased light bending effect, the disk illumination eventually becomes larger. Larger illumination increases the reprocessed component of the blackbody entering to the corona. This results to larger coronal output as well. For higher spin, the luminosity lost in the BH from both disk and corona increases. This increment in different luminosity components (Equation 1) for higher spin forces α to be lower (Equation 3) to bring the disk-corona in equilibrium (L_tot = L_acc). This reduction in α decreases the contribution from seed blackbody which makes the spectra harder with higher spin. The difference in photon indexes for different spins is larger at higher inclinations (not shown in the plots). If we increase the coronal temperature, then the spectra also become harder. Fig. 4b shows the variation of photon index at 60° inclination with electron temperature. Larger electron temperature transfers larger energy through inverse Comptonization, making the coronal output higher and α lower for equilibrium solution, subsequently the spectra harder. However, although higher spin and electron temperature produces relatively harder spectra at the same optical depth, only little improvement is achieved in lowering the Γ_min. The trend of hardening the spectra with higher BH spin remains the same even when the coronal temperature is estimated from GRMHD simulation, resulting larger temperature for higher spin (Kinch et al. 2021).

Fig. 4. Effect of different parameters on the spectral hardness of equilibrium solutions at an inclination of 60°. In each panel, a single parameter is varied while keeping all others fixed at their reference values, as described in Section 3.1. The reference parameter set (black hole spin a = 0, coronal electron temperature Te = 120 keV, and disk albedo = 0.5) is shown by black circles in all panels. Panel (a), (b), and (c) shows the effect of black hole spin, coronal electron temperature, and disk albedo respectively. In general, higher spin, higher electron temperature, and higher albedo lead to harder spectra at a given optical depth. However, among these, only increasing the albedo results in a significant reduction of the minimum photon index Γmin.

Figure 4c shows how the photon index (Γ) of the Comptonized spectrum at 60° inclination varies with albedo. In addition to the previously discussed case of albedo = 0.5, we compute spectra for albedo = 0 (complete absorption) and albedo = 1 (complete reflection), keeping all other parameters fixed to their reference values (Section 3.1). Consistent with earlier studies (e.g., Fig. 4 in Malzac et al. 2005, Fig. 3 in Poutanen et al. 2018), we find that higher albedo leads to systematically harder spectra. For albedo = 1, we halt the increase in optical depth at τ = 0.6, where the resulting spectrum reaches Γ ∼ 1.0. The dependence of spectral hardness on albedo arises primarily from its effect on the reprocessed blackbody radiation from the disk, rather than from intrinsic dissipation. A higher albedo reduces the amount of absorbed and reprocessed flux, thereby lowering the soft seed photon population for Comptonization. This leads to harder spectra at a fixed optical depth and significantly lowers the minimum achievable photon index (Γ_min). In the case of high optical depth (τ > τ_max), the increase in reprocessed blackbody flux–driven by stronger coronal illumination–can no longer be offset by reducing α, as it has already reached zero. This prevents the system from attaining equilibrium. However, when albedo is increased, the reprocessed blackbody component is suppressed, allowing equilibrium to be reached even at higher τ values. For albedo = 1, the reprocessed blackbody contribution vanishes entirely, effectively decoupling the corona from the disk. In this scenario, the seed blackbody photons arise solely from intrinsic dissipation (Equation 5). At large optical depths, the corresponding low α values yield minimal seed photon flux, resulting in only a modest increase in coronal luminosity with increasing τ. This reduced feedback allows a much larger τ_max and, consequently, a much lower Γ_min.

3.3. Escaping and anisotropic fraction

The equilibrium solution for a given disk–corona system is achieved through iteration of the α parameter, which represents the fraction of the total available accretion power dissipated in the disk. Once equilibrium is established, the final converged value of α not only fixes the intrinsic disk and coronal luminosities, but also determines the distribution of all radiative components. In particular, it allows us to compute two important diagnostic quantities: the escaping fraction (f_e) and the anisotropic fraction (f_a). The escaping fraction f_e quantifies the fraction of the total disk luminosity that escapes directly to infinity without undergoing any scattering. The anisotropic fraction f_a is defined as the fraction of scattered coronal luminosity that irradiates the disk. The mathematical definitions of both quantities are provided in Appendix A. Since the MONK computation tracks each individual photon contributing to the global disk–corona equilibrium, we are able to determine the luminosity in each component with high accuracy. From these, we compute f_e and f_a by taking appropriate ratios of the relevant luminosities after the equilibrium solution has been obtained.

For a slab corona with vertical optical depth τ, it is instructive to compare the actual escaping fraction f_e with the naive expectation exp(−τ). We find that f_e is systematically lower than exp(−τ) due to geometric and angular effects. Specifically, photons do not all propagate vertically; those emitted at oblique angles traverse a larger effective optical depth and are thus more likely to be scattered. Moreover, the angular distribution of the emitted photons is not isotropic–limb darkening and the cos θ dependence of the emission introduce directional weighting, further reducing f_e relative to the vertical approximation. The anisotropic fraction f_a also reflects the inherent geometry of the system. Since the seed blackbody photons are emitted from the bottom boundary of the slab corona (i.e., the disk surface) rather than from its mid-plane, the setup naturally introduces an asymmetry in the radiation field. Photons that are scattered backward toward the disk typically undergo more energetic scatterings compared to those scattered forward toward infinity. As a result, a larger fraction of the scattered luminosity tends to illuminate the disk rather than escape, as also highlighted in Haardt (1993). This anisotropy plays a crucial role in shaping the balance of energy between the disk and the corona.

For the equilibrium solutions obtained using the reference parameter set described in Section 3.1, we present the variation of the escaping fraction (f_e) and anisotropic fraction (f_a) as functions of optical depth τ in Fig. 5. To assess the deviation of f_e from analytical expectations, we also plot exp(−2τ). We note that the factor of 2 arises because in MONK, τ is defined as the vertical optical depth from the mid-plane to the boundary of the slab. As anticipated, f_e falls below this analytical estimate due to angular dependence and geometric effects, as discussed earlier. The behavior of f_a with τ is more complex. At low optical depths (τ ≤ 0.025), the illumination of the disk is dominated by singly scattered photons, which tend to be backward-scattered and hence strongly anisotropic. As τ increases in this regime, more photons are scattered once and redirected toward the disk, leading to an increase in f_a. However, at higher optical depths (τ > 0.025), multiple scatterings become significant. These multiple scattered photons are more isotropic in their angular distribution, which slightly reduces f_a. This trend reflects the reduction of anisotropy with increasing scatterings, a phenomenon that also manifests as the so-called ‘anisotropy break’ in the Comptonized spectrum (Petrucci et al. 2000). Although we only show results up to τ = 0.16 – beyond which no equilibrium solution with L_tot = L_acc is possible for our setup – the qualitative behavior of f_a continues to evolve at larger optical depths. For τ ≳ 5 (with the exact threshold depending on the electron temperature), f_a begins to increase monotonically, eventually approaching unity, as expected for a slab that fully traps and reprocesses the radiation.

Fig. 5. Variation of the escaping fraction (fe) and anisotropic fraction (fa) with the coronal optical depth (τ). The equilibrium solutions are computed using the same reference parameter set as in Fig. 3. To facilitate comparison, the theoretical expectation for fe, given by exp(−2τ), is also shown. The factor of 2 arises because τ in MONK is defined as the vertical optical depth from the mid-plane of the slab to its surface.

3.4. Indication from local imbalance

In our computations, we consider an extended slab corona with uniform physical properties – specifically, constant electron temperature and optical depth across all radial locations – positioned just above the accretion disk. With this configuration, we impose a condition of global energy equilibrium, requiring that the total luminosity emitted by the disk-corona system matches the total available accretion power. However, we do not impose any condition of local equilibrium, i.e., energy balance between the disk and corona at each radius. To examine the consistency of the global equilibrium solution, we compare the radial profile of coronal cooling with the available local heating power in the corona. The radial distribution of coronal cooling is computed by tracking the energy transferred from electrons to photons during each scattering event within the corona, recording both the scattering location and the energy exchange. This information is then vertically integrated over each radial annulus – defined identically to the disk annuli – to obtain the total cooling rate per radial bin. The available heating power for the corona at each radius is estimated as (1 − α) times the local accretion energy release given by the Novikov-Thorne (NT) prescription. That is, the fraction (1 − α) of the NT flux is assumed to heat the corona locally, while the remaining α contributes to the intrinsic disk dissipation. As per the NT model, the luminosity drops sharply at the innermost stable circular orbit (ISCO) due to the zero-torque boundary condition (Novikov & Thorne 1973; Zimmerman et al. 2005), and this feature is naturally incorporated in our computation of local accretion power.

The radial variation of the power available in the corona and the corresponding coronal cooling rate is shown in Fig. 6 for τ = 0.15, with all other parameters fixed to the reference set (Section 3.1). Global energy balance ensures that the total area under both curves equals the available coronal power, i.e., (1 − α)×L_acc = 5.6 × 10³⁷ erg s⁻¹, where the equilibrium solution yields α = 0.07 (see Section 3.2.1, Fig. 2). We observe that in most of the inner region of the slab corona, the cooling rate lies below the locally available power, implying that a portion of the dissipated energy remains stored in the corona. In contrast, within a narrow innermost zone and across the entire outer region, the coronal cooling rate exceeds the locally available power. While the discrepancy in the innermost region is relatively small–and sensitive to the radial binning scheme–the excess cooling in the outer regions is substantial. This mismatch can be primarily attributed to our simplifying assumption of a uniform electron temperature and optical depth across the slab. Physically, we expect the electron temperature to decrease with increasing radius, as the corona becomes less compact and less energetic away from the black hole. Such a radial decline in temperature would naturally reduce the cooling rate in the outer regions, potentially bringing it in line with the locally available power. Therefore, to achieve not only global but also local energy balance between the disk and the corona, either (i) the electron temperature must decrease with radius, or (ii) there must be a mechanism–possibly magnetic energy transport–that redistributes energy from the inner to the outer corona to sustain the high temperatures required there. In future work, we will explore radial profiles of decreasing electron temperature to examine whether such configurations can produce self-consistent solutions satisfying both global and local energy balance.

Fig. 6. Comparison of the coronal cooling rate with the available accretion power in the corona as a function of radius, for a coronal optical depth of τ = 0.15. All other parameters are set to their reference values.

4. Discussion

4.1. Static slab corona (possibly) can not explain hard state observations

Assuming a disk albedo of 0.5 and a typical coronal electron temperature, we find that the hardest photon index (Γ_min) achievable with a static slab corona fully covering the disk is approximately 1.7−1.8, depending on the black hole spin, coronal temperature, and inclination. Although the relationship between disk albedo and ionization parameter (ξ) is not straightforward, Fig. 3 of Malzac et al. (2005)–which estimates albedo based on reflected/reprocessed flux above 1 keV–suggests that an albedo of 0.5 corresponds to ξ ∼ 2000 for a constant-density slab. Thus, the Γ_min values obtained in our study would be expected only when ξ ∼ 2000, which is consistent with earlier studies of two-phase accretion flows in local energy balance (Haardt & Maraschi 1993; Malzac et al. 2005; Poutanen et al. 2018). Observationally, typical ionization parameters for the inner disk reflection fall in the range ξ ∼ 10³ − 10⁴, supporting the plausibility of our assumed albedo (Basak et al. 2017; Steiner et al. 2017; Buisson et al. 2019; Chakraborty et al. 2020; Zdziarski et al. 2021; Chand et al. 2024). Even much more self-consistent treatment of disk-corona equilibrium incorporating GRMHD simulation, non-uniform temperature and optical depth of the corona, and energy dependent disk albedo with proper photoionization treatment can produce the hardest power law only ∼1.6 (Liu et al. 2025). However, for most X-ray binaries, the hardest observed photon indices lie in the range Γ ∼ 1.5 − 2.0 (Wu & Gu 2008; Yang et al. 2015; Liu et al. 2023; You et al. 2023), with some requiring Γ < 1.5 (Chakraborty et al. 2020; Datta et al. 2024; Chand et al. 2024). Therefore, while spectra with Γ ∼ 1.6 − 1.7 can be reproduced with a static slab corona above the disk, harder spectra likely require an alternative mechanism or coronal geometry. One possible way to produce harder spectra is by increasing the disk albedo (and hence its ionization), as shown in Fig. 4c. However, this comes at the cost of significantly reducing or even washing out the iron line feature in the reflection spectrum. This might be acceptable if the observed iron line in the hard state is primarily produced by outer disk reflection, while the inner disk reflection contributes mainly to the continuum. Such an interpretation is supported by the top panel of Fig. 5 in Zdziarski et al. (2021) and Fig. 3 in Datta et al. (2024). This scenario is also consistent with the observed weakening of the iron line in the hard state, possibly due to reflection passing through a Comptonizing corona before reaching the observer (Petrucci et al. 2001; Steiner et al. 2016). Conversely, if a single reflection component suffices, or if both inner and outer reflection components are required to reproduce the iron line feature – as suggested by Liu et al. (2023) – then there is no flexibility to arbitrarily increase the disk ionization or albedo. In such cases, the static slab corona cannot explain spectra with Γ < 1.7, effectively ruling it out as a viable geometry for the hardest observed spectral states in XRBs.

The key requirement for producing harder spectra lies in suppressing the soft blackbody component. This can be achieved through either (i) a truncated disk geometry with little or no overlap between the disk and corona (Poutanen et al. 2018), and/or (ii) an outflowing corona (Malzac et al. 2001). In the case of a truncated disk, nonthermal seed photons – such as those from bremsstrahlung or synchrotron emission – become essential for Comptonization if the disk and corona do not overlap; otherwise, the resulting spectrum can become excessively hard. With the framework developed in the current version of MONK, we plan to investigate these alternative geometries in a follow-up study. Beyond spectral hardness, the fractional contribution of different emission components to the final spectrum (e.g., Fig. 3c) provides additional constraints on the disk–corona geometry and associated model parameters.

4.2. Limitations to include atomic transitions for reflection

One of the major limitations of our model is that the current reflection treatment does not include atomic transitions, which prevents us from fitting observed reflection spectra accurately. Although it is in principle possible to couple our Monte Carlo code with an external radiative transfer code to model reflection more realistically – an approach demonstrated by Malzac et al. (2005), Kinch et al. (2019) – this would significantly increase computational complexity and resource requirements, and is therefore beyond the scope of the present work. An alternative approach would be to adopt reflection spectra from pre-computed reflection tables, such as reflionx (Ross & Fabian 2005), XILLVER (García & Kallman 2010; García et al. 2013), or the more recent STOKES tables that include polarization (Podgorný et al. 2022). This methodology has been employed in previous works, including Poutanen et al. (2018) and Krawczynski & Beheshtipour (2022). However, there are several reasons that currently prevent us from incorporating these tables into our framework.

First and foremost, among currently available models, only a single reflection table–refhiden (Ross & Fabian 2007)–includes both intrinsic viscous dissipation and external power law illumination. Even this model is limited: it assumes a minimum blackbody temperature of 0.4 keV, making it suitable only for very high mass accretion rates onto stellar-mass black holes. None of the other widely used reflection tables incorporate intrinsic disk dissipation, which has been shown to significantly affect the reflection features (see Figure 6 in Ross & Fabian 2007). Additionally, these tables do not account for self-irradiation of the disk, a factor that can become important for high-spin black holes and low coronal heights (Mirzaev et al. 2024). Furthermore, when modeling reflection that includes many atomic transitions alongside scattering processes, it becomes computationally challenging to achieve a fully converged solution for the temperature and radiation field within a reasonable number of iterations. This often results in incomplete convergence and a mismatch between the incident and reflected flux, limiting the reliability of these tables to only a narrow range of input parameters (see Section 2.1 of Poutanen et al. 2018; footnote 1 in Zdziarski & De Marco 2020; and Section 3.3 in Krawczynski & Beheshtipour 2022). A final and fundamental limitation arises from the specific assumptions about illumination in these reflection tables. They typically assume a power law spectrum with fixed high- and low-energy cutoffs as the incident radiation. To use such tables in our framework, we would need to approximate the incident radiation at each radial annulus of the disk as a power law with a cutoff. This poses a significant computational challenge: for an optically thin corona, the number of photons diminishes rapidly with the number of scatterings, making it particularly difficult to characterize the high-energy cutoff accurately. Since reflection depends not only on photon energy but also on incidence and emission angles, we would need to model angular distributions of the incoming radiation, dramatically increasing the required photon statistics. Therefore, while incorporating reflection spectra with atomic transitions remains an important future direction, it demands substantial methodological and computational developments beyond the scope of the present study.

4.3. Effect of self-irradiation without corona

If we consider only the disk and its self-irradiation and no corona, then we expect that for equilibrium solution the total luminosity collected at infinity (L_tot in Equation 1) to be L_acc, i.e., α should be equal to 1 as there is no additional luminosity contributed by the corona. However, it does not become strictly true always. As the self-irradiation becomes more significant with higher spin, α deviates more from 1. The reasons for this discrepancy are two fold. First one is the ignorance of second order reflection (reflection of reflected photons) in our computation which becomes more relevant at higher albedo. Second one is the ignorance of work done by the self-irradiation on the disk. We keep the angular momentum profile of the disk as prescribed by Novikov & Thorne (1973) in Kerr geometry even with self-irradiation. Due to this imposition, there can be slight gain or loss of energy of a photon due to reflection with respect to an observer at infinity. Because the gravitational red-shift of a photon depends not only on the location of the disk but also on the direction of the photon propagation. Therefore, due to diffuse nature of reflection and reprocessed radiation, the average direction of photon after reflection or absorption of self-irradiation may change. This effect may lead to a change in energy flux with respect to an observer at infinity even though Thomson scattering through Chandrasekhar’s semi-infinite electron atmosphere conserves energy flux in the comoving frame of the disk. To model the effect of self-irradiation precisely, we need to solve the hydrodynamics of the Keplerian disk incorporating the angular momentum transport by the self-irradiation as incorporated in KERRBB model (Li et al. 2005). However, this effect is always quite small with respect to the flux irradiated to the disk at every radii (Fig. 1 in Li et al. 2005). With no corona, due to self-irradiation we find that α can differ from 1 by 5% in the worst case scenario. The effect of self-irradiation is relevant only when either there is no corona or it is very optically thin. For our main focus of this work, with disk-corona interaction, we can safely ignore this effect as coronal dissipation plays a major role in achieving the equilibrium.

4.4. Comparison with methodology of HM93 (Haardt & Maraschi 1993)

For a two-phase accretion flow (disk and corona) in equilibrium, we present in Appendix A a set of analytical estimates for various luminosity components in the Newtonian framework, following the approach of HM93. However, there is a fundamental difference between our formalism and HM93. We introduce the escaping fraction parameter f_e, and assume that the coronal luminosity is the sum of the disk luminosity which interacts with the hot electrons and the fraction of accretion power available to energize the hot electrons (Equation A.3). On the other hand, HM93 assumes that the coronal luminosity is only the accretion power available in the corona. These two assumptions converge with the approximation of f_e = 1, leading to the same analytical estimations of luminosities in different radiative components through our formalism and HM93 as shown at the end of Appendix A. From Fig. 1b of HM93 it is clear that the hardest photon index achievable for two phase accretion flow in equilibrium is for the optically thinner corona which has higher electron temperature. For very optically thin corona, the escaping fraction of luminosity f_e will become very close to unity and nullify any difference between our presented analytical estimations and HM93.

In our numerical framework, equilibrium solutions are obtained by iterating the parameter α within the MONK code under the constraint L_tot = L_acc. From these simulations, we obtain the fractions of various luminosity components observed at infinity (e.g., as shown in Fig. 3c for an inclination of 60°). We then use the analytic expressions from Appendix A, inputting the corresponding values of α, f_e, f_a, and the assumed albedo, to independently compute the fractional luminosities. These analytic predictions agree very well with the values obtained directly from MONK, particularly for low black hole spin, where general relativistic corrections are minimal and thus not required in the Newtonian estimates.

It is worth noting that the computation of f_e and f_a from the MONK output involves only the luminosities entering and escaping the corona, as well as the scattered luminosity from the corona illuminating the disk. The condition used in the iterative scheme, L_tot = L_acc, does not explicitly involve the disk emission luminosity (e.g., as defined in Equation A.5). Therefore, the good agreement between the full set of luminosity components derived analytically and those obtained from MONK validates the self-consistent iteration scheme used in our approach. The primary purpose of the analytic estimates is to demonstrate that, although our numerical method reaches the equilibrium solution by iterating on the total luminosity L_tot alone (step (ix) in Section 2), the resulting solution naturally aligns with the analytical predictions for individual luminosity components, similar to the framework in HM93. In our case, inverse Comptonization and the additional power provided by the hot corona are handled entirely by MONK, and the equilibrium is found without requiring the analytic luminosity expressions. In contrast, HM93 modeled Comptonization using an amplification factor and relied on analytic luminosity estimates to relate this factor to the spectral slope and the fraction of energy dissipated in the corona.

4.5. Similarities and differences with other models

Dovčiak et al. (2022) developed a physical model, kynsed, to study the global energy balance between the disk and the corona, enabling broadband SED fitting of AGNs within xspec. In this model, self-consistency is achieved by transferring a fraction of the accretion energy to a point-like, lamp-post corona, in accordance with the observed power law component. The disk is then heated via illumination by this corona. In contrast to our approach – where a uniform fraction of accretion energy is transferred from the disk to the corona at all radii – kynsed assumes that only the inner region of the disk, within a specified transition radius, transfers its entire energy to the corona. Reflection is modeled using the xillver tables (García et al. 2013), and the absorbed flux is estimated by subtracting the reflected flux (above 0.1 keV) from the incident illumination. Another key difference lies in the treatment of Comptonization. In kynsed, the corona is assumed to be a point-like lamp-post located along the disk axis at a certain height, and its emission is modeled as a power law with a given normalization – no actual Comptonization is computed. After fitting, the model estimates the size and optical depth of the corona by assuming a spherical geometry. In contrast, MONK performs full Comptonization from an extended corona and, as shown in Zhang et al. (2024), this significantly alters the disk illumination profile compared to that of a lamp-post geometry. However, MONK does not currently support direct spectral fitting of observational data unless a spectral table is generated from its output.

Another contemporary model, kerrC, developed by Krawczynski & Beheshtipour (2022), performs ray tracing in the Kerr metric and models Comptonization by an extended corona – conceptually similar to MONK. However, while it includes reflection from the disk, it does not perform any iterative feedback between the disk and corona and thus does not achieve a final equilibrium solution. A key distinction lies in how reflection is modeled. kerrC implements reflection using Chandrasekhar’s multiple scattering formalism, as done in our work. However, it then rescales the weight of each photon based on the xillver high-density reflection tables (García et al. 2016), under the assumption that the reflection spectrum from a highly ionized slab matches that of a pure electron atmosphere. This approximation enables kerrC to estimate the polarization of reflected emission using the xillver tables, which themselves do not include polarization information. However, recent studies on reflection with polarization suggest that emission from atomic transitions leads to significant depolarization, even in highly ionized slabs (Podgorný et al. 2022). This implies that the scaling approach adopted by kerrC may not be appropriate for accurately modeling polarized reflection.

It should be noted that with the current upgrade of the MONK code, it becomes capable to find equilibrium solution for a given disk-corona system incorporating exact prescription of Comptonization by extended corona in complete relativistic framework. This is novel in nature and no other available code does the similar computation as far as we know. Therefore, although the previous version of MONK code (Zhang et al. 2019) is benchmarked with kynsed (Fig. 6 in Dovčiak et al. 2022) and kerrC (Section 3.5 in Krawczynski & Beheshtipour 2022), it is not possible to compare the results from these codes with the equilibrium solutions achieved through the upgraded version of MONK.

4.6. Possibilities to improve

In this first attempt to find an equilibrium solution with the disk in global balance with an extended corona, several simplifying assumptions have been made. Most notably, we assume a homogeneous corona with uniform optical depth and electron temperature across all radii. However, in reality, the electron temperature is expected to decrease with increasing distance from the black hole. While the density should also decrease outward, the accompanying increase in scale height makes it nontrivial to predict how the coronal optical depth varies with radius. Another simplification is the heating of the disk in addition to the Novikov–Thorne (NT) profile, without solving the disk’s hydrodynamics. This is clearly an idealization. A more realistic treatment would involve using radial and vertical profiles of temperature and optical depth derived from analytical models or simulations of two-phase accretion flows (Cho & Narayan 2022; Bambic et al. 2024), though such models are still limited in availability. We also assume, somewhat ad hoc, that the same fraction of accretion energy is available to the corona at all radii, which may not reflect the actual physical scenario. In Section 3.4, we discuss how enforcing local equilibrium between the corona and the underlying disk at each radius may offer insights into the radial variation of α or the inhomogeneous nature of the corona. Additionally, for an extended corona above the disk, a typical emissivity index of 3 implies that the disk ionization should decrease with radius. This would, in turn, reduce the albedo at larger radii – an effect not incorporated in our current model. The challenge of including energy-dependent albedo is discussed in Section 4.2. All of these simplifying assumptions highlight areas where the current model can be improved to better reflect physical reality.

5. Conclusion

In this work, we further upgraded the existing MONK code to model reflection by treating the disk as a semi-infinite electron atmosphere. The portion of the incident flux that is not reflected is absorbed by the disk, thereby increasing its temperature above the purely Keplerian value. The mutual feedback between the extended corona and the disk enabled us to determine an equilibrium solution for a given disk–corona system within a fully general relativistic framework. Our key findings are summarized below:

• A static slab corona covering the disk cannot account for the hard state observations of XRBs, particularly for photon indices Γ < 1.7, consistent with earlier studies based on local disk–corona equilibrium.

• For a slab coronal geometry in equilibrium with the disk, higher black hole spin and higher electron temperature lead to harder spectra; however, they do not significantly reduce Γ_min, the lowest photon index achievable from a disk–corona system in equilibrium.

• As in previous studies, Γ_min can be significantly reduced by increasing the disk albedo (i.e., its ionization); however, this requires fine-tuning to remain consistent with observations (see Section 4.1).

• Although the disk and corona are in global equilibrium, achieving local balance at each radius requires either an outward flow of energy within the corona or a radial decrease in coronal temperature to reduce the energy transferred through Comptonization (see Section 3.4).

• Further work is needed to model reflection more accurately, including detailed atomic transitions. Other plausible geometries, such as a truncated disk or an outflowing corona, have the potential to produce harder spectra and will be explored in future studies.

Acknowledgments

SRD would like to thank Roberto Taverna and Andrzej Zdziarski for the discussions about reflection and Comptonization. SRD, MB, MD 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: Theoretical estimations of luminosities

If the disk-corona system is in a global equilibrium then the disk and coronal luminosities do not remain independent any more. Here we present the estimation of luminosities expected from different components when the disk-corona is in equilibrium, very similar to what presented in Haardt & Maraschi (1991), HM93. We assume that the α fraction of the available accretion power (L_acc) is available for intrinsic disk dissipation (blackbody radiation) and the rest (1-α) is dissipated in the corona to re-energize hot electrons that scatter the photons. Therefore,

$$\begin{array}{c}\hfill L_{\mathrm{intr}}=\alpha L_{\mathrm{acc}},L_{\mathrm{sc}}=(1-\alpha )L_{\mathrm{acc}}.\end{array}$$(A.1)

Further, the albedo fraction of the illuminated luminosity reflected back from the disk gives

$$\begin{array}{c}\hfill L_{\mathrm{refl}}=\mathrm{albedo}\times L_{\mathrm{ill}}.\end{array}$$(A.2)

Here we assume the anisotropy fraction f_a of the coronal luminosity goes downward and illuminates the disk (η in HM93). For isotropic coronal radiation, f_a = 0.5. In addition, we assume the escaping fraction through corona is f_e. Therefore, f_e fraction of the luminosity from the disk (blackbody emission + reflection) which enters the corona from below escapes through, and rest contributes to the total coronal luminosity. Therefore, total coronal luminosity (L_corona) and the illumination to the disk becomes

$$\begin{array}{c}\hfill L_{\mathrm{corona}}=(1-f_{\mathrm{e}})(L_{\mathrm{disk}}+L_{\mathrm{refl}})+L_{\mathrm{sc}},\end{array}$$(A.3)

and

$$\begin{array}{c}\hfill L_{\mathrm{ill}}=f_{\mathrm{a}}{L}_{\mathrm{corona}}.\end{array}$$(A.4)

The above steps finally lead to the luminosity emitted by the disk, which comes from intrinsic dissipation as well as absorption of illumination:

$$\begin{array}{cc}\hfill L_{\mathrm{disk}}& =\alpha L_{\mathrm{acc}}+(1-\mathrm{albedo})L_{\mathrm{ill}},\hfill \\ \hfill & \Rightarrow L_{\mathrm{disk}}=\frac{\alpha +f_{\mathrm{a}}-\alpha f_{\mathrm{a}}-\mathrm{albedo}(f_{\mathrm{a}}-\alpha f_{\mathrm{a}}{f}_{\mathrm{e}})}{1-f_{\mathrm{a}}+f_{\mathrm{a}}{f}_{\mathrm{e}}}{L}_{\mathrm{acc}}.\hfill \end{array}$$(A.5)

As the disk is completely covered by the slab corona for our case, the observed disk luminosity is

$$\begin{array}{c}\hfill L_{\mathrm{disk},\mathrm{obs}}=f_{\mathrm{e}}{L}_{\mathrm{disk}}.\end{array}$$(A.6)

Similarly, we can reduce the expressions of luminosity from other components as well in terms of α, f_e, f_a and albedo. Summing up all the luminosity components, the total luminosity reaching to the observer at infinity from the disk-corona system is

$$\begin{array}{c}\hfill L_{\mathrm{tot}}=f_{\mathrm{e}}(L_{\mathrm{disk}}+L_{\mathrm{refl}})+(1-f_{\mathrm{a}})L_{\mathrm{corona}}=L_{\mathrm{acc}},\end{array}$$(A.7)

which is the basic assumption for the total disk-corona system to be in global energy balance. For a given disk-corona system (including specifying the albedo) equilibrium between disk and corona provides unique values of α, f_e and f_a. We note that at equilibrium all these three parameters are coupled. Using these values, we can estimate the observed disk luminosity using equation A.6 corresponding to the given disk-corona system. Similarly we also can estimate the observed luminosity from corona as well as due to reflection. We note that the L_tot of analytic estimation in equation A.7 is the same as in equation 1 with the difference of only the luminosity components contributed by photons from self-irradiation or lost in the BH as these components come in the MONK computations solely due to general relativistic effects which is ignored in the simple analytic estimation.

The numerical scheme of achieving equilibrium solution by iterating α (as described in Sect. 2.1) fixes the f_e, f_a self-consistently. In the limit of very very optically thin corona (τ → 0), equilibrium solution gives α → 1 (Fig. 3(a)), f_e → 1 (Fig. 5). Putting these values in equation A.3–A.6 we find that irrespective of actual value of f_a and albedo, L_corona = L_sc = 0 and L_disk, obs = L_disk = L_acc. This gives the expected result that if the corona ceases to exist, the total accretion power is dissipated in the disk and we observe only the disk component.

The analytical estimations presented above is very similar to what presented in HM93 to find the equilibrium solution. However, although the assumption of disk luminosity is the same in both the works (equation 2a in HM93), there is a fundamental difference about the assumption of coronal luminosity. Compared to HM93, we introduced the parameter f_e and assume the L_corona to be the sum of the available accretion power in the corona and the luminosity from the disk interacting with the hot electrons (equation A.3). On the other hand, HM93 assumed that the coronal luminosity is only the accretion power available in the corona, L_corona, HM = L_sc (equation 3b in HM93). We can mimic the convergence of our formalism with HM93 by assuming f_e = 1 without putting any constraint on α.

Now, if we approximate f_e = 1, we find L_corona = L_sc = (1 − α)L_acc, and L_disk = (α+f_a−αf_a−albedo(f_a−αf_a))L_acc. Therefore, with the limit of f_e = 1, the ratio between the two becomes

$$\begin{array}{cc}\hfill \frac{L_{\mathrm{corona}}}{L_{\mathrm{disk}}}& =\frac{1-\alpha }{\alpha +f_{\mathrm{a}}-\alpha f_{\mathrm{a}}-\mathrm{albedo}(f_{\mathrm{a}}-\alpha f_{\mathrm{a}})}\hfill \\ \hfill & =\frac{1-\alpha }{1-(1-\alpha )+f_{\mathrm{a}}(1-\alpha )-\mathrm{albedo}\times f_{\mathrm{a}}(1-\alpha )}\hfill \end{array}$$(A.8)

Following the notation used in HM93, if we replace (1-α) by f, albedo by a, and f_a by η, we get back the same equation (5a) in HM93.

It should be noted that in our formalism of equilibrium solution for a given disk-corona system, f_e is coupled with α. Here to show the convergence of our formalism to HM93 by putting f_e = 1 without putting any constraint on α simply indicates that we assume the same as in HM93, L_corona consists only the accretion power available in the corona (L_sc). Once we make our fundamental assumptions same, we reach to the same equations of luminosities of HM93 as expected.

All Figures

	Fig. 1. Benchmarking the self-irradiated reflection computed using MONK, assuming the disk to be a semi-infinite electron atmosphere, against the results of Taverna et al. (2020) for a Kerr black hole with spin parameter a = 0.998 and an inclination angle of 75°. The data points from Taverna et al. (2020) were extracted using the WebPlotDigitizer..
In the text

Fig. 2. Change in the observed spectrum at an inclination of 60° due to the imposition of global energy equilibrium (Ltot = Lacc). The coronal optical depth is fixed at 0.15, and all other parameters are set to their reference values as described in Section 3.1. Different spectral components are shown using distinct colors. The solid lines correspond to the equilibrium solution obtained through iteration of the disk dissipation fraction α, which converges to a value of 0.07. The dashed lines represent the case with α = 1, where no iteration is performed and the equilibrium condition Ltot = Lacc is not enforced. The dotted black line shows the total spectrum for the α = 1 case, but renormalized such that Ltot = Lacc. The solid magenta line represents a power law fit to the scattered spectrum of the equilibrium solution, yielding a photon index of 1.82. The two vertical lines mark the energy range of 2.5−43.7 keV over which the power law fit is performed.

In the text

Fig. 3. Equilibrium solutions for increasing coronal optical depth, computed using the reference set of parameters. Panel (a): Variation of the disk dissipation fraction α with optical depth, illustrating a reduction in intrinsic disk dissipation as optical depth increases. Panel (b): Fitted photon index as a function of optical depth for different viewing inclinations. Panel (c): Fractional contributions of the disk, corona, and reflection components to the total luminosity (integrated over all energies) at an inclination of 60°.

In the text

Fig. 4. Effect of different parameters on the spectral hardness of equilibrium solutions at an inclination of 60°. In each panel, a single parameter is varied while keeping all others fixed at their reference values, as described in Section 3.1. The reference parameter set (black hole spin a = 0, coronal electron temperature Te = 120 keV, and disk albedo = 0.5) is shown by black circles in all panels. Panel (a), (b), and (c) shows the effect of black hole spin, coronal electron temperature, and disk albedo respectively. In general, higher spin, higher electron temperature, and higher albedo lead to harder spectra at a given optical depth. However, among these, only increasing the albedo results in a significant reduction of the minimum photon index Γmin.

In the text

Fig. 5. Variation of the escaping fraction (fe) and anisotropic fraction (fa) with the coronal optical depth (τ). The equilibrium solutions are computed using the same reference parameter set as in Fig. 3. To facilitate comparison, the theoretical expectation for fe, given by exp(−2τ), is also shown. The factor of 2 arises because τ in MONK is defined as the vertical optical depth from the mid-plane of the slab to its surface.

In the text

	Fig. 6. Comparison of the coronal cooling rate with the available accretion power in the corona as a function of radius, for a coronal optical depth of τ = 0.15. All other parameters are set to their reference values.
In the text