††thanks: Authors contributed equally††thanks: Authors contributed equally

# Supplemental Material for “Giant Spin-Orbit Torque in Cr-based Janus
Transition Metal Dichalcogenides”

Libor Vojáček  Université Grenoble Alpes, CEA, CNRS, IRIG-Spintec, 38000
Grenoble, France  Joaquín Medina Dueñas  ICN2 — Catalan Institute of
Nanoscience and Nanotechnology, CSIC and BIST, Campus UAB, Bellaterra, 08193
Barcelona, Spain  Department of Physics, Universitat Autònoma de Barcelona
(UAB), Campus UAB, Bellaterra, 08193 Barcelona, Spain  Jing Li  Université
Grenoble Alpes, CEA, Leti, F-38054, Grenoble, France  Fatima Ibrahim
Université Grenoble Alpes, CEA, CNRS, IRIG-Spintec, 38000 Grenoble, France
Aurélien Manchon  Aix-Marseille Université, CNRS, CINAM, Marseille 13288,
France  Stephan Roche  ICN2 — Catalan Institute of Nanoscience and
Nanotechnology, CSIC and BIST, Campus UAB, Bellaterra, 08193 Barcelona, Spain
ICREA — Institució Catalana de Recerca i Estudis Avançats, 08010 Barcelona,
Spain  Mairbek Chshiev  Université Grenoble Alpes, CEA, CNRS, IRIG-Spintec,
38000 Grenoble, France  Institut Universitaire de France, 75231 Paris, France
José H. García  ICN2 — Catalan Institute of Nanoscience and Nanotechnology,
CSIC and BIST, Campus UAB, Bellaterra, 08193 Barcelona, Spain

##  S1 Ab initio calculation details

The ab initio ground state calculations were performed using density
functional theory (DFT) as implemented in the Vienna ab initio simulation
package (VASP) kresse_ab_1993; kresse_efficiency_1996 with the GGA-PBE
exchange-correlation functional perdew_generalized_1996-1. The Hubbard $U$
correction in Dudarev’s formulation dudarev_electron-energy-loss_1998 was
applied to the Cr d orbitals with an effective Hubbard correction
$U_{\mathrm{eff}}=3.0$ eV, a value that makes the band structure closely match
the experimentally-measured one zhang_room-temperature_2021. A Cr
pseudopotential with semicore p electrons was chosen and an energy cutoff of
330 eV was used for the plane-wave basis. Brillouin zone was sampled with a
19x19x1 $\Gamma$-centered mesh. The van der Waals interaction was approximated
by the DFT-D3 method grimme_consistent_2010 with the Becke-Johnson damping
grimme_effect_2011. Forces were minimized below $10^{-3}$ eV/$\mathrm{\AA}$
during the full structural relaxation and total energy below $10^{-7}$ eV to
achieve the electronic self-consistency. Spin-orbit coupling was included
except in the relaxation step.

The Wannier tight-binding models were constructed using the Wannier90 package,
with a 22-orbital basis (10 Cr-d and $2\times 6$ chalcogen-p orbitals) and a
large $25\times 25$ real-space supercell that limits the real-space
interactions. A smaller supercell might very well be sufficient, see Sec. S4.
The Wannier energy window was chosen to tightly encompass the 22 bands with
the Cr-d and chalcogen-p character. The frozen energy window has an upper
limit to avoid including an intersecting parasitic s/dz2 band and keeps an
additional margin of 0.2 eV from this band, a value chosen to optimize the
RMSE quality of the CrTe2 Wannier model (see Sec. S4). Disentanglement and
maximal localization were performed with 200 and 100 maximum steps,
respectively. The spin operator in the maximally-localized Wannier function
(MLWF) basis was derived, as described below.

### The spin operator in real space

Accessing the spin expectation values of the interpolated band structure with
Wannier tight-binding models requires deriving the Wannier real-space spin
operator.

Just like any other operator that can be obtained from the ab initio
calculation (for all the ab initio k-points $\bm{q}$), the spin operator in
the basis of ab initio eigenstates $\mathcal{S}_{mn}^{\mathrm{H}}(\textbf{q})$
can be converted into its real-space representation by first applying the
semi-unitary matrix $V_{m^{\prime}m}(\textbf{q})$ to convert $\mathcal{S}$
from the Hamiltonian gauge to the Wannier gauge

$\mathcal{S}_{mn}^{\mathrm{W}}(\bm{q})=\sum_{m^{\prime}n^{\prime}}V_{m^{\prime}m}^{\dagger}(\bm{q})\cdot\mathcal{S}_{m^{\prime}n^{\prime}}^{\mathrm{H}}(\textbf{q})\cdot
V_{n^{\prime}n}(\bm{q}),$ (S1)

where the spin operator $\hat{\mathcal{S}}$ projected onto the basis of the ab
initio eigenstates $\psi$ in its matrix form
$\mathcal{S}_{m^{\prime}n^{\prime}}^{\mathrm{H}}(\textbf{q})=\left\langle\psi_{m^{\prime}\bm{q}}|\hat{\mathcal{S}}|\psi_{n^{\prime}\bm{q}}\right\rangle$
is obtained with the help of the vaspspn module of WannierBerri
tsirkin_high_2021. Primed indices run over the (larger) space of
disentanglement bands.

Follows the direct Fourier sum over the ab initio grid

$\mathcal{S}_{mn}(\bm{R})\equiv\frac{1}{N_{\bm{q}}}\sum_{\bm{q}}\mathcal{S}_{mn}^{\mathrm{W}}(\bm{q})\cdot
e^{-i\bm{q}\cdot\bm{R}}\,,$ (S2)

with $N_{\bm{q}}$ the number of ab initio grid points $\bm{q}$.

Having $\mathcal{S}_{mn}(\bm{R})$ at hand, its interpolation
$\overline{\mathcal{S}}_{mn}^{\mathrm{H}}(\bm{k})$ to an arbitrary k-vector
$\bm{k}$ involves an inverse Fourier sum
$\overline{\mathcal{S}}_{mn}^{\mathrm{W}}(\bm{k})=\sum_{\textbf{R}}\mathcal{S}_{mn}(\bm{R})\cdot
e^{i\bm{k}\cdot\bm{R}}$ over the real-space lattice vectors $\bm{R}$ followed
by a rotation back to the Hamiltonian gauge of the original eigenstates
$\overline{\mathcal{S}}_{mn}^{\mathrm{H}}(\bm{k})=(U^{\dagger}\cdot\mathcal{S}^{\mathrm{W}}\cdot
U)_{mn}$, where $U_{mn}$ is a unitary matrix which diagonalizes the
interpolated Hamiltonian
$\overline{\mathcal{H}}_{mn}^{\mathrm{W}}(\bm{k})=\sum_{\textbf{R}}\mathcal{H}_{mn}(\bm{R})\cdot
e^{i\bm{k}\cdot\bm{R}}$.

##  S2 Electrostatics of CrXTe

The field-like torque is directly linked to the Rashba effect attributed to an
inversion symmetry breaking, which comes from an electrical dipolar field
across the monolayer plane.

### Internal E-fields in CrXTe monolayers

The spontaneous internal E-fields in the Janus CrXTe come from the charge
imbalance, as shown in Fig. S1. This creates a work function difference and an
internal electric field $\sim$ 1–2 V/nm in CrXTe. We calculate these values
from the work function difference $\Delta\phi$ [Fig. S2(a)] divided by the
thickness of the layer, which is derived from the positions of the chalcogen
atom centers extended by their Wigner-Seitz radii.

Figure S1: Electrostatics of CrXTe monolayers. (a,b) The (Bader) charge
transfer $Q$ is different at the X and Te sides, giving rise to (c) a large
internal electric field $\sim 2$ V/nm. The equivalent out-of-plane E-field
over CrTe2 would have to be an enormous
$\epsilon_{\mathrm{CrTe2}}\,2~{}\mathrm{V/nm}\sim 100$ V/nm, demonstrating the
potency of CrXTe.

### Relative permittivity of CrTe2 monolayer

Subtracting the $xy$-averaged potential of the CrTe2 monolayer under an out-
of-plane E-field of 2 V/nm and under no E-field, we can estimate its relative
permittivity from the gradient of the potential inside the layer and in
vacuum, see Fig. S2. In vacuum, the potential gradient corresponds to the
imposed $E_{\mathrm{ext}}=2$ V/nm, which is largely screened inside CrTe2 to
$E_{\mathrm{ind}}=0.038$ V/nm. The $E_{\mathrm{ind}}/E_{\mathrm{ext}}$ ratio
gives the relative permittivity of CrTe2, $\epsilon_{\mathrm{CrTe2}}=52.9.$
This is a prefactor that should be applied to the internal electric fields of
CrXTe

$\mathcal{E}_{\mathrm{ext}}^{\mathrm{equivalent\,for\,CrTe2}}=\epsilon_{\mathrm{CrTe2}}\,\mathcal{E}_{\mathrm{int}}^{\mathrm{CrXTe}}$
(S3)

to get an equivalent external E-field over CrTe2 to achieve a comparable
Rashba splitting.

Figure S2: Electrostatics of CrXTe. (a) The ($xy$-averaged) potential in
CrSeTe is different at the Se and Te sides. This work function difference
translates into an internal electric dipole. (b) The relative permittivity of
CrTe2 monolayer can be calculated from the difference between the
$xy$-averaged potential of CrTe2 monolayer under an external
$E_{\mathrm{ext}}=2$ V/nm and without an external E-field. The potential
gradient corresponds to the E-field magnitude. The ratio of the E-field
outside (in vacuum) and inside the layer corresponds to
$\epsilon_{\mathrm{CrTe2}}$.

### The equivalent external E-field to the Janus internal E-field

In Tab. 1, the ratio between $\chi_{\mathrm{FL}}$ of CrXY and gated CrTe2
gives an estimation of the required $E_{\mathrm{ext}}$ which should be applied
to CrTe2 to obtain a similar $\chi_{\mathrm{FL}}$ in CrTe2 as are
spontaneously present in CrXY. The values are colossal 30 to 60 V/nm. They
correspond well to a simple multiplication of the CrXY internal E-fields
$E_{\mathrm{int}}$ with the relative permittivity of CrTe2
$\epsilon_{\mathrm{CrTe2}}=52.9$, giving a 113 and 66 V/nm.

Table 1: The calculated internal E-fields $E_{\mathrm{int}}$ in CrXY
monolayers and an estimated equivalent external field
$E_{\mathrm{ext}}^{\mathrm{CrTe2}}\simeq\epsilon_{\mathrm{CrTe2}}\cdot
E_{\mathrm{int}}$ that should act on CrTe2 to reach a similar effect. Further,
the calculated field-like torques $\chi_{\mathrm{FL}}$ in Janus CrXY
monolayers and an estimated external E-field that should act on CrTe2 to reach
similar torque magnitudes. Both these estimations result in huge
$E_{\mathrm{ext}}^{\mathrm{CrTe2}}\approx 30$ to 50 V/nm, illustrating the
superiority of Janus CrXY for producing the Rashba effect and field-like
torques.

| $E_{\mathrm{int}}$ | $E_{\mathrm{ext}}^{\mathrm{CrTe2}}\simeq\epsilon_{\mathrm{CrTe2}}\cdot E_{\mathrm{int}}$ | $\chi_{\mathrm{FL}}$ | $E_{\mathrm{ext}}^{\mathrm{CrTe2}}\simeq 2.0\cdot\chi_{\mathrm{FL}}/\chi_{\mathrm{FL}}^{\mathrm{CrTe2\,@\,2.0\,V/nm}}$
---|---|---|---|---
| V/nm | V/nm | 10${}^{3}\hbar$/2e (Ohm m)-1 | V/nm
CrSTe | 2.14 | 113 | 163 | 30
CrSeTe | 1.25 | 66 | 308 | 57

  * •

$\epsilon_{\mathrm{CrTe2}}=52.9;\;\;\;\chi_{\mathrm{FL}}^{\mathrm{CrTe2\,@\,2.0\,V/nm}}=10.9$
10${}^{3}\hbar$/2e (Ohm m)-1

##  S3 Strain (lattice parameter variation)

Although theoretical studies suggest that the in-plane lattice parameter of
CrTe2 is reduced in its monolayer form compared to bulk liu_structural_2022,
experimental measurements zhang_room-temperature_2021; meng_anomalous_2021
suggest that ultrathin films of CrTe2 retain their lattice constant reasonably
close to bulk. The experimentally measured lattice parameter of ultrathin
CrTe2 films $a=3.77\,\mathrm{\AA}$ meng_anomalous_2021 is $\approx 4\%$ larger
than our calculated monolayer value, as shown in Fig. S3(a). Since Janus CrXTe
monolayers have not yet been fabricated and their experimental structure may
also slightly differ from our ab initio predictions, we deem it necessary to
provide our results as a function of strain, which is not necessarily a
physical strain, but rather a variation of the lattice parameter with respect
to the predicted ab initio value.

As shown in Fig. S3(b-c), the consequence of such positive strain is an
increase in both the effective magnetic anisotropy and the Heisenberg exchange
coupling, linked to the Curie temperature.

Figure S3: The magnetic anisotropy and exchange coupling with lattice
parameter variation (”strain”). 0% corresponds to the ab initio-predicted
lattice constant. (a) The experimental lattice parameter of CrTe2
meng_anomalous_2021 is $\approx 4\%$ larger than the ab initio prediction.
Both (b) the effective uniaxial anisotropy and (c) the nearest-neighbor
exchange interaction (and the consequent Ising Curie temperature) largely
increase for higher lattice parameters in CrXTe. The vertical dashed line
corresponds to the $\approx 4\%$ latt. param. difference that might be
expected from the comparison of DFT and experiment of CrTe2 from (a). Inset of
(b) shows the 2x1x1 supercell used for the Heisenberg exchange parameter
calculation.

The appropriate effective magnetic anisotropy in case of 2D magnets
lado_origin_2017

$K_{\mathrm{eff}}=K_{\mathrm{u}}+Z(J_{z}-J_{x})+E_{\mathrm{demag}}$ (S4)

is a sum of the calculated single-ion anisotropy $K_{\mathrm{u}}$, the
(super)exchange interaction anisotropy $Z(J_{z}-J_{x})$ with the number of
nearest neighbors $Z=6$ and the demagnetizing energy
$E_{\mathrm{demag}}=-\mu_{0}M_{\mathrm{s}}^{2}/2$ with its small value of
$\approx-0.07$ and $-0.06$ mJ/m2 for CrSTe and CrSeTe respectively.

Note that the dominant contribution is expected to be the anisotropic
superexchange $Z(J_{z}-J_{x})$, while the single-ion contribution
$K_{\mathrm{u}}$ should be negligible lado_origin_2017.

All these quantities can be calculated from the ab initio total energies,
considering the XXZ Hamiltonian lado_origin_2017

$H=-K_{u}\sum_{i}\left(S_{i}^{z}\right)^{2}-\frac{J_{x}}{2}\underset{i\neq
j}{\sum\sum}\left(S_{i}^{x}S_{j}^{x}+S_{i}^{y}S_{j}^{y}\right)-\frac{J_{z}}{2}\underset{i\neq
j}{\sum\sum}\left(S_{i}^{z}S_{j}^{z}\right)$ (S5)

where $i$ runs over all Cr atoms and $j$ runs over all of their 6 nearest
neighbors (thereby double counting is present and redeemed by the factor
$\frac{1}{2}$ coming with $J$). Constructing a 2x1x1 supercell [see the inset
of Fig. S3(b)] of the CrXTe unit cells with both ferromagnetic (FM) and
antiferromagnetic (AF) configurations along both the $x$ and $z$ directions,
it follows from Eq. S5 that the ab initio ground state energies are
$E^{\mathrm{FM}}_{x}=E_{0}-6S^{2}J_{x}$,
$E^{\mathrm{AF}}_{x}=E_{0}+2S^{2}J_{x}$,
$E^{\mathrm{FM}}_{z}=E_{0}-2S^{2}K_{\mathrm{u}}-6S^{2}J_{z}$ and
$E^{\mathrm{AF}}_{z}=E_{0}-2S^{2}K_{\mathrm{u}}+2S^{2}J_{z}$, from which

$\displaystyle J_{x(z)}$
$\displaystyle=\left(-E^{\mathrm{FM}}_{x(z)}+E^{\mathrm{AF}}_{x(z)}\right)/8S^{2}\,,$
(S6) $\displaystyle K_{\mathrm{u}}$
$\displaystyle=\left(E^{\mathrm{FM}}_{x}+3E^{\mathrm{AF}}_{x}-E^{\mathrm{FM}}_{z}-3E^{\mathrm{AF}}_{z}\right)/8S^{2}\,,$

where $S=3/2$ is the magnitude of the Cr spin lado_origin_2017.

In Fig. S3(c) we also plot the strain-enhanced Curie temperature in the Ising
limit torelli_calculating_2018

$T_{\mathrm{C}}^{\mathrm{Ising}}=\frac{S^{2}\tilde{T_{\mathrm{c}}}}{k_{\mathrm{B}}}J_{x}\,,$
(S7)

where the constant dimensionless $\tilde{T_{\mathrm{c}}}=3.64$ for a hexagonal
lattice torelli_calculating_2018, $k_{\mathrm{B}}$ is the Boltzmann constant.

##  S4 Quality of the Wannier models

We provide a quantitative estimate of the quality of the derived Wannier
tight-binding models, which aim to accurately describe the ab initio band
energies and spin expectation values. For this, the band structure along the
$\mathrm{K-}\Gamma\mathrm{-M-K}$ with 101 points along each segment is
calculated both from ab initio and the derived Wannier TB model and compared
to obtain the RMSE errors as

$\mathrm{RMSE}(E)=\sqrt{\frac{1}{N}\sum_{\mathrm{\langle
kpath\rangle}}(E^{\mathrm{DFT}}-E^{\mathrm{WannTB}})^{2}}$ (S8)

and analogically for the spin expectation values. $N$ is the number of summed
data points. The averaging is performed inside a $E_{\mathrm{F}}\pm 0.5$ eV
energy window.

### Quality convergence with the Wannier supercell size

The real-space interactions need to be limited by a large enough supercell to
converge the RMSE below a chosen threshold, as plotted in Fig. S4. A $25\times
25$ supercell provides sufficient accuracy. The saturation of the RMSE of
energy for CrXTe indicates a (small) constant offset of the whole energy
spectrum due to slightly ($\approx 1$ meV) misaligned Fermi energy values.

Figure S4: The root mean square error of energy and spin expectation values of
the Wannier TB model vs. the original DFT calculation, averaged along the
standard $\mathrm{K-}\Gamma\mathrm{-M-K}$ k-path and within $\pm$ 0.5 eV
window around the Fermi level. The RMSE of energy exponentially decreases with
increasing Wannier supercell (interaction cutoff). The vertical line denotes
the $25\times 25$ supercell, used for further transport calculations.

### Quality of the final models

Figure S5 shows the RMSE of energy and spin for CrTe2 (with and without
external $E$-field), CrSTe, and CrSeTe calculated with different ab initio
magnetization directions
$\vec{m}=(\sin\theta\cos\phi,\,\sin\theta\sin\phi,\,\sin\theta)$. The errors
are fairly small: $\approx 1$ meV for energy and $\approx 0.01$ for spin
expectation values.

Figure S5: The root mean square error of energy and spin expectation values of
the Wannier TB model vs. the original DFT calculation, averaged along the
standard $\mathrm{K-}\Gamma\mathrm{-M-K}$ k-path and within $\pm$ 0.5 eV
window around the Fermi level. The azimuthal and polar angles $\theta$ and
$\phi$ denote calculations with different magnetization directions. The errors
are small: around 1 meV in energy and 0.01 in spin expectation values.

##  S5 Non-Equilibrium Spin Density and Symmetry-Allowed Contributions

Consistent with the linear transport regime, we expand the non-equilibrium
spin density $\bm{S}$ to linear order in terms of its $\mathcal{E}_{j}$
($j=x,y$) dependence. On the other hand, since the magnetization direction
does not constitute a perturbative parameter, we expand $\bm{S}$ up to second
order with respect to $\hat{\bm{m}}$. We express the non-equilibrium spin
density as $\bm{S}=\bm{S}_{0}+\bm{S}_{\text{3m}}$, where the first term groups
the torque contributions which are allowed in arbitrary non-centrosymmetric
systems, while the second represents torque contributions enabled by the
reduced 3m symmetry. The first contribution reads

$\begin{split}\bm{S}_{0}=&\,\chi_{\text{FL}}[1-\xi_{\text{FL}}(m_{x}^{2}+m_{y}^{2})]\hat{\bm{z}}\times\bm{\mathcal{E}}-\chi_{\text{DL}}\hat{\bm{m}}\times(\hat{\bm{z}}\times\bm{\mathcal{E}})-\chi_{\text{DL}}^{z}(\hat{\bm{m}}\cdot\bm{\mathcal{E}})\hat{\bm{z}}-\chi_{\parallel}[\hat{\bm{m}}\cdot(\hat{\bm{z}}\times\bm{\mathcal{E}})]\hat{\bm{m}}-\chi_{\parallel}^{z}[\hat{\bm{m}}\cdot(\hat{\bm{z}}\times\bm{\mathcal{E}})]m_{z}\hat{\bm{z}}\text{
,}\end{split}$ (S9)

showing that all systems allow for many contributions beyond the conventional
field-like and damping-like torques jmedina2024, respectively represented in
the terms proportional to $\chi_{\text{FL}}$ and $\chi_{\text{DL}}$, where the
field-like torque is modulated by the total in-plane magnetization via
$\xi_{\text{FL}}$. $\chi_{\text{DL}}^{z}$ accounts for an anisotropic damping
of the in-plane and out-of-plane magnetization components towards the steady
state along $\hat{\bm{z}}\times\bm{\mathcal{E}}$. $\chi_{\parallel}$
represents a non-equilibrium spin density component parallel to the
magnetization which cannot exert torque, while $\chi_{\parallel}^{z}$
corresponds to an out-of-plane anisotropy of the aforementioned component
which, along with the field-like torque, induces an elliptic precession of the
magnetization about $\hat{\bm{z}}\times\bm{\mathcal{E}}$. A set of
unconventional torques are additionally enabled by the system’s 3m symmetry,
which are represented in

$\begin{split}\bm{S}_{\text{3m}}=&\,\chi_{\text{3m}}[(\mathcal{E}_{x}m_{y}+\mathcal{E}_{y}m_{x})\hat{\bm{x}}+(\mathcal{E}_{x}m_{x}-\mathcal{E}_{y}m_{y})\hat{\bm{y}}]-\chi_{\text{3m}}^{(2)}m_{z}[(-\mathcal{E}_{x}m_{x}+\mathcal{E}_{y}m_{y})\hat{\bm{x}}+(\mathcal{E}_{x}m_{y}+\mathcal{E}_{y}m_{x})\hat{\bm{y}}]\\\
&+\chi_{\text{3m}}^{z}(\mathcal{E}_{x}m_{x}^{2}-\mathcal{E}_{x}m_{y}^{2}-2\mathcal{E}_{y}m_{x}m_{y})\hat{\bm{z}}\text{
.}\end{split}$ (S10)

The term proportional to $\chi_{\text{3m}}$ represents the so-called 3m
torque, while the other second-order terms stem from a contribution of the
form $\hat{\bm{m}}\times\bm{s}_{\text{3m}}$, with $\bm{s}_{\text{3m}}$ the
aforementioned first order 3m spin density. The 3m torque derives from a
current-induced in-plane magnetic anisotropy which, when applying a driving
electric field $\bm{\mathcal{E}}=\mathcal{E}_{x}\hat{\bm{x}}$ perpendicular to
the mirror plane, lies along the $m_{x}=-m_{y}$ axis johansen_current_2019.
Thus, the first order 3m torque $\chi_{\text{3m}}$ accounts for precession
about the anisotropy axis, while the second order torques
$\chi_{\text{3m}}^{(2)}$ and $\chi_{\text{3m}}^{z}$ account for a damping
parallel to it.

##  S6 Transport simulations: complete results

### Strainless

We present the complete spin-orbit torque results for all the studied systems,
at 0% strain. In order to disentangle the symmetry-allowed torque
contributions we compute the non-equilibrium spin density for a set of 18
magnetization directions for each system, consisting of 8 evenly spaced
$\varphi$ points for $\theta=\pi/2$, 4 evenly spaced $\varphi$ points for
$\theta=\pi/4,3\pi/4$, and the two poles $\theta=0,\pi$, with
($\theta,\varphi$) the magnetization polar and azimuthal angles. The spin-
torque conductivities are shown in Fig. S6. We observe that all torque
contributions are much larger in the Janus structures than in the electric
field assisted CrTe2, while they are vanishing in centrosymmetric CrTe2.
Additionally, the results for the longitudinal conductivity show that it
remains of the same order in Janus vs non-Janus systems.

Figure S6: Longitudinal conductivity $\sigma$, in units of
$10^{6}(\Omega\,\text{m})^{-1}$, and all spin-torque conductivities
$\tau_{\alpha}$, in units of $10^{3}\frac{\hbar}{2e}(\Omega\,\text{m})^{-1}$,
for all systems. The longitudinal conductivity of CrTe2 is unaffected by
$\mathcal{E}_{z}$

### Strain-dependent

Fig. S7 shows the field-like, damping-like and 3m torques as a function of
strain. All torques are enhanced with strain, retrieving a perpendicular
magnetic anisotropy for strain $\geq 1\%$.

Figure S7: Field-like, damping-like and 3m spin-torque conductivities, in
units of $10^{3}\frac{\hbar}{2e}(\Omega\,\text{m})^{-1}$ as a function of
strain. Solid and dashed curves correspond to perpendicular and in-plane
magnetic anisotropy respectively. The corresponding lattice parameter $a_{0}$
is indicated in the colorbar.


# Smooth Actual EIRP Control for EMF Compliance
with Minimum Traffic Guarantees

Lorenzo Maggi, Aloïs Herzog, Azra Zejnilagic, Christophe Grangeat The authors
are with Nokia. Emails<EMAIL_ADDRESS>
###### Abstract

To mitigate Electromagnetic Fields (EMF) human exposure from base stations,
international standards bodies define EMF emission requirements that can be
translated into limits on the “actual” Equivalent Isotropic Radiated Power
(EIRP), i.e., averaged over a sliding time window. We aim to enable base
stations to adhere to these constraints while mitigating any impact on user
performance. Specifically, our objectives are to: i) ensure compliance with
actual EIRP control methods defined in IEC 62232:2022, ii) guarantee a minimum
EIRP level, and iii) prevent resource shortages at all times. Initially, we
present exact and conservative algorithms to compute the set of feasible EIRP
consumption values that fulfill constraints i) and ii) at any given time,
referred to as the actual EIRP “budget”. Subsequently, we design a control
method based on Drift-Plus-Penalty theory that preemptively curbs EIRP
consumption only when needed to avoid future resource shortages.

###### Index Terms:

EMF exposure, actual maximum approach, EIRP control, GBR, Drift-Plus-Penalty.

##  I Introduction

The Radio Frequency (RF) Electromagnetic Field (EMF) limits for human exposure
are established by international bodies in [1], [2]. Considering the
variability of EMF emissions from base stations (BTS), the “actual maximum
approach” described in [3] provides methods for assessing EMF compliance using
time averaged, also called _actual_ , EMF emission values. These include the
Equivalent Isotropic Radiated Power (EIRP), that can be controlled during BTS
operations. According to channel modeling studies, setting a threshold on the
actual EIRP at four times below the maximum value has a minimal impact on end-
user performance [4]. However, lower thresholds may be necessary when
employing higher directivity antennas [5]. Among the existing approaches for
actual EMF emission control, [6] proposes a proportional-differentiating power
control procedure that limits the usable bandwidth, while [7] optimizes
precoding for MIMO systems by assuming full channel knowledge.

In this paper we propose an actual EIRP control method that complies with [3].
Our approach aims to minimize the impact on user performance by guaranteeing a
minimum EIRP level and preemptively reducing emissions to prevent resource
shortages. While presented for EIRP control, our method can be adapted to
limit _any_ function of the emitted power.

##  II Problem formulation

The EIRP emitted by an antenna in the azimuth/elevation direction
$(\phi,\theta)$ is the product of the radiated power $P$ and the antenna gain
$G(\phi,\theta)$:

$\mathrm{EIRP}(\phi,\theta):=P\,G(\phi,\theta).$ (1)

In free-space, the power density $S$ measured at distance $R$ from the antenna
is proportional to the EIRP:

$S(\phi,\theta)=\frac{\mathrm{EIRP}(\phi,\theta)}{4\pi R^{2}}\ \
[\mathrm{W}/\mathrm{m}^{2}].$ (2)

We consider a discrete time grid indexed by $t$ and we denote the EIRP
_consumption_ $c_{t}$ as the maximum EIRP over a predefined set of
azimuth/elevation angles $\mathcal{A}$, called _segment_ :

$c_{t}:=\max_{(\phi,\theta)\in\mathcal{A}}\langle\mathrm{EIRP}(\phi,\theta)\rangle_{t}$
(3)

where $\langle.\rangle_{t}$ averages the EIRP over period $t$. Following [3],
the _actual_ EIRP consumption, averaged over a sliding window of $W$ periods,
shall not exceed the configured threshold $\overline{C}$:

$\frac{1}{W}\sum_{i=t-\min(t,W-1)}^{t}\\!\\!\\!c_{i}(\gamma_{i})\leq\overline{C},\qquad\forall\,t\geq
0$ (4)

where consumption $c_{t}$ is capped by the control variable $\gamma_{t}$:

$0\leq c_{t}(\gamma_{t})\leq\gamma_{t},\qquad\forall\,t\geq 0.$ (5)

To enforce (5), various methods can be used, such as reducing the antenna
transmit power [8], the transmission bandwidth [6], or the beamforming gain
[7]. Our solution, controlling $\gamma$, is agnostic to the specific power
limitation method chosen.

To guarantee a _minimum service level_ , we aim to keep the control $\gamma$
above the minimum EIRP level $\rho\overline{C}$:

$\gamma_{t}\geq\rho\overline{C},\qquad\forall\,t\geq 0.$ (6)

where $\rho\in[0,1]$. Fulfilling (6) is crucial especially in the presence of
Guaranteed Bit-Rate (GBR) traffic.

As it turns out, various control policies can jointly fulfill the constraints
(4),(6). The simplest one is a greedy policy that at each period $t$ sets
$\gamma_{t}$ to the actual EIRP _budget_ $\Gamma_{t}$, being the highest
possible consumption such that (4),(6) can be satisfied.

Such greedy policy is optimal in case of low traffic, where
$c_{t}\ll\gamma_{t}$. In this case, the only effect of preemptively curbing
consumption to a value $\gamma_{t}<\Gamma_{t}$ is an increase of the
transmission delay. Instead, under the high traffic regime, setting
$\gamma_{t}=\Gamma_{t}$ leads to budget depletion, after which the control
drops to its minimum value $\rho\overline{C}$. One natural method to avoid
such drops is to increase $\rho$. Yet, this leads to overly conservative
policies: for instance, we observe that when $\rho$ approaches 1, $\Gamma_{t}$
tends to the actual EIRP threshold $\overline{C}$.

Thus, we study a method which decides whether and when preemptively curbing
consumption to a value $\gamma_{t}<\Gamma_{t}$ to avoid any future budget
depletion. We formalize the desired preemptive control behavior as the one
maximizing the average $\alpha$-fairness [9] of the control $\gamma$ over
time, where the $\alpha$-fair function is defined as:

$f(x):=\left\\{\begin{array}[]{ll}\frac{x^{1-\alpha}}{1-\alpha}\quad\mathrm{if\
}\alpha\geq 0,\,\alpha\neq 1\\\ \log(x)\quad\mathrm{if\
}\alpha=1\end{array}\right..$ (7)

Hence, our problem becomes:

$\displaystyle\max_{\gamma}$
$\displaystyle\,\lim_{T\rightarrow\infty}\frac{1}{T}\sum_{t=0}^{T-1}\mathbb{E}\left[f(\gamma_{t})\right]$
(8) $\displaystyle\mathrm{s.t.}$
$\displaystyle\,\eqref{eq:cap},\eqref{eq:EMFconstr},\eqref{eq:rhoC}$

where the expectation $\mathbb{E}$ is with respect to the user demand.

Note that, if $\alpha=1$, $\alpha$-fairness coincides with proportional
fairness; as $\alpha$ grows, it tends to max-min fairness. Since $f(.)$ is
concave, small values of $\gamma$ are severely penalized, thus encouraging
$\gamma$ to evolve smoothly over time. At the same time, since $f(.)$ is
increasing, higher values of the control are preferred. Hence, the policy
solving the problem (8) shall preemptively curb the EIRP consumption only
“when needed”.

##  III Actual EIRP budget

We start tackling our EIRP control problem (8) by characterizing its
feasibility set, i.e., the set of control values $\gamma$ for which
constraints (4)-(6) are fulfilled.

######  Proposition 1.

Define the variable $\Omega_{t}$ as:

$\Omega_{t}:=\max_{0\leq
k\leq\min(t,W-1)}\sum_{i=1}^{k}\left(c_{t-i}-\rho\overline{C}\right).$ (9)

To fulfill constraints (4)-(6), $\gamma_{t}$ must satisfy:

$\rho\overline{C}\leq\gamma_{t}\leq\rho\overline{C}+\overline{C}(1-\rho)W-\Omega_{t}:=\Gamma_{t},\quad\forall\,t\geq
0.$ (10)

We refer to the upper bound $\Gamma_{t}$ of the feasibility set as _budget_.
All proofs can be found in the Appendix.

###  III-A Budget computation

We now study how to compute the term $\Omega_{t}$, defining the budget
$\Gamma_{t}$ via (10). A naive approach would compute all $W$ terms in (9) and
take their maximum. Its complexity is _quadratic_ in the window length $W$,
which may be prohibitive in practice, where typically $W\approx 10^{3}$ and
the period duration is a few hundreds of ms. Thus, we next design more
efficient methods. First, we define the auxiliary variable:

$\Omega_{t}^{n}:=\,\max_{k=0,\dots,n}\sum_{i=1}^{k}(c_{t-i}-\rho\overline{C}),\quad
0\leq n\leq t,\ \forall\,t\geq 0$ (11)

and we call $\ell_{t}^{n}$ the argument of the maximum in (11). Note that
$\Omega_{t}:=\Omega_{t}^{\min(t,W-1)}$ and $\ell_{t}:=\ell_{t}^{\min(t,W-1)}$,
where $\ell_{t}$ is the argument of maximum in (9). We next prove a useful
recursive property for $\Omega_{t}^{n}$, where $[.]^{+}:=\max(.,0)$.

######  Proposition 2.

The following recursive relation holds:

$\Omega_{t+1}^{n+1}=\,\left[\Omega_{t}^{n}+c_{t}-\rho\overline{C}\right]^{+},\quad
n\leq t,\ \forall\,t\geq 0.$ (12)

Proposition 2 suggests that, in order to compute $\Omega_{t}$, one can apply
(12) recursively over the last $W-1$ EIRP consumption samples, as described in
Algorithm 1. The computational complexity of such procedure is linear in $W$.

Algorithm 1 _From-scratch_ computation of _exact_ budget $\Gamma_{t}$

Fix period $t$. Initialize $\Omega_{t}:=0$.
for $i=\min(t,W-1),\dots,1$:

  * Update $\Omega_{t}:=\left[\Omega_{t}+c_{t-i}-\rho\overline{C}\right]^{+}$

return $\Gamma_{t}=\rho\overline{C}+\overline{C}(1-\rho)W-\Omega_{t}$

###  III-B Iterative computation

Algorithm 1 computes $\Omega_{t}$ from scratch at each period $t$, i.e., it
does not reuse computations performed at previous periods. Yet, under some
conditions, $\Omega_{t+1}$ can be computed from $\Omega_{t}$ with _constant_
(and negligible, in practice) complexity. This hinges on the following
properties of variables $\Omega,\ell$.

######  Property 1.

If $\ell_{t}<W-1$, then

$\Omega_{t+1}=\left[\Omega_{t}+c_{t}-\rho\overline{C}\right]^{+}.$ (13)

######  Property 2.

If $\Omega_{t+1}^{n+1}>0$, then $\ell_{t+1}^{n+1}=\ell_{t}^{n}+1$. Else,
$\ell_{t+1}^{n+1}=0$.

######  Property 3.

If $c_{t-i}\geq\rho\overline{C}$ for $0\leq i\leq W-1$, then

$\Omega_{t+1}=\Omega_{t}+c_{t}-c_{t-W+1}.$ (14)

Properties 1-3 suggest that, if the last $W$ consumption samples exceed
$\rho\overline{C}$, then $\Omega_{t+1}$ can be updated via (14). Else, if
$\ell_{t}<W-1$, then $\Omega_{t+1}$ can be derived as (13). In both cases, the
incremental complexity is _constant_ with respect to $W$. Else, if
$\ell_{t}=W-1$, then $\Omega_{t+1}$ is computed from scratch via the original
Algorithm 1. Finally, the value of $\ell_{t}$ is updated via Property 2. This
procedure is formalized in Algorithm 2.

Algorithm 2 _Update_ of _exact_ budget $\Gamma_{t}$

Initialize $\ell_{0}:=0,\,c_{i}:=0,\Omega_{0}:=0\ \forall\,i<0$.
for period $t\geq 0$:

  * if $c_{t-i}\geq\rho\overline{C},\ \forall\,i=0,\dots,W-1$ then:

    * Set $\Omega_{t+1}=\Omega_{t}+c_{t}-c_{t-W+1}$

  * else if $\ell_{t}<W-1$ then: $\Omega_{t+1}=\left[\Omega_{t}+c_{t}-\rho\overline{C}\right]^{+}$

  * else: Compute $\Omega_{t+1}$ via Algorithm 1

  * if $\Omega_{t+1}>0$ then: $\ell_{t+1}=\ell_{t}+1$

  * else: $\ell_{t+1}=0$

  * return $\Gamma_{t}=\rho\overline{C}+\overline{C}(1-\rho)W-\Omega_{t}$

###  III-C Low-complexity conservative approximation

The _worst-case_ complexity of iterative Algorithm 2 is still _linear_ in $W$,
although its _average_ complexity is lower than that of Algorithm 1. The
natural question is whether one can do better, especially if computational
resources are a bottleneck.

We here propose an alternative iterative method that provides a _conservative_
estimate of the EIRP budget with _constant_ complexity, i.e., not depending on
the window length $W$.

First, it is convenient to introduce the variable $\widetilde{\Omega}_{t}$,
which is a variant of the original $\Omega_{t}$ defined in (9), and is
computed as the sum of the consumption excess with respect to
$\rho\overline{C}$ over the sliding window:

$\widetilde{\Omega}_{t}:=\sum_{i=1}^{\min(t,W-1)}\left[c_{t-1}-\rho\overline{C}\right]^{+}.$
(15)

The associated feasibility set, analogous to (10), writes:

$\rho\overline{C}\leq\gamma_{t}\leq\rho\overline{C}+\overline{C}(1-\rho)W-\widetilde{\Omega}_{t}:=\widetilde{\Gamma}_{t},\quad\forall\,t\geq
0.$ (16)

Then, $\widetilde{\Gamma}_{t}$ represents the maximum feasible control under
the assumption that all past consumption values below the minimum guaranteed
level $\rho\overline{C}$ are _overestimated_ as $\rho\overline{C}$.
Consequently, it is expected that the approximated budget is generally _lower_
than the exact budget, as we formalize next.

######  Proposition 3.

$\widetilde{\Gamma}_{t}\leq\Gamma_{t}$, for all periods $t\geq 0$.

Thus, this approximate method produces _overly conservative_ policies,
reducing EIRP consumption more than necessary.

On the positive side, the budget can be updated at each period via a simple
rule with _constant_ complexity, only requiring to add the latest consumption
excess and subtracting the one exiting the new sliding window, as formalized
in Algorithm 3.

Algorithm 3 _Update_ of _conservative_ budget $\widetilde{\Gamma}_{t}$

Initialize $\widetilde{\Omega}_{0}:=0,\ c_{i}:=0,\ \forall\,i<0$.
for period $t\geq 1$:

  * Compute $\widetilde{\Omega}_{t}=\widetilde{\Omega}_{t-1}+\left[c_{t-1}-\rho\overline{C}\right]^{+}-\left[c_{t-W}-\rho\overline{C}\right]^{+}$

  * return $\widetilde{\Gamma}_{t}=\rho\overline{C}+\overline{C}(1-\rho)W-\widetilde{\Omega}_{t}$

Figure 1: Impact of different control strategies on EIRP consumption under
identical input traffic. Parameters:
$W=10,\rho=0.15,\alpha=1,\beta=0.95,V=15$.

##  IV Smooth EIRP control

After characterizing the feasible EIRP control set, we now turn to the
solution of the original problem (8).

Two natural strategies arise. One approach is a greedy strategy, where
$\gamma_{t}$ is simply set to $\Gamma_{t}$ (or
$\gamma_{t}=\widetilde{\Gamma}_{t}$, if the conservative budget is used) for
all $t$’s. Here, EIRP is constrained only when necessary. However, under high
traffic conditions, the budget can deplete, leading to $\gamma$ being pushed
down to its minimum value $\rho\overline{C}$ until the budget refills.
Alternatively, to prevent resource shortages, one can cautiously set
$\gamma_{t}=\overline{C}$ for all $t$’s. In this case, the control is
consistently smooth, but low. While optimal for high-load conditions, this
cautious strategy proves overly conservative when traffic is low and EIRP
control is unnecessary.

Ideally, the optimal strategy would not restrict EIRP during low traffic,
while proactively reducing EIRP as traffic increases or when the budget is
limited to prevent future resource shortages. Figure 1 elaborates further on
such intuitions.

To effectively control $\gamma$, we propose a method based on Drift-Plus-
Penalty (DPP) theory.

###  IV-A Preliminaries on Drift-Plus-Penalty (DPP)

Drift-Plus-Penalty (DPP) is a stochastic optimization technique typically
employed to maximize a network utility while guaranteeing queue stability
[10]. Formally speaking, DPP solves problems of the following form:

$\displaystyle\min_{x_{t}\in A(\omega_{t})}$
$\displaystyle\,\lim_{T\rightarrow\infty}\frac{1}{T}\sum_{t=0}^{T}\mathbb{E}\left[P(x_{t},\omega_{t})\right]:=p^{*}$
(17) $\displaystyle\mathrm{s.t.}$
$\displaystyle\,\lim_{T\rightarrow\infty}\frac{1}{T}\sum_{t=0}^{T-1}\mathbb{E}\left[Y(x_{t},\omega_{t})\right]\leq
0$ (18)

where $x_{t}$ is the control variable, $\\{\omega_{t}\\}_{t}$ are i.i.d.
random variables, $P,Y$ are deterministic functions and $A(\omega_{t})$ is the
feasible control set. By defining the _virtual queue_ $Q$ evolving as
$Q_{t+1}=\left[Q_{t}+Y(x_{t},\omega_{t})\right]^{+}$ for all $t\geq 0$, the
constraint (18) can be equivalently formulated as $Q$ being mean-rate stable,
i.e., $\lim_{t\rightarrow\infty}\mathbb{E}[Q_{t}]/t=0$. The DPP control is
computed as:

$x_{t}^{*}=\operatorname{argmin}_{x\in
A(\omega_{t})}Q_{t}Y(x,\omega_{t})+VP(x,\omega_{t})$ (19)

which jointly reduces the queue Lyapunov drift (defined as the increment
$Q_{t+1}^{2}-Q_{t}^{2}$) and the objective function. The parameter $V>0$
regulates the trade-off between the achieved performance, away from the
optimal $p^{*}$ by $\mathcal{O}(V)$, and the average size $\mathcal{O}(1/V)$
of the virtual queue, as shown in [10].

###  IV-B DPP-based heuristic

We can draw a _partial_ parallel between our problem (8) and the DPP
formulation (17-18) by setting the control $x_{t}:=\gamma_{t}$, the objective
function $P(x_{t},\omega_{t}):=-f(\gamma_{t})$ and the function
$Y(x_{t},\omega_{t}):=c_{t}(\gamma_{t})-\overline{C}$. In this case, (19)
becomes:

$\min_{\gamma\in[\rho\overline{C},\Gamma_{t}]}\overline{Q}_{t}c_{t}(\gamma)-Vf^{(\alpha)}(\gamma)$
(20)

where the virtual queue $\overline{Q}$ evolves as follows:

$\overline{Q}_{t+1}=\left[\overline{Q}_{t}+c_{t}(\gamma_{t})-\overline{C}\right]^{+},\quad\forall\,t\geq
0$ (21)

with $\overline{Q}_{0}=0$. Yet, (20) cannot be computed in practice, since
$c_{t}(\gamma)$ is an _unknown_ function of $\gamma$. Indeed, the variable
$\omega_{t}$ (which can be interpreted as the consumption $c_{t}$ if _no_
power control is applied, i.e., $\gamma_{t}=\infty$) is unknown. We will then
bound (20) from above by replacing $c_{t}(\gamma)$ with $\gamma$ via (5) and
minimize $\overline{Q}_{t}\gamma-Vf^{(\alpha)}(\gamma)$ instead, leading to
the expression:

$\gamma_{t}^{*}=\min\left(\max\left(V/\overline{Q}_{t}^{1/\alpha},\,\rho\overline{C}\right),\,\Gamma_{t}\right),\quad\forall\,t\geq
0.$ (22)

Note that, as traffic increases and EIRP consumption consistently surpasses
the threshold $\overline{C}$, the virtual queue grows, leading to a decrease
in $\gamma$. By preemptively curbing consumption before the budget depletes,
we prevent resource scarcity in the future, as illustrated in Figure 1.

A second discrepancy between formulations (8) and (17)-(18) lies in the
averaging window, which spans $W$ samples in (8) and has an infinite duration
in (17)-(18). It turns out that this leads the virtual queue $\overline{Q}$ to
empty at most every $W$ samples in our finite window scenario, as stated
below.

######  Proposition 4.

For any list $\tau=[t^{\prime},\dots,t+W-1]$ of $W$ successive periods, there
exists a period $t\in\tau$ where $\overline{Q}_{t}=0$.

As the virtual queue empties, the control $\gamma^{*}_{t}$ as defined in (22)
peaks at the budget $\Gamma_{t}$ at most every $W$ samples. This can be
particularly problematic in the case of high traffic, where any value
$\gamma_{t}>\overline{C}$ results in future budget depletion. To avoid this
behavior, we redefine the virtual queue as follows:

$\overline{Q}^{\beta}_{t+1}=\left[\overline{Q}^{\beta}_{t}+c_{t}(\gamma_{t})-\beta\overline{C}\right]^{+},\quad
0\leq\beta<1,\ \forall\,t\geq 0$ (23)

where $\beta$ artificially inflates the virtual queue, thus preventing the
undesirable behavior demonstrated in Proposition 4. We describe our DPP-based
EIRP control method in Algorithm 4.

Algorithm 4 DPP-based power control for EMF compliance

Initialize $\overline{Q}^{\beta}_{0}:=0$.
for period $t\geq 0$:

  * Compute the budget $\Gamma_{t}$ via Algorithm 2

  * return ${\gamma_{t}^{\beta}}^{*}=\min(\max(V/(\overline{Q}^{\beta}_{t})^{1/\alpha},\,\rho\overline{C}),\,\Gamma_{t})$

  * Assess consumption $c_{t}$ and compute $\overline{Q}^{\beta}_{t+1}$ as in (23)

Figure 2: Exact budget $\Gamma$ versus conservative $\widetilde{\Gamma}$ when
traffic demands are generated via a sparse zipf distribution that equals 0
with probability $\ell$.

###  IV-C Numerical evaluations

For our numerical evaluations, we generated a series of traffic demands of
$d_{t}$ bits. Demands are sparse: with probability $\ell$, variable $d_{t}$ is
drawn from a Zipf distribution, else $d_{t}=0$. The control $\gamma_{t}$ caps
the requested EIRP at period $t$. The unserved demand is buffered and shows up
again in the next period.

Figure 2 illustrates the relationship between the exact and the conservative
budget, $\Gamma$ and $\widetilde{\Gamma}$, respectively. In the case of sparse
($\ell\approx 0$) and high-load ($\ell>0.3$) traffic, the two solutions tend
to coincide, as theory suggests. Indeed, as consumption $c$ vanishes, $\Omega$
and $\widetilde{\Omega}$ both tend to zero; conversely, as consumption $c$
increases and is consistently higher than the guaranteed $\rho\overline{C}$,
the approximate excess $(c-\rho\overline{C})^{+}$ equals the actual excess
$c-\rho\overline{C}$, resulting in an exact budget, i.e.,
$\widetilde{\Gamma}=\Gamma$. However, in the middle load regime, the exact
solution $\Gamma$ outperforms the conservative one $\widetilde{\Gamma}$. In
fact, when consumption fluctuates around $\rho\overline{C}$, only the exact
solution is able to account for the past “unused” (positive) consumption
deficit $\rho\overline{C}-c$ to produce the current budget, while
$\widetilde{\Gamma}$ approximates the deficit as null, leading to an overly
conservative budget.

Figure 3 illustrates the relationship between the traffic sparsity $\ell$ and
the optimal parameter $V$, maximizing the proportional fairness ($\alpha=1$)
of the control $\gamma$, as formulated in (8). In the high traffic regime,
preemptively curbing emissions becomes necessary to avoid future budget
depletion, leading to a decrease in the optimal $V$. Conversely, if the
traffic is sparse ($\ell\approx 0$), it is optimal to curb emissions only when
the budget depletes, resulting in a higher $V$, and the DPP policy tends
towards the greedy one. Since $\beta$ and $V$ have opposite effects on the
value of the control $\gamma$, decreasing $\beta$ makes the policy more
conservative, leading to an increase in the optimal value of $V$. Here, in
practice, one can set $\beta$ to a constant value ($\beta\approx 0.95$) and
adjust $V$ according to variations in traffic load.

Figure 3: Value of parameter $V$ maximizing proportional fairness as the
traffic demand sparsity $\ell$ varies.

##  V Conclusions

We have investigated an actual EIRP control method designed to ensure
compliance with EMF exposure limits established in [1], [2] using the actual
maximum approach described in [3]. Our method ensures that the _actual_ EIRP,
averaged over a sliding time window, remains below the configured threshold.
Additionally, it guarantees a minimum EIRP level at all times, which is
crucial for services like GBR traffic. Furthermore, our method proactively
limits emissions to prevent future EIRP budget depletion, which could lead to
network resource shortages.

To update the EIRP budget over time, we have devised both an exact procedure
with linear complexity (Algorithm 2) and a conservative one (Algorithm 3) with
constant complexity. Our actual EIRP control algorithm, rooted in Drift-Plus-
Penalty theory, preemptively limits EIRP consumption based on the status of a
virtual queue, assessing recent consumption overshoot relative to the actual
EIRP threshold $\overline{C}$. It relies on one main parameter $V$, which can
be adjusted as traffic load profile varies.

## Appendix: Proofs

###  V-A Proposition 1

###### Proof.

The first inequality in (10) stems from (6). Next we fix period $t$ and assume
$\gamma_{t}\leq\Gamma_{t}$. We show that all the constraints (4) containing
period $t$ can be satisfied by _some_ control policy. Set
$\gamma_{t^{\prime}}=\rho\overline{C}$ for all $t^{\prime}>t$. Then, since
$c_{t}\leq\gamma_{t}$,

$\displaystyle c_{t}\leq\Gamma_{t}$ $\displaystyle=\,\\!\\!\\!\min_{0\leq
k\leq
W-1}\,\\!\\!\\!W\overline{C}-\\!\\!\\!\sum_{i=1}^{\min(t,k)}\\!\\!\\!c_{t-i}-(W\\!-\\!k\\!-\\!1)\rho\overline{C}$
(24) $\displaystyle\leq\,\\!\\!\\!\min_{0\leq k\leq
W-1}\\!\\!\\!W\overline{C}-\sum_{i=1}^{\min(t,k)}c_{t-i}-\sum_{i=1}^{W-k-1}c_{t+i}$
(25)

which is an equivalent expression of all constraints in (4) containing period
$t$, Conversely, if $\gamma_{t}>\Gamma_{t}$, then by setting
$\gamma_{t^{\prime}}=\rho\overline{C}$ for all $t^{\prime}>t$, at least one
constraint (4) is violated if $c_{t^{\prime}}=\gamma_{t^{\prime}}$ for all
$t^{\prime}>t$. Then, to satisfy (4) one must set
$\gamma_{t^{\prime}}<\rho\overline{C}$ for some $t^{\prime}>t$, which in turn
violates constraint (6). Thus, the upper bound is tight, q.e.d.. ∎

###  V-B Proposition 2

###### Proof.

We define $a_{t,k}:=\sum_{i=1}^{k}(c_{t-i}-\rho\overline{C})$. By observing
that $a_{t+1,k+1}=a_{t,k}+c_{t}-\rho\overline{C}$, we derive that:

$\displaystyle\Omega_{t+1}^{n+1}=\max_{0\leq k\leq n+1}a_{t+1,k}=$
$\displaystyle\,\max\left(a_{t+1,0},\,\max_{0\leq k\leq n}a_{t+1,k+1}\right)$
$\displaystyle=$ $\displaystyle\,\Big{[}\max_{0\leq k\leq
n}a_{t,k}+c_{t}-\rho\overline{C}\Big{]}^{+}\\!\\!\\!\\!$ (26)

which equals $\left[\Omega_{t}^{n}+c_{t}-\rho\overline{C}\right]^{+}$, q.e.d..
∎

###  V-C Property 1

###### Proof.

If $t<W-1$, then $\ell_{t}:=\ell_{t}^{\min(t,W-1)}=\ell_{t}^{t}<W-1$. It stems
from Proposition 2 that:

$\displaystyle\Omega_{t+1}^{\min(t+1,W-1)}=\Omega_{t+1}^{t+1}=$
$\displaystyle\,\left[\Omega_{t}^{t}+c_{t}-\rho\overline{C}\right]^{+}$ (27)
$\displaystyle=$
$\displaystyle\,\big{[}\Omega_{t}^{\min(t,W-1)}+c_{t}-\rho\overline{C}\big{]}^{+}.$
(28)

Next, let $t\geq W-1$. From Proposition 2 we know that
$\Omega_{t+1}^{W-1}=[\Omega_{t}^{W-2}+c_{t}-\rho\overline{C}]^{+}$. If
$\ell_{t}^{W-1}<W-1$, then

$\displaystyle\Omega_{t}^{W-1}:=$ $\displaystyle\,\max_{0\leq k\leq
W-1}a_{t,k}=\,\max_{0\leq k\leq W-2}a_{t,k}=:\,\Omega_{t}^{W-2}.$ (29)

The thesis follows. ∎

###  V-D Property 2

###### Proof.

We define $a_{t,k}$ as in the proof of Proposition 2. Then,
$\ell_{t+1}^{n+1}=\operatorname*{arg\,max}_{0\leq k\leq n+1}a_{t+1,k}$. If
$\Omega_{t+1}^{n+1}=0$, then $\ell_{t+1}^{n+1}=0$ since $a_{t+1,0}=0$. Else,
if $\Omega_{t+1}^{n+1}>0$, we can write:

$\displaystyle\ell_{t+1}^{n+1}=$
$\displaystyle\,\operatorname*{arg\,max}_{1\leq k\leq
n+1}a_{t,k-1}+c_{t}-\rho\overline{C}=\,\operatorname*{arg\,max}_{1\leq k\leq
n+1}a_{t,k-1}$

which equals $\operatorname*{arg\,max}_{0\leq k\leq
n}a_{t,k}+1=\ell_{t}^{n}+1$, q.e.d.. ∎

###  V-E Property 3

###### Proof.

If $c_{t-i}\geq\rho\overline{C}$ for $0\leq i\leq W-1$, then it stems from (9)
that $\Omega_{t}=\sum_{i=1}^{W-1}(c_{t-i}-\rho\overline{C})$ and
$\Omega_{t+1}=\sum_{i=0}^{W-2}(c_{t-i}-\rho\overline{C})$. By subtracting term
by term we obtain (14), q.e.d.. ∎

###  V-F Proposition 3

###### Proof.

Since (16) holds, it suffices to prove that
$\Omega_{t}\leq\widetilde{\Omega}_{t}$. Indeed, $\Omega_{t}\leq\max_{0\leq
k\leq\min(t,W-1)}\sum_{i=1}^{k}\left[c_{t}-\rho\overline{C}\right]^{+}\\!\\!=\widetilde{\Omega}_{t}$.
∎

###  V-G Proposition 4

###### Proof.

Let $t_{1}:=0<t_{2}<\dots$ be the periods at which $\overline{Q}=0$. Then, we
will prove that $t_{i+1}-t_{i}\leq W$ for all $i\geq 0$. Let us suppose that
$\overline{Q}_{t_{i}}=0$. We define
$\Delta_{i}=\min\\{\delta:\,\sum_{j=0}^{\delta}(c_{t_{i}+j}-\overline{C})\leq
0\\}$. It stems from (4) that $\Delta_{i}<W$. Also,
$\overline{Q}_{t_{i}+\Delta_{i}}=\overline{Q}_{t_{i}}+\sum_{j=0}^{\Delta_{i}-1}(c_{t_{i}+j}-\overline{C})$.
Thus,

$\displaystyle\overline{Q}_{t_{i}+\Delta_{i}+1}=$
$\displaystyle\,\left[\overline{Q}_{t_{i}+\Delta_{i}}+(c_{t_{i}+\Delta_{i}}-\overline{C})\right]^{+}$
(30) $\displaystyle=$
$\displaystyle\,\big{[}\overline{Q}_{t_{i}}+\sum_{j=0}^{\Delta_{i}}(c_{t_{i}+j}-\overline{C})\big{]}^{+}=0.$
(31)

Therefore, $t_{i+1}-t_{i}=\Delta_{i}+1\leq W$, q.e.d.. ∎

## References

  * [1] Int. Commission. on Non-Ionizing Radiation Protection (ICNIRP), “Guidelines for limiting exposure to electromagnetic fields (100 kHz to 300 GHz),” _Health physics_ , vol. 118, no. 5, pp. 483–524, 2020.
  * [2] IEEE C95. 1-2019, “IEEE standard for safety levels with respect to human exposure to electric, magnetic, and electromagnetic fields, 0 Hz to 300 GHz,” _IEEE Std C95. 1-2019 (Revision of IEEE Std C95. 1-2005/Incorporates IEEE Std C95. 1-2019/Cor 1-2019)_ , pp. 1–312, 2019.
  * [3] Int. Electrotechnical Commission, “Determination of RF field strength, power density and SAR in the vicinity of radiocommunication base stations for the purpose of evaluating human exposure,” _IEC 62232:2022_ , 2022.
  * [4] P. Baracca, A. Weber, T. Wild, and C. Grangeat, “A statistical approach for RF exposure compliance boundary assessment in massive MIMO systems,” in _WSA 2018; 22nd International ITG Workshop on Smart Antennas_. VDE, 2018, pp. 1–6.
  * [5] M. Rybakowski, K. Bechta, C. Grangeat, and P. Kabacik, “Impact of Beamforming Algorithms on the Actual RF EMF Exposure From Massive MIMO Base Stations,” _IEEE Access_ , 2023.
  * [6] C. Törnevik, T. Wigren, S. Guo, and K. Huisman, “Time averaged power control of a 4G or a 5G radio base station for RF EMF compliance,” _IEEE Access_ , vol. 8, pp. 211 937–211 950, 2020.
  * [7] D. Ying, D. J. Love, and B. M. Hochwald, “Closed-loop precoding and capacity analysis for multiple-antenna wireless systems with user radiation exposure constraints,” _IEEE Transactions on Wireless Communications_ , vol. 14, no. 10, pp. 5859–5870, 2015.
  * [8] S. Mandelli, L. Maggi, B. Zheng, C. Grangeat, and A. Zejnilagic, “EMF Exposure Mitigation via MAC Scheduling,” _arXiv:2404.06830_ , 2024.
  * [9] J. Mo and J. Walrand, “Fair end-to-end window-based congestion control,” _IEEE/ACM Transactions on networking_ , vol. 8, no. 5, pp. 556–567, 2000\.
  * [10] M. J. Neely, “Stochastic network optimization with application to communication and queueing systems,” _Synthesis Lectures on Communication Networks_ , vol. 3, no. 1, pp. 1–211, 2010.


# Emergence of hot corona and truncated disk in simulations of accreting
stellar mass black holes

Rodrigo Nemmen1,2, Artur Vemado1,3, Ivan Almeida1,4, Javier Garcia5,6 and
Pedro N. Motta1
1Instituto de Astronomia, Geofísica e Ciências Atmosféricas, Universidade de
São Paulo, São Paulo, SP 05508-090, Brazil
2Kavli Institute for Particle Astrophysics and Cosmology (KIPAC), Stanford
University, Stanford, CA 94305, USA
3LEEGA, São Paulo, SP 04729-080, Brazil
4 School of Mathematics, Statistics and Physics, Newcastle University, NE1
7RU, UK
5Cahill Center for Astronomy and Astrophysics, California Institute of
Technology, Pasadena, CA 91125, USA
6NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA E-mail:
<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Stellar mass black holes in X-ray binaries (XRBs) are known to display
different states characterized by different spectral and timing properties,
understood in the framework of a hot corona coexisting with a thin accretion
disk whose inner edge is truncated. There are several open questions related
to the nature and properties of the corona, the thin disk, and dynamics behind
the hard state. This motivated us to perform two-dimensional hydrodynamical
simulations of accretion flows onto a $10M_{\odot}$ black hole. We consider a
two-temperature plasma, incorporate radiative cooling with bremmstrahlung,
synchrotron and comptonization losses and approximate the Schwarzschild
spacetime via a pseudo-Newtonian potential. We varied the mass accretion rate
in the range $0.02\leq\dot{M}/\dot{M}_{\rm Edd}\leq 0.35$. Our simulations
show the natural emergence of a colder truncated thin disk embedded in a hot
corona, as required to explain the hard state of XRBs. We found that as
$\dot{M}$ increases, the corona contracts and the inner edge of the thin disk
gets closer to the event horizon. At a critical accretion rate
$0.02\leq\dot{M}_{\text{crit }}/\dot{M}_{\rm Edd}\leq 0.06$, the thin disk
disappears entirely. We discuss how our simulations compare with XRB
observations in the hard state.

###### keywords:

accretion, accretion discs – black hole physics – stars: black holes –
methods: numerical

††pubyear: 2023††pagerange: Emergence of hot corona and truncated disk in
simulations of accreting stellar mass black holes–7

##  1 Introduction

Stellar mass black holes in binary systems (also known as X-ray binaries,
XRBs) are known to display three states characterized by different spectral
and timing properties (e.g. Remillard & McClintock 2006; Fender & Belloni
2012). In the low or hard state the source is fainter and its X-ray spectrum
displays a hard power-law. As it brightens, it enters a high or soft state, in
which the spectrum is characterized by a thermal component. Between these two
states the source goes through an intermediate (or very high state) state
where the thermal and power-law components are comparable. A black hole walks
though all these states in a couple of months, and sometimes days.

It is widely accepted that the hard and thermal states can be understood
separately in terms of a hot corona and a thin disk, respectively (Esin et
al., 1997; Done et al., 2007). The basic qualitative idea is that, in the hard
state, the thin disk is truncated at a certain transition radius $R_{\rm tr}$
inside of which it gives way to a hot corona (i.e. a radiatively inefficient
accretion flow or RIAF). As the mass accretion rate $\dot{M}$ increases,
$R_{\rm tr}$ becomes progressively smaller until the entire accretion flows
consists of a thin disk: the source transitions to the soft state. Physically,
this is due to the cooling time becoming shorter than the accretion timescale
(Narayan & Yi, 1995a).

Two unresolved questions in XRB research stand out. The first: what is the
nature and properties of the X-ray-emitting corona? The second question: what
is the mechanism and dynamics behind state transitions? Observationally, there
is firm evidence that both a hot corona and a truncated thin disk are present
along the state cycles of XRBs. However, it has been challenging to explain
from first principles how a stable corona is formed and its physical
properties (e.g. Schnittman et al. 2013; Yuan & Narayan 2014; Dexter et al.
2021). In addition, the details of the transition from the hard to soft state
occurs are hotly debated. Is it driven by a decrease in the inner edge of the
thin accretion disk (Ingram & Done, 2011; Plant et al., 2014), or a diminution
of the corona (García et al., 2015; Kara et al., 2019; Wang et al., 2021)?
This remains unclear.

What is clear is that radiative cooling in the accretion flow should play a
major role in explaining these transitions. While historically analytical one-
dimensional models have been very useful in beginning to address these
questions (e.g. Chen et al. 1995; Esin et al. 1997), some aspects of accretion
physics such as their dynamics are beyond the scope of such models.

In this work, we use viscous hydrodynamical simulations designed to mimic the
state of the XRB accretion flow in the hard state. We incorporate energy
losses due to radiative cooling and investigate a range of mass accretion
rates near the critical value $\dot{M}\approx 0.01\dot{M}_{\rm Edd}$ at which
a RIAF collapses to a thin disk. Section 2 describes the numerical experiment
and initial condition of the accretion flow. Section 3 describes the results
of our experiment. Section 4 discusses the results and how they can help to
explain XRB observations. Finally, section 5 summarizes the work. The
appendices describe technical details about the methods.

##  2 Summary of methods

We have adopted the code PLUTO111See the data availability section for the
code with specific customizations used in this work. (Mignone et al., 2007) to
solve the equations of fluid dynamics in spherical coordinates (Vemado 2020
for the full equations and details). We approximate the Schwarzschild
spacetime via the pseudo-Newtonian gravitational potential
$\phi=-GM/(r-r_{\text{s}})$ (Paczyńsky & Wiita, 1980) ($r_{s}$ is the
Schwarzschild radius), as the spin parameter does not seem to be important at
first order to understand XRB state transitions (e.g. Remillard & McClintock
2006) and we are not addressing jet formation.

We assume an $\alpha$-disk, i.e. an accretion flow with shear stress given by
the Shakura-Sunyaev $\alpha$ viscosity (Shakura & Sunyaev, 1973) where the
kinematic viscosity is given by $\nu=\alpha r^{1/2}$ (Landau & Lifshitz, 1959;
Stone et al., 1999). Here, $\alpha$ is a constant similar to (but not equal)
the standard Shakura-Sunyaev $\alpha$. We set $\alpha=0.01$ which is a common
choice (e.g. Wu et al. 2016).

A key aspect of our work are the energy losses via radiative cooling. Below we
give a summary of the heating and cooling implementation; we describe it in
more details in Appendix A. The radiative cooling rate is given by $Q^{-}_{\rm
rad}=Q^{-}_{\rm brem}+Q^{-}_{\rm syn}+Q^{-}_{\rm ssc}$ where we include the
relevant processes for stellar-mass BHs: bremsstrahlung emission $Q^{-}_{\rm
brem}$ including both the electron-ion and electron-electron interactions,
optically thin synchrotron emission $Q^{-}_{\rm syn}$ and comptonized self-
synchrotron emission $Q^{-}_{\rm ssc}$. We adopt a series of approximations
for the calculation of the cooling rates. For the bremsstrahlung emission, we
assume that the proton temperature is small enough to neglect its motion and
the electrons follow a Maxwellian energy distribution. For the synchrotron
emission, we adopt the optically thin limit (Pacholczyk, 1970) and assume that
the angle between the velocity vector of the electron and the direction of the
local magnetic field follows an isotropic distribution (Narayan & Yi, 1995b).
In order to obtain the $Q^{-}_{\rm syn}$, we assume that the observer sees a
blackbody source with radius determined by the condition $\nu=\nu_{c}(r)$ (the
self-absorption critical frequency). We neglect nonlocal radiative transfer
effects. For the self-synchrotron Comptonization, we assume Comptonization of
synchrotron radiation at the frequency $\nu_{c}$. We include the cooling
contribution from higher density, optically thick regions of the flow using
the cooling rate computed by Hubeny (1990).

The radiative cooling term $Q^{-}_{\text{rad}}(r,n_{e},T_{e})$ in the energy
equation depends the distance from the BH $r$, the electron temperature
$T_{e}$ and the electron number density $n_{e}$. We generated a lookup table
of $-Q^{-}_{\text{rad}}$ as a function of $n_{e}$, $r$ and $T_{e}$ for the
following parameter ranges, equally spaced in the log space: $10^{11}\leq
n_{e}/\text{cm}^{-3}\leq 10^{20}$, $2.6\leq r/r_{g}\leq 800$ and $10^{6}\leq
T_{e}/{\rm K}\leq 10^{12}$, encompassing the range of typical XRB values (e.g.
Done et al. 2007). This is much faster than computing cooling rates on-the-
fly.

Distances in the simulations are parameterized in terms of the gravitational
radius $r_{g}\equiv GM/c^{2}=M$. We adopt $M=10M_{\odot}$. The grid resolution
is $400\times 400$ elements ($R\times\theta$), logarithmic in $r$ bounded by
$2.6r_{g}<r<800r_{g}$, linear and non-uniform in $\theta$ with ten cells at
$\theta<15^{\circ}$ $\theta>165^{\circ}$ and the remaining uniformly
distributed, in order to increase resolution towards the equatorial plane. The
initial condition consists of a stationary torus in dynamical equilibrium
(Papaloizou & Pringle, 1984) with a constant specific angular momentum. Its
starting radius is chosen as $r_{\text{min}}=150r_{g}$ and the pressure
maximum radius as $r_{0}=200r_{g}$, resulting in an initial torus extending
all the way to $r_{\text{max}}=300r_{g}$. The mass accretion rate is set by
the torus maximum density $\rho_{0}$, which sets its total mass.

We assume the presence of a local, isotropic, randomly tangled magnetic field
characterized by $\beta=P_{\mathrm{gas}}/P_{\mathrm{mag}}=10$ on the
assumption that the flow is dominated by gas pressure (e.g. Wu et al. 2016).
RIAFs are characterized by a two-temperature plasma (e.g. Yuan & Narayan
2014). We approximate the temperature difference using

$T_{\mathrm{e}}=\frac{T_{i}}{2+(r/r_{0})^{-1}},$ (1)

where $r_{0}$ is the radius of maximum density and the ion temperature $T_{i}$
is taken directly from the simulation. Equation 1 is based on height-
integrated solutions (Narayan & Yi, 1995a; Wu et al., 2016). The effective
mean molecular weights of ions and electrons are $\mu_{i}=1.23$ and
$\mu_{e}=1.14$, corresponding to a hydrogen mass fraction $X=0.75$ (Narayan &
Yi, 1995b).

The beginning of each simulation involves a burn-in period during which the
radiative cooling is disabled until the flow achieves a quasi-steady state
characterized by an almost constant accretion rate at $t\sim 5000GM/c^{3}$.
More precisely, radiative cooling is turned on at $t=1.4t_{\text{orb}}$ where
$t_{\text{orb}}=2\pi r/v_{\phi}\approx 1.8\times 10^{4}GM/c^{3}$ is the
orbital time of the accretion flow at $r=100r_{s}$.

##  3 Results

We carried out five simulations corresponding to different values of the mass
accretion rate $\dot{M}$, which are listed in Table 1. The mass accretion
rates quoted in Table 1 are the time-average of $\dot{M}(r,t)=2\pi
r^{2}\int_{0}^{\pi}\rho v_{\mathrm{r}}\sin{\theta}\,d\theta$ at the inner
boundary; we computed the average for the time intervals over which
$\dot{M}(t)$ is relatively constant thus indicating that the accretion flow
reached steady state, as indicated in Figure1. The accretion rates range from
$\dot{m}=0.02$ up to $\dot{m}=0.35$, where $\dot{m}\equiv\dot{M}/\dot{M}_{\rm
Edd}$, $L_{\rm Edd}\equiv 0.1\dot{M}_{\rm Edd}c^{2}$ and $L_{\rm
Edd}=1.3\times 10^{38}(M/M_{\odot})\ {\rm erg\ s}^{-1}$. Since the critical
value that signals the transition between a RIAF and thin disk is
$\dot{m}_{\rm crit}\approx 0.01$ (Yuan & Narayan, 2014), our simulations
should probe the conditions involved in the hard to thermal state transition
of XRBs.

Figure 1: Temporal evolution of the mass accretion rate after the radiative cooling is enabled. The dotted vertical indicate the time range used to estimate $\dot{m}$ in models 4 and 5 (cf. Table 1); the dashed line shows the range used in the remaining simulations.  ID | $\dot{M}/\dot{M}_{\rm Edd}$ | $L/L_{\rm Edd}$ | $r_{\mathrm{in}}/r_{g}$ | Time interval | $\rho_{0}$
---|---|---|---|---|---
|  |  |  | ($r_{g}/c$) | (${\rm g\ cm}^{-3}$)
1 | $0.35\pm 0.09$ | $0.3\pm 0.08$ | $22\pm 5$ | $35815-40930$ | $1.5\times 10^{-6}$
2 | $0.24\pm 0.07$ | $0.2\pm 0.07$ | $26\pm 5$ | $35815-40930$ | $1\times 10^{-6}$
3 | $0.11\pm 0.02$ | $0.09\pm 0.02$ | $41\pm 10$ | $35815-40930$ | $5\times 10^{-7}$
4 | $0.06\pm 0.01$ | $0.046\pm 0.008$ | $55\pm 12$ | $32400-37520$ | $3\times 10^{-7}$
5 | $0.02$ | $0.015\pm 0.008$ | $>100$ | $32400-37520$ | $1\times 10^{-7}$
Table 1: Time-averaged mass accretion rates, Eddington ratios and thin disk
inner radii for the numerical simulations. The ID column indicates the number
assigned to each model. Column 5 shows the time intervals over which the
averages were computed. Column 6 shows the maximum torus density in CGS units.

A snapshot from the simulations can be seen in Figure 2 where the electron
density and temperature maps are shown. Our simulations demonstrate that for
$\dot{m}\geq 0.06$ the initial torus collapses resulting in the formation of a
cold thin disk embedded in a hot corona, as required to explain XRBs. The
dense disk is located at the equatorial plane and is 10-100 times denser than
the torus in non-radiative simulations. The corona in our models typically has
$T_{e}\sim 10^{9-10}K$ while the disk has typical electron temperatures $\sim
10^{7}K$.

Figure 2: Temperature and density maps of the five simulations carried out in
this work. On the left side, the colors indicate the log of the electron
temperature in Kelvin; on the right side, they indicate the log of electron
number density in particles per cubic cm. All snapshots were taken at
$t=34112GM/c^{3}$.

Over the duration of the simulation, the cold and hot phases coexist until the
accretion rate is decreased to $\dot{m}=0.02$. At this value, we could not
find any trace of a cold disk. The only structure remaining is a geometrically
thick, low-density, hot corona. This is expected because as the accretion rate
gets lower the density also decreases, and the gas starts becoming inefficient
in radiating i.e. the cooling time is longer than the viscous time.

We have also found that the coronal size and the inner radius of the thin disk
change noticeably across the different simulations, in response to the
$\dot{m}$ variation. We estimate the time-averaged coronal scale height
$H_{c}$ and the inner radius $R_{\rm in}$ using the same time range adopted to
compute $\dot{m}$ (Table 1). The uncertainties in these estimates reflect the
standard deviation in the time-averages. For more details on the methods for
obtaining $H_{c}$ and $R_{\rm in}$, please refer to Appendices B and C,
respectively.

Figure 5 shows that the corona contracts as the accretion rate increases.
Quantitatively, as $\dot{m}$ increases sixfold from 0.06 to 0.35 the corona
contracts to about sixty per cent of the size at the lowest $\dot{m}$
simulated. This is reflected in both its lateral width and height. We fitted a
linear relation to the dependence of $H_{c}$ on $\dot{m}$ in log-log space and
found that it is approximately given by

$H_{c}=20r_{g}\left(\frac{\dot{m}}{0.1}\right)^{-0.37\pm 0.08}.$ (2)

Figure 5 also displays the inner radius of the thin disk as a function of the
mass accretion rate. We find that the inner radius of the thin disk is
anticorrelated with $\dot{m}$. In other words, as the accretion rate
increases, the disk gets closer to the event horizon. We fitted a power-law to
our data and obtained

$R_{\rm in}=42r_{g}\left(\frac{\dot{m}}{0.1}\right)^{-0.60\pm 0.27}.$ (3)

Figure 3: Corona vertical height and inner radius of the thin disk as a
function of the mass accretion rate $\dot{m}$. The corona and inner radius
contract as $\dot{m}$ increases. Each line indicate a linear fit in log-log
space using $y|x$ BCES regression with bootstrapping taking into account
uncertainties in both axes (e.g. Nemmen et al. 2012). The band around each fit
denotes the $1\sigma$ prediction band.

##  4 Discussion

Our simulations begin with a hot torus which collapses to a colder thin disk
due to thermal instability. A hot corona persists above and below the thin
disk. The corona and the inner radius of the thin disk gradually increase in
size as $\dot{m}$ decreases. At $\dot{m}=0.02$ the thin disk disappears and
there is only a hot accretion flow left. Thus, we constrain the critical
accretion rate at which the thin disk disappears entirely to be in the range
$0.02<\dot{m}_{\text{crit }}\lesssim 0.06$. This interval is consistent with
estimates based on analytical theory which find values in the range
$0.01\lesssim\dot{m}_{\text{crit }}\lesssim 0.1$ (Narayan & Yi, 1995b; Xie &
Yuan, 2012). In one of their radiation general relativistic
magnetohydrodynamic (GRRMHD) simulations (model M3), Dexter et al. (2021) find
tentative collapse of the hot accretion flow at $\dot{m}\sim 10^{-3}$; this is
one order of magnitude lower than our $\dot{m}_{\text{crit}}$.

###  4.1 Comparison with observations

It is worth comparing our values of $R_{\rm in}$ with measurements from
observations of X-ray binaries. Since observers quote the source luminosity
while our simulations are parameterized in terms of $\dot{m}$, we estimate $L$
from $\dot{m}$ for our models using equation $L=\eta\dot{M}c^{2}$ with the
efficiency of conversion of accreted rest-mass energy into radiation given by
$\eta\equiv\eta_{\rm RIAF}+\eta_{\rm thin}$. $\eta_{\rm RIAF}=\eta_{\rm
RIAF}(\dot{m})$ is based on the RIAF model of Xie & Yuan (2012) assuming an
electron turbulent heating efficiency of $\delta=0.1$; $\eta_{\rm
thin}=\eta_{\rm thin}(R_{\rm in},R_{\rm out})$ corresponds to a truncated thin
disk given by equation 4 of Reb et al. (2018) where we fix $R_{\rm
out}=10^{4}r_{g}$.

Since it is widely agreed that in the thermal state the disk is truncated at
or near the innermost stable circular orbit radius, here we compare our
results with XRBs in the hard state. Figure 4 compares the $R_{\rm in}$ values
as a function of $L/L_{\rm Edd}$ from our simulations with observations of GX
339-4 in the hard state. Figure 4 illustrates the disk truncation controversy:
the inferred inner edge of the disk varies widely among different techniques
even for the same object. One of the methods for measuring $R_{\rm in}$ is
based on fitting the disk reflection component observed with spectroscopy
(e.g. García et al. 2015). Two other methods are fitting the continuum
spectrum of the accretion flow (McClintock et al., 2014) and timing analysis
of X-ray lags (e.g. De Marco et al. 2015). The observations form two separate
tracks in Figure 4: the upper track is roughly consistent with $R_{\rm in}\sim
74r_{g}(\lambda/0.1)^{-0.9\pm 0.2}$ ($\lambda\equiv L/L_{\rm Edd}$) while the
lower track follows $R_{\rm in}\sim 2r_{g}(\lambda/0.1)^{-0.6\pm 0.1}$.

Figure 4: Inner radius of the thin accretion disk as a function of the
luminosity in Eddington units. $R_{\rm in}$ is in units of the gravitational
radius. Circles indicate the results from our numerical simulations. Squares,
triangles and hexagons denote measurements based on observations of X-ray
binary GX 339-4 in the hard state (Reis et al., 2008; Miller et al., 2006;
Petrucci et al., 2014; Kolehmainen et al., 2014; García et al., 2015; Plant et
al., 2015; Basak & Zdziarski, 2016; De Marco et al., 2015; Wang-Ji et al.,
2018). The stars show results from the simulations of Takahashi et al. (2016)
and Liska et al. (2022).

From our simulations we obtain $R_{\rm in}\approx
40r_{g}(\lambda/0.1)^{-0.5}$. In other words, the simulations are naturally
predicting a radius-luminosity trend consistent with observations—namely,
$R_{\rm in}$ approaches the BH as the luminosity increases. While the
characteristic simulated $R_{\rm in}$ values are closer to the upper track,
the slope of the radius-luminosity correlation is consistent the lower track.
One common aspect of all the studies predicting the largest truncation (i.e.
the upper track) is that the observations were taken with XMM-Newton’s Epic PN
camera in timing mode, which is prone to calibration issues and pile-up
effects at the very high count rates reached by the source and could
potentially bias their results (García et al., 2015).

Our results indicate that with increasing Eddington ratios the corona
contracts but does not go away, even at luminosities typical of the soft state
($0.1<L/L_{\rm Edd}<0.4$; cf. Table 1). One possibility is that the compact
corona seen in our models at such luminosities could be responsible for the
hard X-ray tail observed in the soft state of several sources (e.g. Sridhar et
al. 2020; Connors et al. 2021). More work is needed to probe the stability
properties and electromagnetic spectrum of this remnant corona.

###  4.2 Caveats

The specific values of $R_{\rm in}$ and $L/L_{\rm Edd}$ from the models should
be taken with care due to the many shortcomings of our numerical simulations
described below.

We compute the synchrotron cooling rate based on a single frequency—the
critical frequency $\nu_{c}$ of self-absorbed synchrotron emission; $\nu_{c}$
is also used as the seed photon frequency for comptonization; comptonization
is assumed to happen only locally. These simplifications were adopted in order
to reduce the computational costs and as a result we could be underpredicting
the amount of radiative cooling.

Since we performed only hydrodynamical simulations, the dynamical effect of
magnetic fields is only incorporated via the assumption of an isotropically-
tangled magnetic field in equipartition with the gas, with a constant $\beta$.
MHD simulations show that it is not constant as assumed here (e.g. Hirose et
al. 2004). Furthermore, viscous stresses are incorporated using an
$\alpha$-viscosity prescription proportional to the density. This means that
most of the dissipation takes place at the equatorial plane which concentrates
most of the flowing mass. MHD simulations indicate that the vertical
distribution of energy dissipation can extend to higher heights (e.g. Huang et
al. 2023).

At accretion rates $\dot{m}\gtrsim 0.1$, radiation pressure becomes important
and is not incorporated in this work. In addition, gravity in our model is
pseudo-Newtonian. This is very good approximation for the Schwarzschild
spacetime. However, the innermost stable circular orbit is closer to the
horizon in the Kerr spacetime for $a_{*}>0$ ($a_{*}$ is the BH spin) compared
to the Schwarzschild spacetime. Indeed, the ISCO around a BH spinning at
$a_{*}=0.998$ is about five times smaller compared to a nonrotating hole. We
suspect that the Schwarzschild spacetime introduces a positive bias in the
values of $R_{\rm in}$, which could explain why our values are systematically
larger than the measured ones.

Our coronal temperatures reach $10^{9-10}K$, hotter than XRB corona
measurements (e.g. García et al. 2015; Chakraborty et al. 2020), suggesting we
might be missing important physics related to cooling. All the above effects
can impact the results, particularly the corona size, thin disk inner radius,
their dependence on $\dot{m}$ and the $\dot{m}_{\rm crit}$.

XRBs display a hysteresis loop in the hardness-luminosity diagram in which
they move from the hard to soft state along a higher luminosity than the soft
to hard (e.g. Remillard & McClintock 2006). This would correspond to a loop in
the $R_{\rm in}-L/L_{\rm Edd}$ diagram, which our simulations are presently
unable to account for. We think additional physics is required beyond the
dependence of cooling on $\dot{m}$. Different causes have been proposed for
the observed hysteresis, including the interplay between advected magnetic
fields and turbulence (Begelman, 2014), warped disks (Nixon & Salvesen, 2014),
the relation between $\alpha$ and the magnetic Prandtl number (Balbus & Henri,
2008) and the disk magnetization (Petrucci et al., 2008).

Our simulations have a duration of at most $t_{\rm sim}\approx 4\times
10^{4}r_{g}/c$. In order to achieve inflow equilibrium, $t_{\rm sim}$ would
need to be of the order of the viscous time for the thin disk (e.g. Frank et
al. 2002) which scaled to the relevant XRB parameters becomes

$\displaystyle t_{\rm visc}(M,\dot{m},r)=3\times
10^{6}\frac{r_{g}}{c}\left(\frac{\alpha}{0.1}\right)^{-0.8}\left(\frac{M}{10M_{\odot}}\right)^{1.2}\times$
$\displaystyle\left(\frac{r}{10r_{g}}\right)^{1.25}\left(\frac{\dot{m}}{0.1}\right)^{-0.3}.$
(4)

At the radius where most of action is happening in our models, $r\approx
50r_{g}$, $t_{\rm visc}$ is about $8\times 10^{6}r_{g}/c$ (or 15 minutes in
physical units). Therefore, our models do not achieve inflow equilibrium. We
would need to evolve them for durations at least a hundred times longer in
order to achieve a steady state solution. For this reason, it is difficult to
ascertain the robustness of our measured truncation radii. As far as we know,
no global simulation in the literature currently has a long enough duration to
ensure that thin accretion disks achieve inflow equilibrium. Future
simulations—with much longer durations to reach $t_{\rm sim}\sim t_{\rm visc}$
and correspondingly larger tori such that the BH does not run out of gas—are
warranted.

###  4.3 Comparison with other simulations

We now compare our work with other numerical simulations of truncated
accretion disks. Hogg & Reynolds (2017) performed MHD simulations of the
transition state of XRBs with a $-1/r$ gravitational potential. Their initial
condition is a thin disk with an ad-hoc cooling function. Due to their choice
of cooling function, their simulation gives $R_{\rm in}\approx 90r_{g}$ at
$\dot{m}=0.01$ by construction. Other numerical works have improved on many of
the issues we pointed previously. Takahashi et al. (2016) performed GRRMHD
simulations of truncated discs around black holes including radiation pressure
in the M1 closure. For $\dot{m}\approx 0.1$, Takahashi et al. find $R_{\rm
in}\approx 30r_{g}$ which is consistent with our results within uncertainties.
Dexter et al. (2021) performed GRRMHD simulations of a MAD flow onto a Kerr BH
with frequency-dependent Monte Carlo radiative transfer, including synchrotron
and bremsstrahlung emission. In one of their simulations, M3, they reported
the tentative formation of a thin disc at $\dot{m}\sim 10^{-3}$. Dexter et al.
did not report whether their collapsed disk is truncated at the inner edge.
Finally, Liska et al. (2022) simulated a XRB with a GRRMHD formalism including
electron thermodynamics and M1 closure, with an initial condition broadly
similar to ours. They also find a two-phase flow at $\dot{m}=0.2$ with the
thin disk truncated at $R_{\rm in}\approx 20r_{g}$. Liska et al. argue that
the inner edge of their radiation pressure-supported thin disk is dictated by
the onset of a magnetically dominated flow.

The comparison of our results with other numerical simulations suggests that
the formation of a thin disk sandwiched by a hot corona is a surprisingly
robust result. Furthermore, despite their vastly different complexity, the
GRRMHD models of Takahashi et al. (2016); Liska et al. (2022) produce
truncated thin disks whose inner edge locations are consistent with the
results of our models within $1\sigma$ uncertainties (cf. Figure 4). The cause
of truncation, however, differs among the simulations. In ours—which does not
include magnetic fields—it is thermal instability; in Liska et al. it is due
to a magnetic fields. There are different routes to truncation; it remains to
be seen which one is taken by nature.

##  5 Summary and conclusions

We performed two-dimensional hydrodynamical simulations of two-temperature
accretion flows onto a $10M_{\odot}$ black hole, incorporating radiative
cooling from bremmstrahlung, synchrotron and comptonized synchrotron to energy
losses. Schwarzschild spacetime is mimicked using a pseudo-Newtonian
potential. Our motivation was to investigate the hard state of X-ray binaries
and the hard-to-soft state transition. As such, we varied the mass accretion
rate in the range $0.02\leq\dot{M}/\dot{M}_{\rm Edd}\leq 0.35$ near the
critical value at which a RIAF is thought to collapse into a radiatively
efficient disk. Our main results are:

(i) Our simulations result in the formation of a cold thin disk embedded in a
hot corona, as required to explain the hard state of XRBs. The typical
electron temperatures are $10^{9-10}K$ for the corona and $10^{7}K$ for the
thin disk.

(ii) By varying the mass accretion rate in our simulations, we find that the
corona contracts as $\dot{m}$ increases following the relation
$H_{c}=20r_{g}(\dot{m}/0.1)^{-0.37}$ where $H_{c}$ is the corona scale height
and $\dot{m}$ is in Eddington units.

(iii) The inner edge of the thin disk approaches the event horizon as
$\dot{m}$ increases following the relation $R_{\rm
in}=42r_{g}\left(\dot{m}/0.1\right)^{-0.60}$.

(iv) At $\dot{m}=0.02$ the thin disk disappears and there is only a hot
accretion flow left. We constrain the critical accretion rate at which the
thin disk disappears entirely to be in the range $0.02<\dot{m}_{\text{crit
}}\lesssim 0.06$.

These new results are important because they shed light on the details of the
hard state and the hard-to-soft state transition. When a stellar BH emerges
from quiescence it has a hard X-ray spectrum produced by a centrally located
and compact hot corona, which then changes to a soft spectrum dominated by the
thin disk emission (e.g. Done et al. 2007; Fender & Belloni 2012). Much debate
persists over how this transition occurs. In particular, it is not clear
whether the transition is driven by a reduction in the inner radius of the
disk or by a shrinking corona (Plant et al., 2014; Fabian et al., 2014; García
et al., 2015). For example, the X-ray observations of MAXI J1820+070 indicate
that the culprit is the reduction in the spatial extent of the corona (Kara et
al., 2019). Contrary to simple expectations, however, our models suggest that
the truncation radius decreases and the corona shrinks as the XRB evolves
during the hard state raising.

We find that the thin disk is embedded in the hot corona rather than the
latter being above the disk as in a simple “lamp post” picture (Done et al.,
2007). Our values of $R_{\mathrm{in}}$—much like those obtained in current
GRRMHD models—are systematically larger than many reflection spectroscopy XRB
measurements. As BH accretion models become increasingly more sophisticated,
it remains to be seen whether this discrepancy between theory and observations
will be resolved, and which of the two mechanisms—thermal or magnetic—are
responsible for the disk truncation.

## Data Availability

Please refer to https://github.com/black-hole-group/pluto, branch cooling,
commit ce495cc for the specific code and radiative cooling implementations
used in this work. The initial condition files are available at
https://bitbucket.org/nemmen/pluto-simulations in the directory
hd/torus_2D_pn_cool. The simulation binary data files, cooling table and
analysis notebooks used to generate the figures in this work are available on
Figshare at http://tiny.cc/xrb-hydro-figshare.

## Acknowledgements

We acknowledge useful discussions with Roger Blandford, Dan Wilkins, Elisabete
de Gouveia Dal Pino, Raimundo Lopes de Oliveira Filho, Sergio Martin Alvarez,
Chris Done and Andrzej Zdziarski. We are grateful to the anonymous referee for
the insightful comments and suggestions, which enhanced the quality of this
paper. This work was supported by the following agencies: FAPESP (Fundação de
Amparo à Pesquisa do Estado de São Paulo) through grants 2017/25710-1,
2017/01461-2, 2019/10054-7, 2021/01563-5, 2022/07456-9 and 2022/10460-8; NASA
through the Fermi Guest Investigator Program (Cycle 16). Research developed
with the help of HPC resources provided by the Information Technology
Superintendence of the University of São Paulo. IA acknowledge funding from an
United Kingdom Research and Innovation grant (code: MR/V022830/1). The Black
Hole Group used GPUs generously donated by NVIDIA under the GPU Grant Program.

## References

  * Balbus & Henri (2008) Balbus S. A., Henri P., 2008, ApJ, 674, 408
  * Basak & Zdziarski (2016) Basak R., Zdziarski A. A., 2016, MNRAS, 458, 2199
  * Begelman (2014) Begelman M. C., 2014, preprint,  (arXiv:1410.8132)
  * Chakraborty et al. (2020) Chakraborty S., Navale N., Ratheesh A., Bhattacharyya S., 2020, MNRAS, 498, 5873
  * Chen et al. (1995) Chen X., Abramowicz M. A., Lasota J.-P., Narayan R., Yi I., 1995, ApJ, 443, L61
  * Connors et al. (2021) Connors R. M. T., et al., 2021, ApJ, 909, 146
  * De Marco et al. (2015) De Marco B., Ponti G., Muñoz-Darias T., Nandra K., 2015, ApJ, 814, 50
  * Dermer et al. (1991) Dermer C. D., Liang E. P., Canfield E., 1991, ApJ, 369, 410
  * Dexter et al. (2021) Dexter J., Scepi N., Begelman M. C., 2021, ApJ, 919, L20
  * Done et al. (2007) Done C., Gierliński M., Kubota A., 2007, A&ARv, 15, 1
  * Esin et al. (1997) Esin A. A., McClintock J. E., Narayan R., 1997, ApJ, 489, 865
  * Fabian et al. (2014) Fabian A. C., Parker M. L., Wilkins D. R., Miller J. M., Kara E., Reynolds C. S., Dauser T., 2014, MNRAS, 439, 2307
  * Fender & Belloni (2012) Fender R., Belloni T., 2012, Science, 337, 540
  * Frank et al. (2002) Frank J., King A., Raine D. J., 2002, Accretion Power in Astrophysics, 3rd edn. Cambridge University Press
  * García et al. (2015) García J. A., Steiner J. F., McClintock J. E., Remillard R. A., Grinberg V., Dauser T., 2015, ApJ, 813, 84
  * Hirose et al. (2004) Hirose S., Krolik J. H., De Villiers J.-P., Hawley J. F., 2004, ApJ, 606, 1083
  * Hogg & Reynolds (2017) Hogg J. D., Reynolds C. S., 2017, ApJ, 843, 80
  * Huang et al. (2023) Huang J., Jiang Y.-F., Feng H., Davis S. W., Stone J. M., Middleton M. J., 2023, ApJ, 945, 57
  * Hubeny (1990) Hubeny I., 1990, ApJ, 351, 632
  * Ingram & Done (2011) Ingram A., Done C., 2011, MNRAS, 415, 2323
  * Kara et al. (2019) Kara E., et al., 2019, Nature, 565, 198
  * Kolehmainen et al. (2014) Kolehmainen M., Done C., Díaz Trigo M., 2014, MNRAS, 437, 316
  * Landau & Lifshitz (1959) Landau L. D., Lifshitz E. M., 1959, Fluid mechanics
  * Liska et al. (2022) Liska M. T. P., Musoke G., Tchekhovskoy A., Porth O., Beloborodov A. M., 2022, ApJ, 935, L1
  * McClintock et al. (2014) McClintock J. E., Narayan R., Steiner J. F., 2014, Space Science Reviews, 183, 295
  * Mignone et al. (2007) Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, ApJS, 170, 228
  * Miller et al. (2006) Miller J. M., Homan J., Steeghs D., Rupen M., Hunstead R. W., Wijnands R., Charles P. A., Fabian A. C., 2006, ApJ, 653, 525
  * Narayan & Yi (1995a) Narayan R., Yi I., 1995a, ApJ, 444, 231
  * Narayan & Yi (1995b) Narayan R., Yi I., 1995b, ApJ, 452, 710
  * Nemmen et al. (2012) Nemmen R. S., Georganopoulos M., Guiriec S., Meyer E. T., Gehrels N., Sambruna R. M., 2012, Science, 338, 1445
  * Nixon & Salvesen (2014) Nixon C., Salvesen G., 2014, MNRAS, 437, 3994
  * Pacholczyk (1970) Pacholczyk A. G., 1970, Radio astrophysics. Nonthermal processes in galactic and extragalactic sources
  * Paczyńsky & Wiita (1980) Paczyńsky B., Wiita P. J., 1980, A&A, 88, 23
  * Papaloizou & Pringle (1984) Papaloizou J. C. B., Pringle J. E., 1984, MNRAS, 208, 721
  * Petrucci et al. (2008) Petrucci P.-O., Ferreira J., Henri G., Pelletier G., 2008, MNRAS, 385, L88
  * Petrucci et al. (2014) Petrucci P. O., Cabanac C., Corbel S., Koerding E., Fender R., 2014, A&A, 564, A37
  * Piran et al. (2015) Piran T., Sądowski A., Tchekhovskoy A., 2015, MNRAS, 453, 157
  * Plant et al. (2014) Plant D. S., Fender R. P., Ponti G., Muñoz-Darias T., Coriat M., 2014, MNRAS, 442, 1767
  * Plant et al. (2015) Plant D. S., Fender R. P., Ponti G., Muñoz-Darias T., Coriat M., 2015, A&A, 573, A120
  * Reb et al. (2018) Reb L., Fernández-Ontiveros J. A., Prieto M. A., Dolag K., 2018, MNRAS, 478, L122
  * Reis et al. (2008) Reis R. C., Fabian A. C., Ross R. R., Miniutti G., Miller J. M., Reynolds C., 2008, MNRAS, 387, 1489
  * Remillard & McClintock (2006) Remillard R. A., McClintock J. E., 2006, ARA&A, 44, 49
  * Schnittman et al. (2013) Schnittman J. D., Krolik J. H., Noble S. C., 2013, ApJ, 769, 156
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, A&A, 24, 337
  * Sridhar et al. (2020) Sridhar N., García J. A., Steiner J. F., Connors R. M. T., Grinberg V., Harrison F. A., 2020, ApJ, 890, 53
  * Stepney & Guilbert (1983) Stepney S., Guilbert P. W., 1983, MNRAS, 204, 1269
  * Stone et al. (1999) Stone J. M., Pringle J. E., Begelman M. C., 1999, MNRAS, 310, 1002
  * Takahashi et al. (2016) Takahashi H. R., Ohsuga K., Kawashima T., Sekiguchi Y., 2016, ApJ, 826, 23
  * Vemado (2020) Vemado A., 2020, Master’s thesis, Universidade de Sao Paulo, https://doi.org/10.11606/D.14.2020.tde-14052020-172856
  * Wang-Ji et al. (2018) Wang-Ji J., et al., 2018, The Astrophysical Journal, 855, 61
  * Wang et al. (2021) Wang J., et al., 2021, ApJ, 910, L3
  * Wu et al. (2016) Wu M.-C., Xie F.-G., Yuan Y.-F., Gan Z., 2016, MNRAS, 459, 1543
  * Xie & Yuan (2012) Xie F.-G., Yuan F., 2012, MNRAS, 427, 1580
  * Yuan & Narayan (2014) Yuan F., Narayan R., 2014, ARA&A, 52, 529

##  Appendix A Heating and cooling

Here, we describe technical details of the heating and cooling implementation.
Two source terms need to be added to the equations in order to incorporate (i)
heating and angular momentum removal due to magnetically-driven turbulent
viscosity and (ii) energy losses due to radiative cooling.

We assume an $\alpha$-disk (Shakura & Sunyaev, 1973) and adopt the “K-model”
viscosity prescription in Stone et al. (1999), which corresponds to the
kinematic viscosity $\nu=\alpha r^{1/2}$ where $\alpha$ is a constant similar
to (but not equal) the standard Shakura-Sunyaev $\alpha$. We set $\alpha=0.01$
which is a common choice (e.g. Wu et al. 2016).

The radiative cooling rate is given by

$Q^{-}_{\rm rad}=Q^{-}_{\rm brem}+Q^{-}_{\rm syn}+Q^{-}_{\rm ssc}$ (5)

where we include the relevant processes for stellar-mass BHs: bremsstrahlung
emission $Q^{-}_{\rm brem}$ including both the electron-ion and electron-
electron interactions, optically thin synchrotron emission $Q^{-}_{\rm syn}$
and comptonized self-synchrotron emission $Q^{-}_{\rm ssc}$. We adopt a series
of approximations for the calculation of the cooling rates, as follows.

###  A.1 Bremsstrahlung

The bremsstrahlung emission is given by

$Q^{-}_{\rm brem}=Q^{-}_{ei}+Q^{-}_{ee}$ (6)

where the subscript $ei$ and $ee$ denote the electron-ion and the electron-
electron rates. In this approach, we assume that the proton temperature is
small enough to neglect its motion and the electrons follow a Maxwellian
energy distribution. The spectral emissivity is integrated and the fit gives
(Stepney & Guilbert, 1983)

$Q^{-}_{ei}=1.48\times 10^{-22}n_{e}^{2}\times F(\theta_{e})\ {\rm erg\
cm}^{-3}{\rm s}^{-1},$ (7)

where $F(\theta_{e})$ is given by different expressions depending on the value
of $\theta_{e}\equiv kT_{e}/m_{e}c^{2}$:

$F(\theta_{e})=\begin{cases}\frac{9\theta_{e}}{2\pi}\left[\ln\left(1.123\theta_{e}+0.48\right)+1.5\right]&\text{if
}\theta_{e}<1,\\\
4\left(\frac{2\theta_{e}}{\pi^{3}}\right)^{1/2}\left(1+1.781\theta_{e}^{1.34}\right)&\text{if
}\theta_{e}>1.\end{cases}$ (8)

For the electron-electron rate, Stepney & Guilbert (1983) gives the photon
spectrum and integrate it with an approximated cross-section, resulting in

$Q^{-}_{ee}=\begin{cases}2.56\times
10^{-22}n_{e}^{2}\theta_{e}^{3/2}\left(1+1.1\theta_{e}+\theta_{e}^{2}-1.25\theta_{e}^{5/2}\right)&\text{if
}\theta_{e}<1,\\\ 3.40\times
10^{-22}n_{e}^{2}\theta_{e}\left(\ln{1.123\theta_{e}}+1.28\right)&\text{if
}\theta_{e}>1,\end{cases}$ (9)

in units of ${\rm erg\ cm}^{-3}{\rm s}^{-1}$. With these approximations, the
formulas above incur in errors of less than 5 per cent (Stepney & Guilbert,
1983). Some of the constants here are different from Stepney & Guilbert (1983)
to ensure smoothness across $\theta_{e}=1$ (Narayan & Yi, 1995b).

###  A.2 Synchrotron

For the synchrotron emission, we adopt the optically thin limit with a
relativistic Maxwellian distribution of electrons (Pacholczyk, 1970). We
assume that the angle between the electron velocity vector and the direction
of the local magnetic field is isotropically distributed; this incurs an error
of less than 2 per cent (Narayan & Yi, 1995b). We estimate the self-absorption
critical frequency $\nu_{c}$ by assuming that at each radius $r$, the
synchrotron emission occurs over a volume $2H\Delta s$ where $H$ is the
vertical scale height ($H=Rc_{s}/v_{K}$), $c_{s}$ is the sound speed, $v_{K}$
is the Keplerian velocity and $\Delta s$ is the surface area following Wu et
al. (2016). Matching the synchrotron and Rayleigh-Jeans blackbody emission
results in a transcendental equation for $x_{m}$,

$\displaystyle\exp{\left(1.8899x_{m}^{1/3}\right)}=2.49\times
10^{-10}\frac{12\pi
n_{e}H}{B}\frac{1}{\theta_{e}^{3}K_{2}(1/\theta_{e})}\times$ (10)
$\displaystyle\left(\frac{1}{x_{m}^{7/6}}+\frac{0.40}{x_{m}^{17/12}}+\frac{0.5316}{x_{m}^{5/3}}\right),$
(11)

where $K_{2}$ is the modified Bessel’s function. This equation is solved at
each $H$ (i.e. each $r$) to obtain $x_{m}$. It is then used to calculate
$\nu_{c}$ from

$\nu_{c}=\frac{3}{2}\nu_{o}\theta_{e}^{2}x_{m}$ (12)

where $\nu_{o}=2.8\times 10^{6}B$ Hz with $B$ is Gauss.

We obtain the cooling rate per unit volume by matching the total cooling over
a shell with the net flux reaching an observer at infinity and assuming that
at each frequency $\nu$ the observer sees a blackbody with radius determined
by the condition $\nu=\nu_{c}(r)$. This results in

$Q^{-}_{\rm
syn}=\frac{2\pi}{3c^{2}}kT_{e}(r)\frac{d\nu_{c}^{3}(r)}{dr}\approx\frac{2\pi}{3c^{2}}kT_{e}(r)\frac{\nu_{c}^{3}(r)}{r}.$
(13)

Equation 13 is written in a form that explicitly ensures that the integral of
$Q^{-}_{\rm syn}$ over the entire flow is equal to the total cooling radiation
reaching infinity.

###  A.3 Comptonization

We neglect nonlocal radiative transfer effects, i.e. photons being generated
in one place and scattered in another. The Compton cooling is mainly
characterized by the Comptonized energy enhancement factor $\eta$. This
parameter is defined as the average change in energy of a photon between
injection and escape. An approximate description for $\eta$ is given by Dermer
et al. (1991) as

$\eta=1+\frac{P(A-1)}{1-PA}\left[1-\left(\frac{x}{3\theta_{e}}\right)^{-1-\frac{\ln{P}}{\ln{A}}}\right]$
(14)

where $x\equiv h\nu/(m_{e}c^{2})$, $P$ is the probability of photon
scattering,

$P=1-e^{-\tau_{\rm s}},$ (15)

$\tau_{\rm s}=2n_{e}\sigma_{T}H$ is the scattering optical depth and $A$ is
the mean amplification factor in the energy of a scattered photon, considering
a Maxwellian velocity distribution of temperature $\theta_{e}$ for the
electrons,

$A=1+4\theta_{e}+16\theta_{e}^{2}.$ (16)

In the synchrotron cooling model, the escaping radiation is emitted mostly at
the local self-absorption cutoff frequency $\nu_{c}$. The Comptonization of
this radiation gives an additional cooling rate of

$Q^{-}_{\rm ssc}=Q^{-}_{\rm syn}\left[\eta(\nu_{c},\tau_{s})-1\right].$ (17)

###  A.4 Optically thick regime

Since we have the formation of dense regions in our simulations, it is
important to compute the radiative cooling of the higher density, colder,
optically thick regions of the accretion flow associated with the thin disk.
We adopt the net cooling rate of Hubeny (1990),

$Q^{-}_{\rm rad}=\frac{4\sigma
T_{e}^{4}}{H}\frac{1}{3\tau/2+\sqrt{3}+1/\tau_{\rm abs}}.$ (18)

Here, $\tau=\tau_{\rm abs}+\tau_{s}$ is the total optical depth and $\tau_{\rm
abs}$ can be approximated as

$\tau_{\rm abs}=\frac{H}{4\sigma T_{e}^{4}}\left(Q^{-}_{\rm brem}+Q^{-}_{\rm
syn}+Q^{-}_{\rm ssc}\right).$ (19)

(cf. equation 3.33 from Narayan & Yi 1995b). Note that this expression is
valid for both optically thick and thin limits, because when $\tau\gg 1$,
equation 18 results in the appropriate blackbody limit. On the other hand,
when $\tau\ll 1$, equation 18 becomes equation 5.

##  Appendix B Corona thickness

Before quantifying the coronal vertical thickness $H_{c}$, we first convert
the array with the numerical values of density from the native polar
coordinates adopted in PLUTO to cylindrical coordinates $R-z$—a procedure we
call regridding for which GPU-accelerated routines written by the group are
used for faster calculations222The routines we developed for analyzing PLUTO
dump files are available at https://bitbucket.org/nemmen/mickey/.. The density
values are then obtained as a function of $R$.

At each value of $R$, the distribution of density as function of $z$ is a
superposition of two Gaussians centered at $z=0$, the disk midplane: one with
low amplitude and small width corresponding to the thin disk (if it forms),
and a second broader Gaussian with large amplitude describing the RIAF.

We compute the time-averaged density map for each model within the time
intervals displayed in Table 1. Given the time-averaged density map
$\bar{\rho}_{t}=\rho(R,z)$, for the points where the condition
$\log{(\rho/\rho_{0})}\geq-0.5$ is valid the density is replaced by the
threshold value, i.e. $\log\rho\rightarrow-0.5$. This procedure is necessary
because as the thin disk is much denser than the corona, the distribution of
densities of the thin disk will dominate in the estimates of vertical
thickness. Therefore, the denser part corresponding to the thin disk must be
removed. One example of thin disk removal is displayed in Figure 5.

Figure 5: Result of removing the thin disk for estimating the corona
properties. Here, the time-averaged density of model 3 is displayed.

At each value of $R$ we fit a Gaussian function to the $\rho(z)$ distribution,
obtaining a list of accretion flow scale heights $H^{\prime}(R)=\sigma(R)$
where $\sigma$ is the Gaussian width at a given $R$. Finally, we estimate the
coronal thickness as $H_{c}=2\max[H^{\prime}(R)]$.

##  Appendix C Inner radius of thin disk

In order to estimate the inner radius of the thin disk $R_{\rm in}$, we
perform the same coordinate transformation described in Appendix B and obtain
the spacetime density distribution $\rho=\rho(R,z;t)$. At each value of $R$,
we fit a Gaussian function in order to characterize the vertical density
distribution. The Gaussian width corresponds to the accretion flow scale
height, $H^{\prime}=\sigma$—not to be confused with the coronal thickness
$H_{c}$. At each simulation time step, this gives us $H^{\prime}/R$ as a
function of $R$ (Figure 6).

When the flow collapses into a thin disk, the density increases by a factor of
$10-100$ times the RIAF equatorial density. The Gaussian width in this case
corresponds to the vertical density distribution of the thin disk instead of
the RIAF. We define the thin accretion disk as the cold and dense gas with
$H^{\prime}/R<0.015$. This is consistent with the range $10^{-3}\lesssim
H/R\lesssim 0.01$ expected in accretion disk theory $10^{-3}-0.01$ (Piran et
al., 2015; Frank et al., 2002). The thin disks in our simulations settle to
vertical equilibrium with $H^{\prime}/R>4\times 10^{-3}$. Our results are not
very sensitive to the particular value of the threshold thickness as long as
it is not much larger than $0.01$.

We notice in Figure 6 that the transition between the hot accretion flow
($H^{\prime}/R\approx 0.15$) and the cold thin disk ($H^{\prime}/R\lesssim
0.01$) is not sharp. We define $R_{\rm in}$ as the value of $R$ at which the
condition $H^{\prime}/R<0.015$ is attained based on the moving average
computed using 80 points starting from smaller $R$. This definition ensures
that the inner radius is located where the oscillations in the value of
$H^{\prime}/R$ reaches a minimum (cf. dashed line in Figure 6).

Figure 6: An example of the dependency of the disk thickness $H^{\prime}/R$ on
the cylindrical radius $R/R_{\mathrm{s}}$, where $R_{s}$ is the Schwarzschild
radius. The dashed line indicates the resulting inner radius based on the
method outlined in the text.

Figure 7 shows a scatter plot of $R_{\rm in}$ versus $\dot{m}$ for the four
models where a thin disk emerged. This figure illustrates the variability of
the inner radius. For a given model, the final inner radius is computed as the
median of $R_{\rm in}$ using the time range indicated in Table 1. The
uncertainty that we quote in our $R_{\rm in}$ estimates correspond to the
standard deviation.

Figure 7: Distribution of inner radii for the simulations where a thin disk
emerged, as a function of the mass accretion rate in Eddington units. Each
point corresponds to the value of $R_{\rm in}$ computed using the procedure
outlined in Appendix C for a given timestep.