#  $R$-matrix analysis of 22Ne($\alpha$, n)25Mg and 22Ne($\alpha$,
$\gamma$)26Mg reaction

Rajkumar Santra  Department of Nuclear and Atomic Physics, Tata Institute of
Fundamental Research, Mumbai-400005, India.

###### Abstract

The 22Ne($\alpha$, n)25Mg and its competing channel 22Ne($\alpha$,
$\gamma$)26Mg has an major influence on neutron flux in weak s-process
nucleosynthesis path in low mass AGB stars and massive stars of mass (M$\geq$
10M⊙). So the ratio rate of this two competing reaction control the neutron
flux in weak s-process nucleosynthesis. Various experiment has been performed
to study the properties of nuclear states of 26Mg to evaluate rate of
22Ne+$\alpha$ reaction rate and corresponding rate from these studies vary by
up to a factor of 500 in the astrophysical relevant temperature. The recent
evaluation by Philip et al. of 22Ne($\alpha$, n)25Mg reaction rate using most
recent nuclear data of 26Mg from number of sources shows similar result with
previous estimation for 22Ne($\alpha$, $\gamma$)26Mg but got lower rate for
22Ne($\alpha$, n)25Mg reaction due to updated nuclear data. Also Philip et al.
suggested that rate based on full $R$-matrix modeling will required to take
into accounted the interference effects between distant levels and sub-
threshold resonance. In present work full $R$-matrix calculation has been
performed for 22Ne($\alpha$, n)25Mg and 22Ne($\alpha$, $\gamma$)26Mg reaction
based on fitting the 22Ne($\alpha$, n)25Mg reaction data of Jaeger et al; in
energy range 0.8 to 1.45 MeV and updated nuclear data of 26Mg states. The
$R$-matrix fitting for the 22Ne($\alpha$, n)25Mg reaction nicely explain
experimental data in 0.8 to 1.45 MeV energy range by changing spin, parity of
Ex = 11.784 and 11.63 MeV states from 1- to 0+.

###### keywords:

The resonance capture, $R$-matrix analysis

Almost half of the elements heavier than iron are synthesis via slow (over
time scales of thousands of years) neutron capture reaction on stable isotope
in s-process nucleosynthesis path [1]. The s-process manly activated in
Asymptotic Giant Branch (AGB) stars, seeded by 56Fe iron. The 13C($\alpha$,
n)16O and 22Ne($\alpha$, n)25Mg reaction are main source of neutrons for
s-process. For low mass AGB stars with $1M_{\odot}\leq M\leq 4M_{\odot}$ both
this two reaction activated He inter-shell during the Thermal Pulses
contribute as a smaller neutron source and responsible for synthesis elements
of Atomic masses A $\approx$ 90 - 209. But for massive stars ($M\geq
8M_{\odot}$) only 22Ne($\alpha$, n)25Mg reaction is main source of neutron
which synthesis isotopes of mass A $\approx$ 60 - 90 is so call weak s-process
nucleosynthesis. Also +ev Q-value reaction 22Ne($\alpha$, $\gamma$)26Mg is act
as a competing reaction that decide neutron release from 22Ne($\alpha$, n)25Mg
reaction. So constraining the rate of this two reaction has high sensitive for
weak s-process nucleosynthesis.

So many direct [2] and [4, 3, 5] indirect experimental measurement has been
performed to extracted the spin, parity, partial widths
($\Gamma_{\alpha},\Gamma_{n},\Gamma_{\gamma}$) to constraining the rate of
22Ne($\alpha$, n)25Mg and 22Ne($\alpha$, $\gamma$)26Mg reactions. Recently
rate of this two reaction has been re-evaluate [6] based upon updated nuclear
data from a number of sources and they suggested that an $R$-matrix modeling
will required due to lack of uncertainty of spin, parity of so many relevant
states of 26Mg as well as interference effects same partials waves between two
states.

In this context an multilevel R-matrix analysis has been performed on
available cross-section data of 22Ne($\alpha$, n)25Mg reaction including
interference effect in energy range E${}_{\alpha}^{c.m}$ $\approx$ 0.8 to 1.45
MeV range and extrapolate up to 0.57 MeV energy. Main aim is to constraining
spin-parity and study interference effect of states of 26Mg in excitation
energy Ex $\approx$ 11.319 to 11.828 MeV energy range.

##  1 R-matrix analysis of 22Ne($\alpha$, n)25Mg reaction

The behavior of excitation function for 22Ne($\alpha$, n)25Mg reaction was
mainly determined by the resonance capture process through the several
resonance states of compound nucleus 26Mg. In this work an $R$-matrix
calculation has been performed to describe capture data and constrain the
spin, parity and width of the states of 26Mg.

Figure 1:  $R$-matrix calculation with literature reported resonance parameters. Filled symbols represent the direct measurement data taken from Ref. [2] Table 1: Summary of resonance parameters used in $R$-matrix calculation with literature reported values [6] for comparison with direct measurement data of 22Ne($\alpha$, n)25Mg reaction. Ex | Er | Jπ |  $\Gamma_{\alpha}$(eV) |  $\Gamma_{n}$(keV)
---|---|---|---|---
(MeV) | (MeV) | Literature [6] | Literature[6] | Literature[6]
11.828 | 1.214 | 2+ | 0.18 | 1.1
11.7847 | 1.169 | 1- | 8.0 $\times$ 10-3 | 24.5
11.749 | 1.1456 | 1- | 0.02 | 64
11.63 | 1.016 | 1- | 2.4$\times$10-4 | 13.5
11.526 | 0.9112 | 1- | 4.3$\times$10-4 | 1.8
11.508 | 0.894 | 1- | 1.2$\times$10-4 | 1.27
11.461 | 0.847 | 3- | 7.9$\times$10-6 | 6.55
11.441 | 0.827 | 3- | 5.5$\times$10-6 | 1.47
11.3196 | 0.7056 | 1- | 5.5$\times$10-6 | 0.132
11.272 | 0.65 | 3- | 9.2$\times$10-8 | 1.81
11.258 | 0.644 | 2+ | 1.0$\times$10-6 | 0.41
11.171 | 0.557 | 2+ | 1.9$\times$10-8 | 0.03
11.169 | 0.552 | 3- | 4.4$\times$10-10 | 1.94
11.163 | 0.549 | 2+ | 2.7$\times$10-9 | 5.31
11.112 | 0.498 | 2+ | 4.3$\times$10-10 | 2.095
11.084 | 0.470 | 2+ | 5.7$\times$10-11 | -
10.9491 | 0.3351 | 1- | 3.0$\times$10-14 | 30

$R$-matrix modeling of 22Ne($\alpha$, n)25Mg reaction has been performed using
AZURE2 code [7]. This code was developed based on the theory developed of Lane
and Thomas [8] and of Vogt [9]. In $R$-matrix formalism, radial space is
divided into two disting regions –an internal region extended up to a radius
R${}_{c}\approx r_{0}(A_{p}^{1/3}+A_{t}^{1/3})$ and an external region above
Rc. A choice of radius for entrance and exit channels is needed for the model
calculation. Rc of the two channels have been obtained through $\chi^{2}$
minimization. However, as channel radius is not a free parameter in the model,
we performed a grid search on the channel radius by changing the value in
small steps and varying the parameters to get the fit. The chosen channel
radii values are 5.37 $fm$ for the 22Ne + $\alpha$ channel and 4.21 $fm$ for
the 25Mg + n channel. During calculation in AZURE2, the energy resolution of
the system also accounted.

Table 2: Summary of resonance state parameters of 26Mg used for $R$-matrix calculation obtained from Litarature and $R$-matrix fit of 22Ne($\alpha$, n)25Mg data with Eα=0.8 to 1.3 MeV. Ex | Er | Jπ |  $\Gamma_{\alpha}$(eV) |  $\Gamma_{n}$(keV)
---|---|---|---|---
(MeV) | (MeV) |  |  |
10.9491 | 0.335 | 3- | 3.0$\times$10-14[4] | 0.125
11.084 | 0.470 | 2+ | 5.7$\times$10-11[4] | 0.175
11.112 | 0.498 | 2+ | 4.3$\times$10-10[13] | 2.095 [14]
11.163 | 0.549 | 2+ | 2.7$\times$10-9[13] | 5.31[14]
11.169 | 0.552 | 3- | 4.4$\times$10-10[13] | 1.94[14]
11.171 | 0.557 | 2+ | 1.9$\times$10-8[13] | 0.01[14]
11.272 | 0.65 | 2+ | 1.0$\times$10-6[6] | 0.410[14]
11.278 | 0.71 | 3- | 9.2$\times$10-8[6] | 1.81[14]
11.319 | 0.705 | 1- | 5.1429$\times$10-5 | 0.452
11.4401 | 0.827 | 3- | 4.59$\times$10-6 | 0.700
11.458 | 0.847 | 3- | 9.0014$\times$10-6 | 15.588
11.506 | 0.894 | 1- | 1.399$\times$10-4 | 15.346
11.525 | 0.911 | 1- | 2.529$\times$10-4 | 0.5209
11.63 | 1.016 | 0+ | 9.106$\times$10-3 | 14.341
11.759 | 1.145 | 1- | 0.0502 | 139.935
11.784 | 1.169 | 0+ | 23.952 $\times$ 10-3 | 17.226
11.8276 | 1.214 | 2+ | 0.210657 | 1.144
Table 3: Summary of estimated ANC from literature reported spectroscopic factors of 26Mg bound states Ex | l | Jπ | S | cb | C
---|---|---|---|---|---
(MeV) |  |  |  | fm-1/2 | fm-1/2
g.s | 0 | $0^{+}$ | 0.418[3] | 4.57$\times$10+4 | 2.95$\times$10+4
1.81 | 2 | 2+ | 0.0915[11] | 2.35$\times$10+4 | 7.10$\times$10+3
2.94 | 2 | 2+ | 0.0211[11] | 2.16$\times$10+4 | 3.13$\times$10+3
3.59 | 0 | 0+ | 0.291[11] | 3.61$\times$10+4 | 1.94$\times$10+4
4.32 | 4 | 4+ | 0.046[11] | 3.16$\times$10+4 | 6.78$\times$10+2
4.33 | 2 | 2+ | 0.12[3] | 2.03$\times$10+4 | 7.032$\times$10+3
4.84 | 2 | 2+ | 0.11[3] | 2.02$\times$10+4 | 6.69$\times$10+3
4.9 | 4 | 4+ | 0.2[3] | 3.07$\times$10+4 | 1.37$\times$10+3
4.97 | 0 | 0+ | 0.281[11] | 3.62$\times$10+4 | 1.91$\times$10+4
5.47 | 4 | 4+ | 0.05[3] | 3.04$\times$10+4 | 6.79$\times$10+2
5.71 | 4 | 4+ | 0.04[3] | 9.80$\times$10+4 | 1.96$\times$10+3
6.745 | 2 | 2+ | 0.22[3] | 4.17$\times$10+4 | 1.95$\times$10+4
6.876 | 3 | 3- | 0.17[3] | 2.81$\times$10+4 | 1.15$\times$10+4
7.348 | 3 | 3- | 0.5[3] | 3.14$\times$10+4 | 2.22$\times$10+4
8.036 | 2 | 2+ | 0.2[3] | 6.93$\times$10+4 | 3.09$\times$10+4
8.937 | 2 | 2+ | 0.14[3] | 1.91$\times$10+5 | 7.14$\times$10+4
9.325 | 2 | 2+ | 0.38[3] | 4.73$\times$10+5 | 2.91$\times$10+5
9.371 | 4 | 4+ | 0.38[3] | 9.2$\times$10+4 | 5.67$\times$10+4
9.856 | 2 | 2+ | 0.1[3] | 6.81$\times$10+6 | 2.153$\times$10+6
10.573 | 1 | 1- | 1 | 1.92$\times$10+33 | 1.92$\times$10+33
Table 4: Summary of bound and resonance state parameters of 26Mg used for $R$-matrix calculation obtained from Literature and $R$-matrix fit of Jaeger et. al., data 22Ne($\alpha$, n)25Mg data with Eα=0.8 to 1.3 MeV. Ex | Er | Jπ |  $\Gamma_{\alpha}$(eV)/ANC(fm-1/2) |  $\Gamma_{n}$(keV)/ANC(fm-1/2) |  $\Gamma^{E^{f}}_{\gamma}$(eV)
---|---|---|---|---|---
(MeV) | (MeV) |  |  |  |
g.s | - | $2^{+}$ | 2.95$\times$10+4 | - |
1.81 | - | 2+ | 7.10$\times$10+3 | - |
2.94 | 0.470 | 2+ | 3.13$\times$10+3 | - |
3.59 | 0.470 | 0+ | 1.94$\times$10+4 | - |
4.32 | 0.470 | 4+ | 6.78$\times$10+2 | - |
4.33 | - | 2+ | 7.032$\times$10+3 | - |
4.84 | - | 2+ | 6.69$\times$10+3 | - |
4.9 | - | 4+ | 1.37$\times$10+3 | - |
4.97 | - | 0+ | 1.91$\times$10+4 | - |
5.47 | - | 4+ | 6.79$\times$10+2 | - |
5.71 | - | 4+ | 1.96$\times$10+3 | - |
6.745 | - | 2+ | 1.95$\times$10+4 | - |
6.876 | - | 3- | 1.15$\times$10+4 | - |
7.348 | - | 3- | 2.22$\times$10+4 | - |
8.036 | - | 2+ | 3.09$\times$10+4 | - |
8.937 | - | 2+ | 7.14$\times$10+4 | - |
9.325 | - | 2+ | 2.91$\times$10+5 | - |
9.371 | - | 4+ | 5.67$\times$10+4 | - |
9.856 | - | 2+ | 2.153$\times$10+6 | - |
10.573 | - | 1- | 1.92$\times$10+33 | - | 0.094g.s [12]
|  |  |  |  | 0.1064.97[12]
10.696 | - | 4+ | 3.5$\times$10-46 [17] | - | 31.81 [17]
10.717 | - | 1- | 2.53$\times$10-36[16] | - | 3(15)(g.s) [6]
10.805 | - | 1- | 1.11$\times$10-22 [15] | - | 0.16g.s[12]
|  |  |  |  | 0.561.81[12]
10.819 | - | 0+ | 3.18$\times$10-21 [15] | - | 3g.s[6]
10.943 | - | 2+ | 5.0$\times$10-11(UL) | 30 [18] | 6.5g.s [18]
10.9491 | 0.3351 | 1- | 3.0$\times$10-14[4] | 0.125 | 0.26g.s[12]
|  |  |  |  | 1.071.81 [12]
|  |  |  |  | 0.252.93 [12]
|  |  |  |  | 0.093.58 [12]
|  |  |  |  | 0.204.33 [12]
11.084 | 0.470 | 2+ | 5.7$\times$10-11 [4] | 0.175 | 34.31 [6]
11.112 | 0.498 | 2+ | 4.3$\times$10-10[13] | 2.095 [14] | 1.71.81[14]
11.163 | 0.549 | 2+ | 2.7$\times$10-9[13] | 5.31[14] | 2.8[14]
11.169 | 0.552 | 3- | 4.4$\times$10-10[13] | 1.94[14] | 3.31.81[14]
11.171 | 0.557 | 2+ | 1.9$\times$10-8[13] | 0.01[14] | 5 5.47[14]
11.272 | 0.65 | 2+ | 1.0$\times$10-6[6] | 0.410[14] | 2.20.0[14]
11.278 | 0.71 | 3- | 9.2$\times$10-8[6] | 1.81[14] | 0.31.81[14])
11.301 | 0.71 | 2+ | 1.53$\times$10-5[6] | - | 31.81[6]
11.319 | 0.7056 | 1- | 5.1429$\times$10-5 | 0.452 | 3g.s[6]
11.3277 | 0.7056 | 1- | 1.8$\times$10-6[6] | - | 2.2g.s[14]
11.343 | 0.7056 | 2+ | 1.1$\times$10-6[6] | - | 1g.s[14]
11.4401 | 0.827 | 3- | 4.59$\times$10-6 | 0.700 | 31.81[6]
11.458 | 0.847 | 3- | 9.0014$\times$10-6 | 15.588 | 31.81[6]
11.5002 | 0.894 | 1- | 1.95$\times$10-1[6] | - | 3g.s[6]
11.506 | 0.894 | 1- | 1.399$\times$10-4 | 15.346 | 3g.s[6]
11.525 | 0.9112 | 1- | 2.529$\times$10-4 | 0.5209 | 3g.s[6]
11.63 | 1.016 | 0+ | 9.106$\times$10-3 | 14.341 | 3g.s [6]
11.759 | 1.1456 | 1- | 0.0502 | 139.935 | 3g.s [6]
11.784 | 1.169 | 0+ | 23.952 $\times$ 10-3 | 17.226 | 31.81[6]
11.8276 | 1.214 | 2+ | 0.210657 | 1.144 | 31.81[6]

In this present work, initially R-matrix calculation has been performed
22Ne($\alpha$, n)25Mg using the recently updated spin, parity and partial
widths of resonance states (Ex=10.9 to 12.82 MeV) reported in Ref. [6] for the
observed resonances in 22Ne($\alpha$, n)25Mg reaction for
E${}_{\alpha}^{c.m}$= 0.57 to 1.45 MeV range and the resonance states
parameters of 26Mg used in calculation are listed in Table 1 and results of
$R$-matrix calculation (without and with experimental resolution correction)
for cross-section of 22Ne($\alpha$, n)25Mg reaction compare with experimental
data of Ref [2] as shown in Fig. 1. With including experimental resolution
correction in $R$-matrix calculation still data are not well reproduced.

Now to better explain experimental data, we adjust resonance energy and widths
of populated all the states with excitation energy range 11.319 to 12.828 MeV
(E${}_{\alpha}^{c.m}$=0.7 to 1.45 MeV) and by adjusting this parameters we
reproduced the data form 0.7 to 0.98 MeV energy region but fails to reproduced
data 0.98 to 1.2 MeV energy region due to strong destructive interference
between same Jπ states. Now spin, parity of Ex=11.63 and 11.784 MeV resonances
are changing form 1- to 0+ and $\Gamma_{\alpha}$, $\Gamma_{n}$ left as free
parameters for all the states with excitation energy range 11.319 to 12.828
MeV. The resultant $R$-matrix fit nicely explain the experimental data in 0.8
to 1.45 MeV energy region as shown in Fig. 2. The fitted parameters are listed
in Table 4.

Figure 2:  (Color Online) $R$-matrix fit with changing spin and parity shown
in Table 2. Filled symbols represent the direct measurement data taken from
Ref. [2] Figure 3:  (Color Online) $R$-matrix calculation for 22Ne($\alpha$,
$\gamma$)26Mg capture reaction.

##  2 Estimation of ANC of bound states 26Mg for direct capture calculation of
22Ne($\alpha$, $\gamma$)26Mg reaction

The Asymptotic Normalization Coefﬁcients(ANC) define the amplitude of the tail
of bound state wave-function for two-body configurations. In low energy
radiative capture reaction, capture occurs via the tail of bound state wave
function and the ANC determined magnitude of direct capture cross-sections. It
is related to the spectroscopic factor of the two-body conﬁguration as

$C=\sqrt{S}\times c_{b}$ (1)

Where S is the spectroscopic factor of two body configurations. cb is single
particle ANC. The spectroscopic factors of different bound states of 26Mg are
taken from literature [3, 11] which was determined from 22Ne(6Li, d)26Mg
transfer reactions. The cb are calculated using the code FRESCO (ver. 3.2)
[10] code. In FRESCO calculation $\alpha$ binding potentials parameters are
use same vale for which spectroscopic factors reported in [3, 11]. The
resultant $\alpha$-single particle ANC and total ANC for different bound
states of 26Mg are listed in Table 3.

##  3 R-matrix calculation of 22Ne($\alpha$, n)25Mg and 22Ne($\alpha$,
$\gamma$)25Mg reaction

A full R-matrix calculation has been performed for 22Ne($\alpha$, n)25Mg with
E${}_{\alpha}^{c.m}$ = 0.57 to 1.45 MeV energy range and 22Ne($\alpha$,
$\gamma$)25Mg E${}_{\alpha}^{c.m}$ = 0 to 1.45 MeV to evaluate reaction rate
of this two reaction at Gamow window region. In this calculation Ex=10.69 to
11.827 MeV excitation energy all natural parity resonance states of 26Mg (0+,
1-, 3- 2+, 4+) are consider to evaluate resonance capture contribution of
22Ne($\alpha$, n)25Mg and 22Ne($\alpha$, $\gamma$)25Mg reaction. The resonance
states of Ex=11.319 to 11.827 MeV (observed in measurement of 22Ne($\alpha$,
n)25Mg by Jaeger et al. [2]) parameters(Ex, Jπ, $\Gamma_{\alpha}$,
$\Gamma_{n}$) are taken by $R$-matrix fitting Jaeger et al. data listed in
Table 2. For the case of other resonance states parameter’s (Ex, Jπ,
$\Gamma_{\alpha}$, $\Gamma_{n}$, $\Gamma_{\gamma}$) taken from updated nuclear
data of 26Mg that are got by indirect measurement or theoretical calculation.
The non resonant capture cross-section( direct capture + sub-threshold state)
in 22Ne($\alpha$, $\gamma$)26 reaction are calculated using ANC,s for
different bound states of 26Mg that are extra from 22Ne(6Li, d)26Mg reaction
listed in Table 3. The states of 26Mg and their parameters used in this
calculation are listed in Table 4. The final calculated S-factors of
22Ne($\alpha$, n)25Mg and 22Ne($\alpha$, $\gamma$)25Mg reaction are shown in
Fig. 4.

Figure 4:  $R$-matrix calculation to the S-factor for 22Ne($\alpha$, n)25Mg
22Ne($\alpha$, $\gamma$)26Mg capture reaction.

##  4 Reaction rate of 22Ne($\alpha$, n)25Mg reaction

The reaction rate of 22Ne($\alpha$, n)25Mg reaction has been reevaluated
numerically using AZURE2 code for temperature range 0.01 to 10 GK. The present
rate ratio with respect to Longland et al., [17] median rate are shown in Fig.
5

Figure 5:  The 22Ne($\alpha$, n)25Mg rate of present work relative to Longland
et al., meadian rate. Dash line corresponds to upper and lower limits.

##  5 Conclusion

A $R$-matrix analysis in 22Ne($\alpha$, n)25Mg reaction was performed to the
description of direct measurement data [2]. The present calculation with
including experimental resolution corrections with recently evaluated
resonance parameters [6] for 26Mg are not well describe measurement data [2]
due to strong destructive interference between same Jπ states and as well as
energy location of some resonances. With changing energy location, decay
widths of some resonances and spin, parity of Ex = 11.784 and 11.63 MeV states
from 1- to 0+ the experimental data are well describe for Eα 0.8 to 1.45 MeV
energy region. The extrapolation of $R$-matrix calculations with indirectly
measured resonance parameters well off with respect to highly uncertain data
of Jaeger et. al.,[2] for 0.57 to 0.8 MeV range. The cross sections for
22Ne($\alpha$, $\gamma$)26Mg are also estimated simultaneously with the
resulting parameters and taking the gamma partial widths from the literature.
The present $R$-matrix calculation shows that 22Ne($\alpha$, n)25Mg is
dominant with respect to 22Ne($\alpha$, $\gamma$)26Mg reaction over the hole
energy range from 0.57 to 1.2 MeV. It is also observed that the R-matrix
calculation yields higher cross sections for 22Ne($\alpha$, n)25Mg in the
Gamow energy window compared to the previous estimates of Philip et al.,[6].
The present reaction rate of 22Ne($\alpha$, n)25Mg reaction evaluated from
$R$-matrix calculation is higher by order an of magnitude as compare to Philip
et al., reported at 0.2 to 0.3 GK temperature.

## References

  * [1] M. Pignatari et al., Astrophys. J. 710, 1557 (2010).
  * [2] M. Jaeger et al., Phys. Rev. Lett. 87, 202501 (2001).
  * [3] S. Ota et al., Phys. Rev. C 104, 055806 (2021).
  * [4] H. Jayatissa et al., Physics Letters B 802 (2020) 135267.
  * [5] U. Giesen et al., Nuclear Physics A561 (1993) 95-111.
  * [6] Philip Adsley et al., Phys. Rev. C 103, 015805 (2021).
  * [7] R. E. Azuma, E. Uberseder, E. C. Simpson, and C. R. Brune et al., Phys. Rev. C 81, 045805 (2010).
  * [8] A. Lane and R. Thomas, Rev. Mod. Phys. 30, 257 (1958).
  * [9] E. Vogt, Rev. Mod. Phys. 34, 723 (1962).
  * [10] I. J. Thompson, Comput. Phys. Rep. 7,167 (1988).
  * [11] N. Anantaraman et al., Nuclear Physics A 279 (1977) 474-492.
  * [12] R. J. DEBOER et al. PHYSICAL REVIEW C 82, 025802 (2010).
  * [13] S. Ota et al., Physics Letters B802(2020)135256.
  * [14] C. Massimi et al., Physics Letters B 768 (2017) 1–6.
  * [15] G. Lotay et al, Eur. Phys. J. A (2019) 55: 109.
  * [16] R. TALWAR et al. PHYSICAL REVIEW C 93, 055803 (2016).
  * [17] R. LONGLAND, C. ILIADIS, AND A. I. KARAKAS PHYSICAL REVIEW C 85, 065809 (2012).
  * [18] C. MASSIMI et al. PHYSICAL REVIEW C 85, 044615 (2012).


# Observational properties of puffy disks: radiative GRMHD spectra of mildly
sub-Eddington accretion

Maciek Wielgus,1,2,3 Debora Lančová,4,5 Odele Straub,6 Włodek Kluźniak,5
Ramesh Narayan,3,2 David Abarca,5 Agata Różańska,5 Frederic Vincent,7 Gabriel
Török,4 and Marek Abramowicz4,5,8
1Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn,
Germany
2Black Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA
02138, USA
3Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street,
Cambridge, MA 02138, USA
4Research Center for Computational Physics and Data Processing; Institute of
Physics, Silesian University in Opava, Czech Republic
5Nicolaus Copernicus Astronomical Centre, Polish Academy of Sciences, Bartycka
18, 00-716 Warsaw, Poland
6Max Planck Institute for Extraterrestrial Physics, Gießenbachstraße 1,
D-85748 Garching, Germany
7LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Universités,
UPMC Univ. Paris 06, Univ. de Paris, Sorbonne Paris Cité,
5 place Jules Janssen, 92195 Meudon, France
8Department of Physics, Göteborg University, SE-412-96 Göteborg, Sweden
E-mail<EMAIL_ADDRESS>(MW)E-mail<EMAIL_ADDRESS>(DL)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Numerical general relativistic radiative magnetohydrodynamic simulations of
accretion disks around a stellar mass black hole with a luminosity above 0.5
of the Eddington value reveal their stratified, elevated vertical structure.
We refer to these thermally stable numerical solutions as puffy disks. Above a
dense and geometrically thin core of dimensionless thickness $h/r\sim 0.1$,
crudely resembling a classic thin accretion disk, a puffed-up, geometrically
thick layer of lower density and $h/r\sim 1.0$ is formed. We discuss the
observational properties of puffy disks, in particular the geometrical
obscuration of the inner disk by the elevated puffy region at higher observing
inclinations, and collimation of the radiation along the accretion disk spin
axis, which may explain the apparent super-Eddington luminosity of some X-ray
objects. We also present synthetic spectra of puffy disks, and show that they
are qualitatively similar to those of a Comptonized thin disk. We demonstrate
that the existing xspec spectral fitting models provide good fits to synthetic
observations of puffy disks, but cannot correctly recover the input black hole
spin. The puffy region remains optically thick to scattering; in its spectral
properties the puffy disk roughly resembles that of a warm corona sandwiching
the disk core. We suggest that puffy disks may correspond to X-ray binary
systems of luminosities above 0.3 of the Eddington luminosity in the
intermediate spectral states.

###### keywords:

X-rays: binaries – accretion, accretion discs – radiative transfer –
scattering – black hole physics – software: simulations

††pubyear: 2022††pagerange: Observational properties of puffy disks: radiative
GRMHD spectra of mildly sub-Eddington accretion–Observational properties of
puffy disks: radiative GRMHD spectra of mildly sub-Eddington accretion

##  1 Introduction

Our understanding of accretion onto stellar mass black holes in X-ray binary
systems has been built by five decades of observations following the Uhuru
satellite mission (Giacconi et al., 1971), and theoretical developments
starting with the $\alpha$-disk solution of Shakura & Sunyaev (1973). The
former are reviewed in, e.g., Levine et al. (1996); Remillard & McClintock
(2006) and the latter in, e.g., Done et al. (2007); Abramowicz & Fragile
(2013).

Most known black hole X-ray binaries (BHXBs) undergo a cycle of outbursts,
driven by a mass accretion rate modulation caused by an outer disk instability
(Homan et al., 2001; Lasota, 2001; Bagińska et al., 2021). During outburst, in
a luminous ($0.05<L/L_{\rm Edd}<0.3$) soft (thermal) spectral state, the
observations are consistent with multi-blackbody spectra expected from the
classic analytic thin $\alpha$-disk model in its relativistic version (NT
disk; Novikov & Thorne, 1973; Page & Thorne, 1974). However, in the radiation
pressure dominated regime $p_{\rm rad}\gg p_{\rm gas}$, expected in the inner
part of the accretion disk, the analytic model is viscously (Lightman &
Eardley, 1974) and thermally (Shakura & Sunyaev, 1976) unstable. Since the
observations indicate stability of BHXB accretion disks on relevant
timescales, one must conclude that while the thin disk model is effective in
predicting spectra consistent with observations in the applicable luminosity
regime (Li et al., 2005; McClintock et al., 2014), it is not a self-consistent
and correct description of the physical reality. For larger mass accretion
rates stability can be achieved through dynamical advection, which removes the
excess thermal energy, as found for the (vertically integrated) slim disk
models (Abramowicz et al., 1988; Sądowski, 2009), but the puzzle remains for
the observationally interesting case of mildly sub-Eddington luminosities.
Several stabilizing mechanisms have been suggested (e.g., Różańska et al.,
1999; Oda et al., 2009; Ciesielski et al., 2012; Zhu & Narayan, 2013), and
magnetic fields seem to be particularly promising.

At luminosities above $\sim 0.3\,L_{\rm Edd}$ the observed spectra of BHXBs
are no longer consistently represented by the thin disk models (McClintock et
al., 2006), as the disk geometry begins to deviate from the geometrical
thinness assumption. Attempts to utilize the more general theoretical
framework of the slim disk model do not fully address the inconsistency
(Straub et al., 2011). Unlike thin (or slim) disks, high mass accretion rate
flows may in reality include the presence of a warm corona (Zhang et al.,
2000; Gronkiewicz & Różańska, 2020), increasing the system energetic output in
soft X-ray, and/or outflows (Ohsuga et al., 2009), modifying the dynamical
structure of the system. The impact of the upper disk layers (disk atmosphere)
is typically solved for in separation from the radial disk structure (Davis &
Hubeny, 2006; Różańska et al., 2011; Sądowski et al., 2011). In numerical
simulations, the three dimensional structure of the accretion disk can be
captured self-consistently.111Apart from simulations, fully two or three
dimensional (analytical) solutions for accretion disks have been obtained only
for $\alpha$-disks, e.g., Kluzniak & Kita (2000); Regev & Gitelman (2002).

A numerical model of a magnetically stabilized, radiation pressure dominated,
high luminosity accretion disk in the general relativistic radiative
magnetohydrodynamic (GRRMHD) framework has been presented by Sądowski (2016).
We have further developed this model in the sub-Eddington case (Lančová et
al., 2019), discussing the morphology of the obtained solutions: a dense,
geometrically thin disk core is sandwiched by an elevated region of gas with
density lower by about $2$ orders of magnitude that nevertheless remains
optically thick to scattering. Because of this puffed-up region, supported
radially by Keplerian rotation and vertically predominantly by the magnetic
pressure, we refer to these solutions as _puffy disks_.

GRRMHD simulations promise an approach in which the relevant physics (magnetic
fields, outflows, cooling, turbulence, etc.) is self-consistently taken into
account, at the cost of large computational complexity (e.g., Ohsuga &
Mineshige, 2011; Sądowski et al., 2013, 2014; Sądowski et al., 2017). In this
paper we discuss the observational properties of puffy disks at luminosity of
$0.6-1.5$ $L_{\rm Edd}$, inferred from GRRMHD simulations presented in
Sądowski (2016) and Lančová et al. (2019). We discuss how these numerical
models could yield a theoretical framework for interpretation of the BHXBs
spectra beyond the $\sim 0.3\,L_{\rm Edd}$ thin disk limit.

##  2 Puffy disks

Figure 1: Structure of a puffy disk. The white dashed line in both panels
indicates the vertical thickness (density scale-height) of the dense disk core
(region C). Left: Gas temperature $T_{\rm{gas}}$ with contours of the optical
depth $\tau=1$ (photosphere, separating regions P and F), $\tau=10$ and
$\tau=100$ (white continuous lines). Right: gas density map, arrows indicate
the gas velocity in $r-\theta$ plane. The white contour denotes the position
of the photosphere ($\tau\,=\,1$) and the dashed yellow line indicates the
border between the inflowing (below) and outflowing (above) regions.

###  2.1 Numerical setup and conventions

We perform GRRMHD simulations of accretion onto a 10 $M_{\odot}$ Schwarzschild
(non-spinning) black hole using KORAL,222https://github.com/achael/koral_lite
a highly parallelized GRRMHD code (Sądowski et al., 2013, 2014; Sądowski et
al., 2017). The details of the numerical setup have been described in Sądowski
(2016). The radiative part of the simulation setup accounts for the electron
scattering which dominates the total opacity in the considered regime, as well
as free-free absorption and thermal Comptonization (following Sądowski &
Narayan, 2015). At this stage frequency-averaged (bolometric) opacities are
used, with the aim of capturing the radiation feedback on the disk structure
and dynamics. We then use the relativistic radiative postprocessing code
HEROIC (Zhu et al., 2015; Narayan et al., 2016) to compute puffy disk spectra
and images, solving the complete radiative transfer problem on the converged,
time-averaged output of GRRMHD models (gas density, temperature, velocity, and
magnetic field maps in the inner part of the simulation volume). The
contributions from the outer, unconverged region of the simulation were masked
with an extrapolation procedure (Narayan et al., 2017) employed for spherical
radii $r>25\,GM/c^{2}$. HEROIC solves for the radiation field in both
optically thick and optically thin regimes, including the scattering effects.

We consider three postprocessed GRRMHD simulation runs, parametrizing the mass
accretion rate through the black hole event horizon in Eddington units333The
Eddington luminosity is defined as $L_{\rm Edd}=4\pi c\,GMm_{\rm
H}/\sigma_{\rm T}=1.26\,\times\,10^{38}\,(M/M_{\odot})$ erg/s, where $G$ is
the gravitational constant, $c$ is the speed of light, $M$ is the mass of the
gravitating body, $m_{\rm H}$ and $M_{\odot}$ are the hydrogen atom and solar
mass, respectively, and $\sigma_{\rm T}$ is the Thompson cross-section. with
$\dot{m}=0.6,0.9,1.5$ (Sądowski, 2016; Lančová et al., 2019). We take

$\dot{m}\equiv\frac{\dot{M}}{\dot{M}_{\rm
Edd}}=\frac{0.057\dot{M}c^{2}}{{L}_{\rm Edd}}\ ,$ (1)

where we defined the Eddington mass accretion rate as

$\dot{M}_{\rm Edd}={{L}_{\rm Edd}/(\eta c^{2})}={{L}_{\rm
Edd}}/{(0.057c^{2})}\ .$ (2)

In Eq. 2, $\eta$ denotes the radiative efficiency of accretion, which we set
to the thin disk value in the Schwarzschild metric, $\eta=0.057$ (Novikov &
Thorne, 1973), as a unit convention defining $\dot{M}_{\rm Edd}$.

The actual radiative efficiency in the simulated accretion process, defined
through

$\eta={{L}/(\dot{M}c^{2})},$ (3)

can be determined by computing $\dot{M}$ through the horizon and $L$ through a
fiducial surface. Note that in presence of advection sphere-integrated
luminosity may not be equal to the integrated radiative power produced locally
in the disk, causing some additional ambiguity in the efficiency definition.
In Lančová et al. (2019) we reported $L=0.36\,L_{\rm Edd}$ for $\dot{m}=0.6$,
measured in the GRRMHD simulation. This implied a significant decrease of
radiative efficiency with respect to the NT disk, which we interpreted as a
consequence of advection, prominent in the simulation. With the postprocessed
radiative transfer we found, somewhat surprisingly, the efficiency of a sub-
Eddington puffy disk in the Schwarzschild metric to be consistent with that of
the NT disk within the measurement uncertainties. While the discrepancy may be
partly explained with the imperfections of the luminosity measurement in the
GRRMHD simulation (we integrated the radiative flux over a 25 $M$ radius
sphere in the optically thin region, ignoring contributions from larger radii
or optically thick outflows), we currently suspect that this is related to
limitations of the M1 closure scheme for the radiation tensor employed in
KORAL (Sądowski et al., 2013). Under the M1 approximation radiation propagates
in a single direction in the optically thin regime, and a large part of the
radiation released in the inner flow consequently appears to be swallowed by
the black hole in the GRRMHD simulation, whereas it may possibly escape to
infinity when a more general radiative transfer problem is solved with HEROIC.

###  2.2 Properties of puffy disks

Our numerical simulations predict a stratified structure in a mildly sub-
Eddington accretion disk. In Fig. 1 we indicate three zones of the emerging
geometry. The zone denoted by C represents the high density core of the disk,
with a density scale-height thickness of $h/r\sim 0.1$, which contains most of
the gas. It is supported above the midplane of the disk primarily by magnetic
pressure which dominates over radiation pressure, and both strongly dominate
over gas pressure (Lančová et al., 2019). The white dashed line in Fig. 1
shows the density scale height of the disk core. Zone P corresponds to the
puffy region. While being optically thick to scattering and delimited by the
photosphere at $h/r\sim 1$ (white solid line in the right panel of Fig. 1),
this region is characterized by a gas density that is about 2 orders of
magnitude lower than in the disk core. The gas in the puffy region is
vigorously accreting onto the black hole. The turbulent dissipation heats up
the puffy region and the temperature reaches several times $10^{7}$ K at the
photosphere, see Fig. 2.

In comparing various disk models it is convenient to model the temperature as
a power-law function of the cylindrical radius $R$ (Kubota & Makishima, 2004).
In the puffy disk, the radial temperature dependence varies from layer to
layer in the different zones of the accreting structure. The photospheric
temperature scales as $R^{-1.2}$, see Fig. 2. While this is significantly
steeper than for a NT disk, where $T\propto R^{-0.75}$ (Shakura & Sunyaev,
1973; Novikov & Thorne, 1973; Page & Thorne, 1974), the photosphere (surface
of last scattering) of a puffy disk is also much hotter in the inner disk
region than the NT disk surface. Moreover, the puffy disk photosphere extends
below the innermost circular stable orbit (ISCO; $R=6M$ for the Schwarzschild
spacetime) with a consistent continuation of the power law temperature
dependence, whereas NT surface temperature drops to zero at ISCO. The
temperature of the disk core, both in the equatorial plane and at the density
scale-height, grows rapidly in the plunging region ($R<R_{\rm ISCO}$) as a
consequence of vertical compression, but changes more slowly with radius for
$R>R_{\rm ISCO}$. As can be seen in Fig. 2, these temperature profiles are
lower but comparable with the NT disk equatorial temperature $T\propto
R^{-0.375}$ for $R>R_{\rm ISCO}$. At $R>15M$ optically thick outflows from the
outer layer of the puffy zone develop (see the right panel of Fig. 1). The
third zone, denoted by F in Fig. 1 is a low density, low optical depth, hot
($10^{8}-10^{9}$ K) funnel zone. The optically thin material in the funnel,
driven by collimated radiation, is outflowing above the stagnation surface
located at $z\approx r+7M$.

The morphology of the puffy disk is somewhat similar to the predictions of
Zhang et al. (2000) based on X-ray spectral modelling of BHXBs, and of
Begelman & Pringle (2007), based on arguments about active galactic nuclei
disks. Conclusions similar to ours about the vertical structure of the disk,
including the importance of magnetic pressure have been reached by Mishra et
al. (2020), based on non-relativistic simulations.

Figure 2: Radial profiles of temperature in the $\dot{m}\,=\,0.6$ puffy disk
for the equatorial plane, disk core scale-height, and the photosphere. The
photospheric temperature is reasonably well approximated as a power law of the
cylindrical radius $R$ with an exponent of $p=-1.22$. The disk core heats up
predominantly very close to the black hole, with slow temperature variations
for $R>R_{\rm ISCO}$. Analytic NT disk temperature profiles for the same
system parameters are shown for comparison. Figure 3: Puffy disk appearance
given as $L_{\Omega}$, frequency-integrated isotropic luminosity per solid
angle $\left(r_{g}/D\right)^{2}$, where $r_{g}=GM/c^{2}$ and $D$ represents
the observer’s distance to the source, using $M=10M_{\odot}$, and $D=10$ kpc.
We show puffy disk images computed for different inclinations $\theta$, for
mass accretion rates of 1.5, 0.9 and 0.6 $\dot{M}_{\rm Edd}$ (first three
columns), compared to the NT disk of 0.6 $\dot{M}_{\rm Edd}$ (last column).
The white dashed lines in the bottom row of the puffy disk images correspond
to $10^{31}$ erg s${}^{-1}(r_{g}/D)^{-2}$ contours – at the inclination of 70∘
the inner disk is no longer visible in the images. All results were computed
in the Schwarzschild metric and correspond to the
$\left(240\,r_{g}/D\right)^{2}$ field of view.

##  3 Beaming and obscuration

The presence of an optically and geometrically thick puffy region has
important consequences for the radiative output received by the observer. We
illustrate this in Fig. 3, where ray-traced bolometric (frequency-integrated)
images of puffy accretion disks (computed with HEROIC; Narayan et al., 2016)
are contrasted with NT thin disk images (computed with GYOTO; Vincent et al.,
2011) for several inclinations $\theta$. For the thin disk no vertical
structure was assumed. While there seem to be no perspectives in the
foreseeable future of observing resolved structures of accretion disks in
BHXBs444The field of view of panels in Fig. 3 corresponds to 2.4 ns of arc,
four orders of magnitude less than the highest resolution achieved by global
radio-interferometric arrays (EHT Collaboration et al., 2019)., this exercise
is very helpful in recognizing the influence of inclination and elucidating
differences between the NT thin disk and the GRRMHD puffy disk solutions. At
low inclinations ($\theta\leq 30^{\circ}$) observers see the puffy disk system
through the low optical depth funnel region. The radiation then becomes
collimated in the funnel surrounded by the geometrically and optically thick
puffed-up accretion disk in a way similar to the classic analytic thick disk
solutions (Abramowicz et al., 1978), and it is enhanced by the reflections
from the optically thick funnel walls (Sikora, 1981). These effects are self-
consistently handled by the HEROIC code. Notice that the mostly purely
geometrical beaming effect discussed here is different from relativistic
beaming (the Doppler effect related to the emitter approaching with a
relativistic velocity), as observed in jets. The latter has relatively little
importance in the case of puffy disks, since the outflowing gas inside the
funnel has a velocity of only $\sim 0.1\,c$, and a very low density.

The luminosity, $L$, is computed directly from the simulated spectra, at a
fiducial distance $D$:

$\displaystyle L$ $\displaystyle=2\pi
D^{2}\int\displaylimits_{0}^{\pi}\left(\int\displaylimits_{0}^{\infty}F_{\nu}\,{\rm
d}\nu\right)\sin\theta{\rm d}\theta$ (4)
$\displaystyle=\frac{1}{2}\int\displaylimits_{0}^{\pi}\left(\int\displaylimits_{0}^{\infty}L_{\nu}\,{\rm
d}\nu\right)\sin\theta{\rm d}\theta\ .$ (5)

The collimation of radiation in the funnel causes an increase of the so called
isotropic luminosity, $L_{\rm iso}$, that would be inferred from observations
at particular inclination $\theta$, with

$L_{\rm iso}=4\pi D^{2}F=\int\displaylimits_{\rm image}L_{\Omega}\,{\rm
d}\Omega=\int\displaylimits_{0}^{\infty}L_{\nu}\,{\rm d}\nu\ .$ (6)

At the axis ($\theta=0$), the increase corresponds to a factor of $1/b=L_{\rm
iso}/L$ relatively to the actual luminosity (e.g., King et al., 2001). In Eq.
6, $F=\int F_{\nu}{\rm d}\nu$ is the flux seen by the observer at distance
$D$. In fact, the flux $F_{\nu}$, luminosity per solid angle $L_{\Omega}$, and
luminosity per frequency $L_{\nu}$, all depend on the observer’s inclination
$\theta$, as can be seen in Fig. 3, and thus $L_{\rm iso}$ is a function of
the same angle (see Fig. 4).

In our simulations we measure a puffy disk beaming factor of $b\approx
0.25-0.30$, without a strong indication of dependence on $\dot{m}$, as the
general three-zone puffy disk structure does not change much between
$\dot{m}=0.6$ and $\dot{m}=1.5$. The lack of dependence of the beaming factor
on the mass accretion rate contradicts the $b\propto\dot{m}^{-2}$ relation
expected from King (2009), which was derived on the assumption that the disk
height increases with the mass accretion rate. We attribute this discrepancy
to the roughly constant $h/r$ ratio (regardless of the accretion rate), which
is related to advection limiting the growth of the disk thickness with
$\dot{m}$ (Lasota et al., 2016; Wielgus et al., 2016).

At inclination $\theta=10^{\circ}$, Fig. 3 shows that puffy and NT disks are
similar in appearance; however, the puffy disk is brighter overall because of
the beaming effect. Because of the low optical depth in the funnel, even the
innermost region of the accretion flow in the immediate vicinity of the event
horizon is visible (the central dark spot in the puffy disk images). This
central gap in the image is much more prominent in the NT model, in which
radiation is only emitted at radii larger than the ISCO. At larger
inclinations a different effect starts to be relevant—the puffed-up optically
thick disk starts to obscure the view of the inner part of the accretion flow.
In Fig. 3 we see that at an inclination of 50∘ the puffy disk image is already
dominated by the funnel wall opposite to the observer, with the core of the
system obscured. At 70∘ inclination the inner regions of the puffy disk are
not visible and the radiative flux decreases by about 2 orders of magnitude.
At the same time, the impact of the inclination on the NT disk appearance is
much less prominent, as the disk remains geometrically thin, supporting
neither efficient beaming nor obscuration effects.

Interestingly, the development of an elevated disk structure during the high
mass accretion rate phase, has been proposed as an interpretation of the X-ray
spectra of the large inclination source V404 Cygni (Motta et al., 2017). The
enhanced complexity of the X-ray outburst light curves of high inclination
sources was also identified as hinting at a time-variable obscuration (Narayan
& McClintock, 2005). Although low luminosity numerical solutions of stable
disk are not yet available, we imagine that a transition between a thin and a
puffy disk may occur at some accretion rate $\sim 0.3\dot{M}_{\rm Edd}$.

We further quantify the collimation / obscuration effects in Fig. 4, where we
compare the beaming pattern of NT thin disks and puffy disks for the mass
accretion rates of $\dot{m}=0.6,0.9,1.5$. In general, beaming is maximal at
low inclinations and decreases for higher viewing angles and the beaming
factor does not strongly depend on the mass accretion rate. The face-on
inclination “beaming” factor555In reality, for the thin disk the beaming
pattern is related to the usual non-isotropic emission from the photosphere,
the flux decreasing with increasing inclination angle. The unrelated effect of
Doppler beaming affects the NT spectra at large viewing angles. of the NT disk
is $b_{\rm NT}\approx 0.7$ which corresponds to an effect about $2.5$ times
less prominent than in the case of the puffy disk. Hence, an NT disk with a
mass accretion rate of 1.5 $\dot{M}_{\rm Edd}$ appears as a 2.1 $L_{\rm Edd}$
luminosity source when viewed along the axis. A puffy disk of the same mass
accretion rate, however, implies an isotropic luminosity of 4.4 $L_{\rm Edd}$.

Figure 4: Bolometric beaming pattern comparison between a puffy disk and the
NT thin disk, for mass accretion rates of $0.6,0.9,$ and $1.5\,\dot{M}_{\rm
Edd}$ and spin $a_{*}=0$. Puffy disks collimate radiation along the axis of
symmetry with a beaming factor $b\approx 0.3$, and the radiation is strongly
attenuated for large inclinations.

A mildly sub-Eddington puffy disk system viewed at low inclination would be
conventionally classified as an ultra-luminous X-ray source (ULX; Kaaret et
al., 2017). We have already noted in Section 2.1 that for sub-Eddington
accretion, the puffy disk luminosity is proportional to the mass accretion
rate, that is, the efficiency does not change appreciably with $\dot{m}$ for
$\dot{m}\leq 1.5$. Assuming approximate proportionality also between
luminosity and black hole mass (with the caveat that, in fact, radiative
simulations do not strictly scale with the mass of the black hole), we could
give an approximate formula for beamed observed isotropic luminosity of a
puffy disk

$L_{\rm iso}\sim 5\times
10^{39}\left(\frac{M}{10M_{\odot}}\right)\left(\frac{\dot{M}}{\dot{M}_{\rm
Edd}}\right)\left(\frac{b}{0.25}\right)^{-1}\frac{\rm erg}{\rm s}.$ (7)

This may allow an interpretation of those ULXs that comprise a BHXB in the
puffy disk accretion framework.

Figure 5: Properties of puffy disk spectra for a $M=10M_{\odot}$ black hole.
Left: isotropic radiative power per logarithmic frequency interval, as a
function of the observer’s inclination, from $10^{\circ}$ (nearly face-on) to
$80^{\circ}$ (nearly edge on). Right: isotropic radiative power per
logarithmic frequency interval for three GRRMHD simulations, corresponding to
different mass accretion rates $\dot{m}$.

##  4 Spectral properties of puffy disks

Traditionally, BHXBs and ULXs are studied with phenomenological models, i.e.,
by composing a total model from as many components as needed to emulate the
observed spectrum: a thin disk, a corona, an iron line, a reflection bump,
etc. Given the comparatively simple and analytic nature of the model
components, this has the crucial advantage that it allows to identify and
characterize different disk properties and effects very efficiently. However,
different underlying disk structures necessarily modify the contribution of
other components, and this can lead to a misreading of the importance of added
effects. Within the GRRMHD framework, we can self-consistently solve for the
structure of the accretion system, taking into account the radiation feedback,
influence of the magnetic field, vertical stratification, and outflows. All
these effects influence the spectra in a highly non-trivial and non-linear way
that cannot be fully reproduced by adding and/or multiplying a series of
analytic model components. Synthetic spectra from GRRMHD simulations of hard
state BHXBs were recently discussed by Dexter et al. (2021).

We present synthetic puffy disk spectra computed using the HEROIC code in
Figs. 5-6. In the left panel we explore the influence of the inclination,
expanding on the total luminosity discussion from Section 3. The inclination
$\theta$ has a strong impact on the puffy disk spectra because of the beaming
at low inclinations, and obscuration at high inclinations. Between
$\theta=10^{\circ}$ and $\theta=80^{\circ}$ the maximum value of the spectral
energy distribution (SED) changes by two orders of magnitude. However, the
location of the spectral peak is only mildly affected, and remains around 3
keV. At low inclinations the SED develops a hump at around 20 keV, related to
the Comptonized radiation produced in the innermost part of the accretion
flow.

The right panel of Fig. 5 explores the differences between our three GRRMHD
simulations for a range of mass accretion rates and two specific inclinations.
Apart from total radiative output scaling, they indicate similar features. A
small shift of the peak towards higher energies is visible in the
$\dot{m}=0.6$ model. This is related to the lower optical depth of the system.

Figure 6: Comparison of puffy disk spectra to multi-color blackbody spectra of
thin (kerrbb) and slim (slimbh) disk models of the same mass $M=10M_{\odot}$,
mass accretion rate $\dot{m}=0.6$, and observer’s inclination $\theta$.

Finally, in Fig. 6 we contrast the puffy disk spectra with those of a thin NT
disk and a slim disk in their xspec implementations (kerrbb and slimbh,
respectively). These are two relativistic, height-integrated $\alpha$-disk
models based on Novikov & Thorne (1973) and Sądowski (2009), respectively.
They both assume a geometrically thin and optically thick disk that consists
of a sum of annular black body emitters. This allows one to determine the
total disk luminosity by radially integrating $L(r)\propto T^{4}_{\rm
eff}(r)$. Both models parametrize the turbulent viscosity by assuming stress
to be proportional to pressure via a constant $\alpha$ (Shakura & Sunyaev,
1973). The slim disk differs from the thin NT disk on account of its
additional advective cooling which transports an increasing fraction of the
disk thermal energy towards the black hole as the mass accretion rate rises.

Slim and thin disk spectra appear similar at the moderate accretion rate of
$0.6M_{\rm Edd}$. Neither of them accounts for Compton scattering which can be
modelled separately (see Section 5). This causes a difference with respect to
the puffy disk, where Comptonization is self-consistently accounted for.

Figure 7: Comparison of NT kerrbb (red lines) and puffy disk (gray lines)
spectra at an intermediate inclination of $\theta=48^{\circ}$ for a mass
accretion rate $\dot{m}=0.6$. The absolute value of the difference between the
models is shown by dashed lines. In the range of 10-40 keV the puffy disk
spectrum can be approximated with a power law corresponding to a photon index
$\Gamma=4$ (blue line).  Figure 8: Luminosity-hardness diagram for the
synthetic puffy disk spectra. Vertical axis represents the 2-10 keV band
isotropic luminosity in the Eddington units. Hardness is defined as the ratio
of photon counts in the 4-10 keV and 2-4 keV bands. For comparison, the
location of the kerrbb and slimbh spectra for $\dot{m}=0.6$ are shown as well;
for small and moderate inclinations they both indicate less luminous and
softer spectra than a puffy disk at the same mass accretion rate.

For intermediate inclinations, around $\theta=48^{\circ}$, where neither
beaming nor obscuration are dominant, the puffy disk peak of SED is similar to
that of the NT disk model. The comparison between NT and puffy disk spectra at
$\theta=48^{\circ}$ is shown in Fig. 7. There is a deficit of power in the
puffy spectrum with respect to the NT disk in the range of 0.1-1.0 keV and an
excess of puffy disk power above 3 keV. The puffy disk spectrum above 10 keV
is described reasonably well by a power law with a $\nu L_{\nu}$ slope equal
to $-2$, corresponding to a photon index $\Gamma=4$.

The effect of inclination on the luminosity and hardness of the spectra for
different disk models is shown in Fig. 8. As discussed in Section 3, the
observer’s inclination has a huge impact on the observed puffy disk
luminosity. The hardness of puffy disk spectra diminishes between
$30^{\circ}$-60∘ inclination because of the obscuration of the inner (beamed
hard X-ray) region of the disk. In comparison, both the thin and the slim disk
are very soft at low inclinations as they completely lack the hard Compton
component. They become harder under increasing viewing angles due to
relativistic beaming in directions close to the plane of the disk, primarily
affecting the inner region of the disk (e.g., Kulkarni et al., 2011).
Moreover, at $\dot{m}=0.6$ slim disk spectra tend be harder than those of an
NT thin disk, owing to advection that traps thermal energy in the disk,
transports it inwards and releases it at smaller radii as harder emission. The
non-monotonic relation between the spectral hardness and mass accretion rate
for the puffy disks, seen in Fig. 8, is more difficult to understand and can
be attributed to an interplay between the temperature of the disk, which
depends on the mass accretion rate, and the obscuration.

##  5 Modelling puffy disk spectra

Figure 9:  Synthetic puffy disk spectra fitted by standard xspec models. The
synthetic puffy disk spectra (black crosses) are calculated for a hypothetical
BHXB at $D=10$ kpc, with a black hole mass $M=10M_{\odot}$, a mass accretion
rate $\dot{m}=0.6$, and inclinations $\theta=10^{\circ}$ (left panel) and
$\theta=60^{\circ}$ (right panel). They are fitted with a disk model plus a
Compton component: slimbh+nthcomps (solid lines) and kerrbb+nthcompk (dashed
lines).

During their outbursts, X-ray binaries pass a series of different disk states
and typically spend a comparatively long time in a highly luminous, disk
dominated state with a characteristic black body dominated SED. Strictly
thermal dominated spectra, i.e., spectra with a disk-to-total flux fraction of
above 75% (in the energy range 2-20 keV), a weak non-thermal power continuum,
and very weak variability (rms < 1%), can be used to extract the black hole
spin (continuum fitting method, e.g., Remillard & McClintock, 2006; Done et
al., 2007; McClintock et al., 2014).

If we were to observe spectra from an astrophysical object consistent with our
puffy disk model, how would they be interpreted by the standard spectral
fitting tools? To address this question we generated synthetic test spectra
from a puffy disk simulation in the energy range 0.1-100 keV, and fitted them
with the xspec software package (Arnaud, 1996). The synthetic puffy spectra
correspond to a GRRMHD simulation of a Schwarzschild black hole of spin
$a_{*}=0$, mass accretion rate $\dot{m}=0.6$ (as defined in Eq. 1), and mass
$M=10M_{\odot}$, at a distance D = 10 kpc, seen from two different
inclinations, $\theta=10^{\circ}$ and $\theta=60^{\circ}$. We fit these
spectra with a fully relativistic (including gravitational bending of photon
trajectories) thin disk model, kerrbb (Li et al., 2005), and a fully
relativistic, advective slim disk model, slimbh (Sądowski, 2009; Straub et
al., 2011). Each disk component is accompanied by a thermal Comptonization
component nthcomp (Zdziarski et al., 1996; Życki et al., 1999) which accounts
for seed photons originating in the disk that are upscattered to higher
energies by a homogeneous corona.

For each disk model the dynamical BHXB parameters, mass $M$, distance $D$, and
inclination $\theta$, are assumed to be known and fixed, while the black hole
spin $a_{*}$ and total (spherically-integrated, inclination-independent) disk
luminosity $L_{\rm disk}$ (in kerrbb, equivalently, the mass accretion rate)
are fitted for. The viscosity parameter is fixed at $\alpha=0.1$ in each
model. In the Compton component, the temperature of the seed photons, $T_{\rm
bb}$, and the electrons in the corona, $T_{\rm e}$, as well as the photon
index, $\Gamma$, are fitted for. Usually, in order to reduce the number of
free parameters, $T_{\rm bb}$ is fixed at some “reasonable” value. Here, we
leave it free to assess how well the $T_{\rm bb}$ parameter matches the puffy
disk temperature, as presented in Fig. 2.

| puffy $\theta=10^{\circ}$ | puffy $\theta=60^{\circ}$
---|---|---
| slimbh | kerrbb | slimbh | kerrbb
$a_{*}$ | $0.82^{+0.08}_{-0.16}$ | $0.86^{+0.08}_{-0.14}$ | $0.55^{+0.11}_{-0.29}$ | $0.48^{+0.12}_{-0.41}$
$L_{\rm disk}/L_{\rm Edd}$ | $0.73^{+0.15}_{-0.16}$ | $0.79^{+0.1}_{-0.1}$ | $0.21^{+0.19}_{-0.24}$ | $0.17^{+0.05}_{-0.03}$
$\Gamma$ | $2.85^{+0.17}_{-0.20}$ | $2.72^{+0.24}_{-0.24}$ | $3.62^{+0.27}_{-0.29}$ | $3.67^{+0.34}_{-0.27}$
$T_{\rm e}$ (keV) | $14.95^{+5.31}_{-2.68}$ | $13.05^{+4.62}_{-2.43}$ | $19.51^{+17.34}_{-5.75}$ | $19.26^{+21.81}_{-5.58}$
$T_{\rm bb}$ (keV) | $0.50^{+0.08}_{-0.07}$ | $0.67^{+0.22}_{-0.15}$ | $1.01^{+0.57}_{-0.31}$ | $0.97^{+0.48}_{-0.29}$
$\chi^{2}/dof$ | 9.32/29 | 8.51/29 | 11.58/29 | 27.97/29
Table 1: xspec fits to synthetic puffy disk spectra with $a_{*}=0$ and
$\dot{m}=0.6$. Each of the disk models, slimbh and kerrbb, is accompanied by
the thermal Comptonization model nthcomp.

Given the phenomenological nature of the models, they are both able to fit the
puffy disk spectra with high fidelity (low $\chi^{2}$ per degree of freedom,
$\chi^{2}/dof$), as shown in Fig. 9 and Tab. 1. We assumed a 10% error budget
on the simulated spectral measurements. While the fits are well behaved, the
models formally overfit our simulated data due to the assumed uncertainties.
Nevertheless, the fits demonstrate that all quantities can be well
constrained. Furthermore, typical values for BHXBs are obtained, albeit with
results that are quite different from the true values of $a_{*}=0$ and
$\dot{m}=0.6$.

The highly beamed low inclination ($\theta=10^{\circ}$) puffy spectrum is
preferentially fitted with a high black hole spin $a_{*}\approx 0.7-0.8$ for
both disk models, and a (thermal) disk luminosity that is slightly higher than
the true luminosity of the puffy disk, but much lower than its isotropic
luminosity at this inclination angle. The Comptonization component fits a
photon index of about $\Gamma$ = 2.8, a disk temperature of $T_{\rm bb}\approx
0.6\,$keV, and an electron temperature of approximately $T_{\rm e}$ = 14 keV.
While the thin and slim disk models require a cranked up spin parameter to fit
the low-energy part of puffy spectra, the Compton model recovers the shape of
the high-energy emission with typical values for a warm corona (see, e.g.,
Różańska et al., 2015). The dim, inclined puffy disk ($\theta=60^{\circ}$) is
fitted with moderate black hole spins and very low disk luminosities, combined
with a warm corona of $\Gamma\approx 3.6$, $T_{\rm bb}\approx$ 1 keV, and
$T_{\rm e}\approx 19$ keV.

The measurement of spin is impaired by the large Compton contributions, which
is a known obstacle to the continuum fitting method and one of the reasons why
only the strictly thermal dominated states can be effectively used to measure
black hole spins, as pointed out by, e.g., Gierliński & Done (2004); Done &
Kubota (2006); McClintock et al. (2006); Middleton et al. (2006). Only limited
attempts have been made to extend the continuum fitting method towards
modelling spectra with an appreciable power law component, e.g., Steiner et
al. (2009).

In Tab. 1 we show the parameters of the best fits with uncertainties
corresponding to 90% confidence intervals. The classical disk models gravely
overestimate the black hole spin, with respect to $a_{*}=0$ that was assumed
in the underlying GRRMHD simulation. This overestimate reflects the high
temperatures in the puffy disk, seen in Fig. 2. The mass accretion of the
fitted models correspond to $\dot{m}\approx 0.4$ for low inclination and
$\dot{m}\approx 0.1$ for large inclination (with $\dot{m}$ defined using the
zero-spin efficiency convention from Eq. 1), i.e., the fitted $\dot{M}$ value
is only equal to 2/3 and 1/6, respectively, of the true value in the
simulation. The mismatch with respect to the true value of $\dot{m}=0.6$ is
explained by the increased disk efficiency at high black hole spin, while the
variation with the inclination angle of the fitted mass accretion value by
beaming/obscuration effects as well.

Puffy disk spectra do not describe the strictly thermal dominated state in
BHXBs; they are not strongly dominated by thermal radiation but comprise of a
considerable amount of Compton up-scattered photons (see Fig. 5). Their global
(0.1-100 keV) blackbody-to-total flux ratio, as measured by xspec’s thin and
slim disk models, lies for $\theta=10^{\circ}$ at 54% and for
$\theta=60^{\circ}$ at 83%, whereas in the 2-20 keV range the disk
contribution is just 44% and 72%, respectively. In combination with their
rather steep high-energy slopes (see $\Gamma$ in Tab. 1) and high
normalizations these properties can be related to the spectra observed in two
accretion disk states. On the one hand they resemble the so-called
intermediate states in BHXBs. These are complex states associated to
transitions between the hard and the soft state and the occurrence of quasi-
periodic variability (Méndez & van der Klis, 1997; Remillard & McClintock,
2006). Although puffy disks can exhibit, depending on the viewing angle, a
considerably larger (even super-Eddington) apparent luminosity than average
BHXBs, a spectral resemblance is expected due to the presence of a Comptonized
corona. On the other hand, puffy disk spectra show morphological similarities
to certain ULXs, namely those that show a soft ultraluminous state. There is
increasing evidence that at least some ULXs are powered by accreting stellar
mass black holes that produce three specific types of ultraluminous disk
states: a single-peaked broadened ultra-soft disk, a two-component hard
ultraluminous, and a soft ultraluminous state. The latter is typically
described as a cool disk with a warm corona component (see Gladstone et al.,
2009; Sutton et al., 2013; Mondal et al., 2021), which resembles our puffy
disk spectra.

##  6 Discussion

In this paper we discussed the observational signatures, and in particular the
spectra, of the puffy disk—a numerical model of stable accretion in the range
of luminosities (currently) above $0.5L_{\rm Edd}$. This is the regime where
standard geometrically thin disk models fail to correctly predict the source
spectrum, which happens at luminosities above $0.3L_{\rm Edd}$. The puffy disk
solution captures the physics of such luminous accretion flows. In our
approach, radiation is incorporated into the GRMHD simulations, allowing for a
self-consistent development of a vertically stratified, magnetic and radiation
pressure dominated, accretion structure. A low density “puffy” region,
Keplerian, rapidly accreting, and geometrically and optically thick, develops
above the denser disk core zone; spectroscopically, it resembles a “warm
corona”.

The observable luminosity of puffy disks depends much more strongly on the
observer’s inclination to the system axis than is the case for standard
geometrically thin disk models. At high inclinations the emission is obscured
by the geometrically thick corona blocking the view of the central region of
the disk. The development of an obscuring geometrically thick structure during
the luminous thermal state may explain unusual variability of BHXB sources,
e.g., V404 Cygni. At low inclinations the emission is geometrically collimated
by the puffy region, increasing the observed radiative flux. Such a strong
beaming effect is absent in thin disks, and may be relevant for the
interpretation of at least some ULX sources.

We investigated synthetic spectra of puffy disks and attempted to fit them
with standard fitting models available in xspec. We found that a comptonized
NT disk fits the puffy disk spectra well in terms of the statistical fit
quality metrics. However, for our model of a non-spinning black hole the
currently available fitting tools erroneously find a large positive spin
value. This offers a theoretical interpretation for known problems with
fitting BHXB spectra in luminous states, although we can not currently provide
a quantitative characterization of these effects in the Kerr metric. For
instance, for LMC X-3 Straub et al. (2011) and Steiner et al. (2014) reported
that the spin estimate of $a_{*}\approx 0.21$ decreases with luminosity due to
the incompleteness of physical models (lack of Comptonization) at luminosities
larger than about 0.3 $L_{\rm Edd}$. For the Schwarzschild puffy disk solution
we can only see an increase of the spin estimate with respect to the ground
truth value of $a_{*}=0$, and indeed this is what happens as the standard
xspec models try to emulate the puffy disk spectra by tuning the free
parameters.

The puffy spectra resemble intermediate states between soft and hard emission
states of BHXBs. The complex intermediate states are associated with quasi-
periodic variability and there currently exists no established accretion disk
model to describe these states. We have demonstrated that numerically
calculated puffy disk spectra, which self-consistently model accretion disks
with a significant Comptonized atmosphere (warm corona), could qualitatively
describe these intermediate states. It remains to be seen if GRRMHD framework
could provide a quantitatively correct interpretation of the real
astrophysical sources, extending the region of applicability of the continuum
fitting spin estimation method.

It is currently computationally too expensive to produce a fine-mesh grid in
the parameter space of puffy disks, to calculate a table model that could be
implemented for public use in the xspec software package. For this reason, we
cannot at present explicitly fit puffy disk models to observational X-ray
spectra. However, this is clearly one of the directions that should be pursued
in the future. Currently, we do not have puffy disk solutions for the Kerr
metric with $a_{*}\neq 0$, nor for $L<0.6L_{\rm Edd}$. The latter solutions
are urgently needed, as there as yet exist no known stable solutions for disk
accretion at rates corresponding to the luminosity of many BHXB sources.

## Acknowledgements

The authors thank Andrew Chael, Brandon Curd, and Aleksander Sądowski for
support with the KORAL code, as well as Michal Bursa, Silke Britzen, and
Monika Mościbrodzka for valuable comments and discussions. The computations in
this work were supported by the PLGrid Infrastructure through which access to
the Prometheus supercomputer, located at ACK Cyfronet AGH in Kraków, was
provided. This work was supported in part by the Polish NCN grant
2019/33/B/ST9/01564 and the Black Hole Initiative at Harvard University, which
is funded by grants the John Templeton Foundation and the Gordon and Betty
Moore Foundation to Harvard University. M.A. acknowledges the Polish NCN grant
2015/19/B/ST9/01099, while A.R. – the NCN grant 2021/41/B/ST9/04110. M.A. and
G.T. the Czech Science Foundation grant No. GX21-06825X, and D.L. the internal
grant of SU SGS/13/2019 and the Student Grant Foundation of the Silesian
University in Opava, Grant No. $\mathrm{SGF/1/2020}$, which has been carried
out within the EU OPSRE project entitled “Improving the quality of the
internal grant scheme of the Silesian University in Opava”, reg. number:
$\mathrm{CZ.02.2.69/0.0/0.0/18\\_054/0014696}$. R.N. acknowledges the support
of the NSF grants OISE-1743747 and AST-1816420. We also thank Alexandra
Elbakyan for her contributions to the open science initiative.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author.

## References

  * Abramowicz & Fragile (2013) Abramowicz M. A., Fragile P. C., 2013, Living Reviews in Relativity, 16, 1
  * Abramowicz et al. (1978) Abramowicz M., Jaroszynski M., Sikora M., 1978, A&A, 63, 221
  * Abramowicz et al. (1988) Abramowicz M. A., Czerny B., Lasota J. P., Szuszkiewicz E., 1988, ApJ, 332, 646
  * Arnaud (1996) Arnaud K. A., 1996, in Jacoby G. H., Barnes J., eds, Astronomical Society of the Pacific Conference Series Vol. 101, Astronomical Data Analysis Software and Systems V. p. 17
  * Bagińska et al. (2021) Bagińska P., Różańska A., Czerny B., Janiuk A., 2021, ApJ, 912, 110
  * Begelman & Pringle (2007) Begelman M. C., Pringle J. E., 2007, MNRAS, 375, 1070
  * Ciesielski et al. (2012) Ciesielski A., Wielgus M., Kluźniak W., Sądowski A., Abramowicz M., Lasota J. P., Rebusco P., 2012, A&A, 538, A148
  * Davis & Hubeny (2006) Davis S. W., Hubeny I., 2006, ApJS, 164, 530
  * Dexter et al. (2021) Dexter J., Scepi N., Begelman M. C., 2021, ApJ, 919, L20
  * Done & Kubota (2006) Done C., Kubota A., 2006, MNRAS, 371, 1216
  * Done et al. (2007) Done C., Gierliński M., Kubota A., 2007, A&ARv, 15, 1
  * EHT Collaboration et al. (2019) EHT Collaboration et al., 2019, ApJ, 875, L1
  * Giacconi et al. (1971) Giacconi R., Kellogg E., Gorenstein P., Gursky H., Tananbaum H., 1971, ApJ, 165, L27
  * Gierliński & Done (2004) Gierliński M., Done C., 2004, MNRAS, 347, 885
  * Gladstone et al. (2009) Gladstone J. C., Roberts T. P., Done C., 2009, MNRAS, 397, 1836
  * Gronkiewicz & Różańska (2020) Gronkiewicz D., Różańska A., 2020, A&A, 633, A35
  * Homan et al. (2001) Homan J., Wijnands R., van der Klis M., Belloni T., van Paradijs J., Klein-Wolt M., Fender R., Méndez M., 2001, ApJS, 132, 377
  * Kaaret et al. (2017) Kaaret P., Feng H., Roberts T. P., 2017, ARA&A, 55, 303
  * King (2009) King A. R., 2009, MNRAS, 393, L41
  * King et al. (2001) King A. R., Davies M. B., Ward M. J., Fabbiano G., Elvis M., 2001, ApJ, 552, L109
  * Kluzniak & Kita (2000) Kluzniak W., Kita D., 2000, arXiv e-prints, pp astro–ph/0006266
  * Kubota & Makishima (2004) Kubota A., Makishima K., 2004, ApJ, 601, 428
  * Kulkarni et al. (2011) Kulkarni A. K., et al., 2011, MNRAS, 414, 1183
  * Lančová et al. (2019) Lančová D., et al., 2019, ApJ, 884, L37
  * Lasota (2001) Lasota J.-P., 2001, New Astron. Rev., 45, 449
  * Lasota et al. (2016) Lasota J. P., Vieira R. S. S., Sądowski A., Narayan R., Abramowicz M. A., 2016, A&A, 587, A13
  * Levine et al. (1996) Levine A. M., Bradt H., Cui W., Jernigan J. G., Morgan E. H., Remillard R., Shirey R. E., Smith D. A., 1996, ApJ, 469, L33
  * Li et al. (2005) Li L.-X., Zimmerman E. R., Narayan R., McClintock J. E., 2005, ApJS, 157, 335
  * Lightman & Eardley (1974) Lightman A. P., Eardley D. M., 1974, ApJ, 187, L1
  * McClintock et al. (2006) McClintock J. E., Shafee R., Narayan R., Remillard R. A., Davis S. W., Li L.-X., 2006, ApJ, 652, 518
  * McClintock et al. (2014) McClintock J. E., Narayan R., Steiner J. F., 2014, Space Sci. Rev., 183, 295
  * Méndez & van der Klis (1997) Méndez M., van der Klis M., 1997, ApJ, 479, 926
  * Middleton et al. (2006) Middleton M., Done C., Gierliński M., Davis S. W., 2006, MNRAS, 373, 1004
  * Mishra et al. (2020) Mishra B., Begelman M. C., Armitage P. J., Simon J. B., 2020, MNRAS, 492, 1855
  * Mondal et al. (2021) Mondal S., Różańska A., Bagińska P., Markowitz A., De Marco B., 2021, A&A, 651, A54
  * Motta et al. (2017) Motta S. E., Kajava J. J. E., Sánchez-Fernández C., Giustini M., Kuulkers E., 2017, MNRAS, 468, 981
  * Narayan & McClintock (2005) Narayan R., McClintock J. E., 2005, ApJ, 623, 1017
  * Narayan et al. (2016) Narayan R., Zhu Y., Psaltis D., Sądowski A., 2016, MNRAS, 457, 608
  * Narayan et al. (2017) Narayan R., Saądowski A., Soria R., 2017, MNRAS, 469, 2997
  * Novikov & Thorne (1973) Novikov I. D., Thorne K. S., 1973, in Dewitt C., Dewitt B. S., eds, Black Holes (Les Astres Occlus). pp 343–450
  * Oda et al. (2009) Oda H., Machida M., Nakamura K. E., Matsumoto R., 2009, ApJ, 697, 16
  * Ohsuga & Mineshige (2011) Ohsuga K., Mineshige S., 2011, ApJ, 736, 2
  * Ohsuga et al. (2009) Ohsuga K., Mineshige S., Mori M., Kato Y., 2009, PASJ, 61, L7
  * Page & Thorne (1974) Page D. N., Thorne K. S., 1974, ApJ, 191, 499
  * Regev & Gitelman (2002) Regev O., Gitelman L., 2002, A&A, 396, 623
  * Remillard & McClintock (2006) Remillard R. A., McClintock J. E., 2006, ARA&A, 44, 49
  * Różańska et al. (1999) Różańska A., Czerny B., Życki P. T., Pojmański G., 1999, MNRAS, 305, 481
  * Różańska et al. (2011) Różańska A., Madej J., Konorski P., Sądowski A., 2011, A&A, 527, A47
  * Różańska et al. (2015) Różańska A., Malzac J., Belmont R., Czerny B., Petrucci P. O., 2015, A&A, 580, A77
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 500, 33
  * Shakura & Sunyaev (1976) Shakura N. I., Sunyaev R. A., 1976, MNRAS, 175, 613
  * Sikora (1981) Sikora M., 1981, MNRAS, 196, 257
  * Sądowski (2016) Sądowski A., 2016, MNRAS, 459, 4397
  * Sądowski et al. (2013) Sądowski A., Narayan R., Tchekhovskoy A., Zhu Y., 2013, MNRAS, 429, 3533
  * Sądowski et al. (2014) Sądowski A., Narayan R., McKinney J. C., Tchekhovskoy A., 2014, MNRAS, 439, 503
  * Sądowski (2009) Sądowski A., 2009, ApJS, 183, 171
  * Sądowski & Narayan (2015) Sądowski A., Narayan R., 2015, MNRAS, 454, 2372
  * Sądowski et al. (2011) Sądowski A., Abramowicz M., Bursa M., Kluźniak W., Lasota J. P., Różańska A., 2011, A&A, 527, A17
  * Sądowski et al. (2017) Sądowski A., Wielgus M., Narayan R., Abarca D., McKinney J. C., Chael A., 2017, MNRAS, 466, 705
  * Steiner et al. (2009) Steiner J. F., McClintock J. E., Remillard R. A., Narayan R., Gou L., 2009, ApJ, 701, L83
  * Steiner et al. (2014) Steiner J. F., McClintock J. E., Orosz J. A., Remillard R. A., Bailyn C. D., Kolehmainen M., Straub O., 2014, ApJ, 793, L29
  * Straub et al. (2011) Straub O., et al., 2011, A&A, 533, A67
  * Sutton et al. (2013) Sutton A. D., Roberts T. P., Middleton M. J., 2013, MNRAS, 435, 1758
  * Vincent et al. (2011) Vincent F. H., Paumard T., Gourgoulhon E., Perrin G., 2011, Classical and Quantum Gravity, 28, 225011
  * Wielgus et al. (2016) Wielgus M., Yan W., Lasota J. P., Abramowicz M. A., 2016, A&A, 587, A38
  * Zdziarski et al. (1996) Zdziarski A. A., Johnson W. N., Magdziarz P., 1996, MNRAS, 283, 193
  * Zhang et al. (2000) Zhang S. N., Cui W., Chen W., Yao Y., Zhang X., Sun X., Wu X.-B., Xu H., 2000, Science, 287, 1239
  * Zhu & Narayan (2013) Zhu Y., Narayan R., 2013, MNRAS, 434, 2262
  * Zhu et al. (2015) Zhu Y., Narayan R., Sądowski A., Psaltis D., 2015, MNRAS, 451, 1661
  * Życki et al. (1999) Życki P. T., Done C., Smith D. A., 1999, MNRAS, 309, 561