11institutetext: College of Mathematics and Physics, Hunan University of Arts
and Sciences, 3150 Dongting Dadao, Changde City, Hunan Province 415000, China
22institutetext: Hunan Province Key Laboratory Integration and Optical
Manufacturing Technology, 3150 Dongting Dadao, Changde City, Hunan Province
415000, China 33institutetext: Center for Astrophysics, Guangzhou University,
230 West Ring Road, Guangzhou, Guangdong Province 510006, China

# Accretion of the degenerate Fermi gas onto a Reissner-Nordström black hole

Ping Li 1122 Jiang-he Yang 1133 Siwei Xu 1122

(Received: date / Revised version: date)

###### Abstract

We investigate the stationary, spherically symmetric accretion of a degenerate
relativistic Fermi gas onto a Reissner-Nordström black hole. The accretion
theory is based on the Boyer-Lindquist coordinates and the Fermi gas follows
Fermi-Dirac statistics at infinity. We have derived the expression for the
particle current density, the stress energy-momentum tensor, and three
accretion rates. As the charged particle falls into the black hole, both the
mass and the charge of the black hole increase. Consequently, the mass
accretion rate and charge accretion rate are proportional to the particle
accretion rate. We have also provided analytical results at infinity and
numerical results within a finite range for these quantities. Our results
indicate that the accretion rate decreases as the charge of the black hole
increases, suggests that the presence of naked singularities is avoided in
black hole accretion theory. In this paper, we also discovered that the
accretion model of Vlasov gas behaves as an anisotropic fluid containing two
perfect-fluid components. One component represents the isotropic fluid of
Fermi gas, while the other represents a null fluid. When using the Boyer-
Lindquist coordinate system, we observed that the contribution from the null
fluid persists even at infinity, which led to the radial pressure always
smaller than the tangential pressure. Therefore, it’s not appropriate to treat
the accretion model as a perfect fluid at infinity.

###### pacs:

04.20.Jbdiscribing text of that key and 04.50.Kddiscribing text of that key
and 04.70.Bwdiscribing text of that key and 02.30.Jrdiscribing text of that
key

††offprints:

##  1 Introduction

In 2017, Rioseco and Sarbach proposed an accretion model Rioseco2017-1  for
Schwarzschild black holes that gained considerable attention. They assumed
that the gas, composed of a large number of particles, is in local
thermodynamic equilibrium at infinity, with rare collisions between particles
and negligible self-gravity, also known as the collisionless Vlasov gas. The
conditions of rare collisions was supported by observations from the Event
Horizon Telescope Collaboration EHTL1 ; EHTL5 ; EHTL12 ; EHTL16 . As the
particles travel along the geodesic into finite range, they form steady-state,
spherically symmetric accretion flows. The researchers solved the relativistic
Vlasov equation in Schwarzschild geometry and obtained the condition satisfied
by the distribution function with an appropriate symmetry. In the case of the
Maxwell-Jüttner distribution, they derived the current density, the stress
energy-momentum tensor, and the mass accretion rate. In another paper
Rioseco2017-2 , they also performed numerical computations of these quantities
and found that the radial pressure and tangential pressure approximately equal
each other at infinity, suggesting the accretion model can be treated as an
isotropic perfect fluid. However, as an observer approaches the horizon, the
tangential pressure increases steadily while the radial pressure decreases. At
the horizon, the tangential pressure becomes larger than the radial pressure.
This provides a partial explanation for the low accretion rate problem, also
known as the Bondi-Hoyle-Lyttleton accretion problem Hoyle1939 ; Lyttleton1940
; Bondi1944 ; Chaverra2015 .

Since then, many studies have expanded this model to encompass more general
scenarios. The Rioseco and Sarbach’s accretion model for the Reissner-
Nordström black hole were generalized by Ciśslik and Mach Cieslik2020 . While
in Refs. Mach2021-1 ; Mach2021-2 ; Mach2022 , authors delved into the
accretion of relativistic Vlasov gas onto a moving Schwarzschild black hole.
The results are using to model the growth of primordial black holes in the
early Universe. Gamboa et. al. also explored the accretion of Vlasov gas onto
a Schwarzschild black hole at a sphere of finite radius Gamboa2021 .
Additionally, beyond static spherical black holes in general relativity,
Rioseco and Sarbach’s results were extended to black holes in modified
gravities. Liao and Liu looked into the accretion of collisionless Vlasov gas
onto a Bardeen regular black hole Liao2022 , while Cai and Yang studied the
accretion of Vlasov gas onto black holes in bumblebee gravity Cai2022 . There
are also numerous studies examining accretion onto Kerr black holes. Ref.
Rioseco2018  illustrated that a collisionless, relativistic kinetic gas
configuration propagating in the equatorial plane of a Kerr black hole would
ultimately settle down to a stationary, axisymmetric configuration surrounding
the black hole. In the slowly rotating case, Ref. Andersson2018  proved the
decay of bounded energy and integrated energy for massless Vlasov gas in Kerr
spacetime. Subsequently, Cieślik et. al. conducted a comprehensive study of
the accretion of Vlasov gas onto a Kerr black hole in the equatorial plane
Cieslik2022 . Ref. Mach2023  discussed the equatorial accretion on a moving
Kerr black hole, while recently, a full $3+1$-dimensional accretion theory of
Vlaso gas onto a Kerr black hole was established Li2023 . Moreover, some
studies have started using numerical simulations to model Vlasov gas accretion
processes. Cieślik, Mach, and Odrzywołek performed Monte Carlo simulations of
stationary planar accretion of Vlasov gas onto a Schwarzschild black hole in
Refs. Cieslik2023 ; Cieslik2024 .

The process of a black hole accreting a red giant is a significant event in
astronomical observations Pacucci2022 ; Kareem2023 . Within red giants,
degenerate electrons are abundant and play a crucial role in their properties
Deheuvels2022 . When a black hole at the center accretes the red giant, the
degenerate electrons would be pulled into the black hole’s event horizon.
Taking quantum effects into account, electrons are categorized as Fermions and
adhere to the Fermi-Dirac statistics. This paper delves into the investigation
of the accretion of Fermi gas onto black holes. Although true astrophysical
plasmas are nearly electrically neutral, the theoretical importance of
studying the accretion of Fermions is evident. According to general
relativity, all black holes are classified as Kerr-Newman black holes, defined
by their mass $M$, charge $Q$, and angular momentum $a$. In astronomy, the
focus is also on the three properties of black holes. In order to alter the
charge of a black hole, the particle falling into the black hole must carry a
charge, too. As charged particles fall into the black hole, both the black
hole’s mass and charge change, as does the charge-to-mass ratio of the black
hole. However, as shown in Ref. Wald1974 , the particles that would violate
the horizon will not be captured by the black hole.

This paper is organized as follows: In Sec.II, we review the trajectory of
charged particles in the Reissner-Nordström black hole. In Sec.III, we discuss
the general Vlasov equation in the Boyer-Lindquist coordinates of the
Reissner-Nordström background. In Sec.IV, we choose the Fermi-Dirac statistic
as the distribution function at infinity. In Sec.V, we derive the expressions
of the particle current density and the stress energy momentum tensor, both
for the absorption part and the scattering part. In Sec.VI, we compare these
quantities with those of an isotropic perfect fluid and found that our model
describes an anisotropic fluid with two perfect-fluid components. In Sec.VII,
we present three accretion rates: the mass accretion rate, the charge
accretion rate, and the energy accretion rate. In Sec.VIII, we calculate the
asymptotic behavior of observable quantities using Taylor series. In Sec.IX,
we numerically calculate the corresponding quantities in a finite range. The
last section is a detailed conclusion of this paper. In this paper, we use
units $G=c=1$.

##  2 The motion of charged particle in Reissner-Nordström black hole: Review

In this paper, we consider the background to be Reissner-Nordström black
holes, which have the form

$\displaystyle ds^{2}$
$\displaystyle=-\left(1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}\right)dt^{2}+\left(1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}\right)^{-1}dr^{2}$
$\displaystyle+r^{2}d\theta^{2}+r^{2}\sin^{2}\theta d\varphi^{2},$ (1)

where $M$ is the mass and $Q$ is the charge of black hole. The outer horizon
is $r_{+}=M+\sqrt{M^{2}-Q^{2}}$. The electromagnetic potential has components
$A_{\mu}=(-\frac{Q}{r},0,0,0)$. The charged test particle should no longer
describe a geodesic in the Reissner-Nordström metric Chandrasekhar1983 . The
motion of charged particles traveling in the Reissner-Nordström geometry is
determined by the Lagrangian

$\displaystyle\mathcal{L}$
$\displaystyle=\frac{1}{2}g_{\mu\nu}\frac{dx^{\mu}}{d\lambda}\frac{dx^{\nu}}{d\lambda}+\frac{q}{m}A_{\mu}\frac{dx^{\mu}}{d\lambda}$
$\displaystyle=-\frac{1}{2}\left(1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}\right)\dot{t}^{2}+\frac{1}{2}\left(1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}\right)^{-1}\dot{r}^{2}$
$\displaystyle+\frac{1}{2}r^{2}\dot{\theta}^{2}+\frac{1}{2}r^{2}\sin^{2}\theta\dot{\varphi}^{2}-q\frac{Q}{r}\dot{t},$
(2)

where $\lambda$ is the affine parameter and the dot denotes the differential
with respect to $\lambda$. $m$ is the mass and $q$ is the charge of particle.

The canonical momentums $p_{\mu}$ are defined by

$\displaystyle p_{t}$
$\displaystyle=-\frac{\partial\mathcal{L}}{\partial\dot{t}}=\left(1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}\right)\dot{t}+q\frac{Q}{r},$
(3) $\displaystyle p_{r}$
$\displaystyle=\frac{\partial\mathcal{L}}{\partial\dot{r}}=\left(1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}\right)^{-1}\dot{r},$
(4) $\displaystyle p_{\theta}$
$\displaystyle=\frac{\partial\mathcal{L}}{\partial\dot{\theta}}=r^{2}\dot{\theta},$
(5) $\displaystyle p_{\varphi}$
$\displaystyle=\frac{\partial\mathcal{L}}{\partial\dot{\varphi}}=r^{2}\sin^{2}\theta\dot{\varphi}.$
(6)

Through the Legendre transformation, one can obtain the Hamiltonian

$\displaystyle\mathcal{H}$ $\displaystyle=p_{\mu}\dot{x}^{\mu}-\mathcal{L}$
$\displaystyle=\frac{1}{2}g_{\mu\nu}\frac{dx^{\mu}}{d\lambda}\frac{dx^{\nu}}{d\lambda}.$
(7)

The equations of motion are given by

$\frac{dx^{\mu}}{d\lambda}=\frac{\partial\mathcal{H}}{\partial
p_{\mu}},\quad\frac{dp_{\mu}}{d\lambda}=-\frac{\partial\mathcal{H}}{\partial
x^{\mu}}.$ (8)

Since the Hamiltonian does not contain $t$ and $\varphi$, the energy
$E=-p_{t}$ and the angular momentum $l_{z}=p_{\varphi}$ in the $z$-direction
are constants, which are

$\dot{p}_{t}=-\frac{\partial\mathcal{H}}{\partial
t}=0,\quad\dot{p}_{\varphi}=\frac{\partial\mathcal{H}}{\partial\varphi}=0.$
(9)

The Hamiltonian is conserved, defining another constant: the rest mass $m$

$m^{2}=-2\mathcal{H}.$ (10)

The last constant in the spherical symmetric spacetime is the angular momentum
$l^{2}$. To sum up, all four constants in this motion are

$\displaystyle m^{2}=$
$\displaystyle-2\mathcal{H}=-g_{\mu\nu}\frac{dx^{\mu}}{d\lambda}\frac{dx^{\nu}}{d\lambda},$
(11) $\displaystyle E=$
$\displaystyle-\left(1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}\right)\dot{t}-q\frac{Q}{r},$
(12) $\displaystyle l^{2}=$ $\displaystyle
r^{4}\left((\frac{d\theta}{d\lambda})^{2}+\sin^{2}\theta(\frac{d\varphi}{d\lambda})^{2}\right),$
(13) $\displaystyle l_{z}=$ $\displaystyle
r^{2}\sin^{2}\theta\frac{d\varphi}{d\lambda}.$ (14)

Using constants, the equations of motion (8) can be written as the decoupled,
first-order equations

$\displaystyle\dot{t}$
$\displaystyle=-\left(1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}\right)^{-1}\left(E+q\frac{Q}{r}\right),$
(15) $\displaystyle\dot{\varphi}$
$\displaystyle=\frac{l_{z}}{r^{2}\sin^{2}\theta},$ (16)
$\displaystyle\dot{r}^{2}$
$\displaystyle=-\left(1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}\right)(m^{2}+\frac{l^{2}}{r^{2}})+(E+q\frac{Q}{r})^{2},$
(17) $\displaystyle\dot{\theta}^{2}$
$\displaystyle=\frac{1}{r^{4}}\left(l^{2}-\sin^{-2}\theta l_{z}^{2}\right).$
(18)

For convenience to write, we define

$\displaystyle R$ $\displaystyle\equiv
r^{4}\dot{r}^{2}=(E^{2}-m^{2})r^{4}+2(m^{2}M+EqQ)r^{3}$
$\displaystyle+((q^{2}-m^{2})Q^{2}-l^{2})r^{2}+2l^{2}Mr-l^{2}Q^{2},$ (19)

which can be understood as the equivalent kinetic energy in the radial
direction.

When the charged particles incident from infinity towards the Reissner-
Nordström black hole, the motions are bounded into a plane since the angular
momentum $l^{2}$ is conserved. Similar to the geodesic motion in the
Schwarzschild spacetime, we only consider the plane motion in
$\theta=\frac{\pi}{2}$. At infinity, the equivalent kinetic energy $R>0$
indicates $E>m$. The equivalent kinetic energy $R$ also varied with respect to
$r$. Since the kinetic energy can not be negative, particles only appeared in
the range $R\geq 0$. Let’s define $r_{0}$ as the biggest real root of equation
$R(r_{0})=0$. A test particle can travel to $r=r_{0}$. There exist three cases
to discuss, also see Li2023 .

  1. 1.

$r_{0}<r_{+}$: The test particle falls into the black hole.

  2. 2.

$r_{0}>r_{+}$, and $R<0$ for $r<r_{0}$: The test particle cannot enter the
range $r<r_{0}$. It will be scattered by the black hole and travel to
infinity. In this case, the position $r=r_{0}$ is the closest location, also
known as the perihelion $r_{p}=r_{0}$, when particle traveling through the
black hole.

  3. 3.

$r_{0}>r_{+}$, and $R>0$ for $r<r_{0}$: The test particle cannot travel to
$r<r_{0}$ but keep getting close to $r=r_{0}$ and never reach there. The
position $r=r_{0}$ stands for an unstable orbit $r_{c}=r_{0}$.

Case 3 is a critical situation lying between case 1 and case 2. Since $R(r)>0$
is true for both $r>r_{0}$ and $r<r_{0}$, $R(r_{0})=0$ is the minimum of
function $R(r)$. Then, case 3 satisfies

$\displaystyle R(r_{c},l_{c})$ $\displaystyle=0,$ (20)
$\displaystyle\partial_{r}R(r_{c},l_{c})$ $\displaystyle=0.$ (21)

The above equations (20) and (21) are the key equations to classify the
trajectories of charged particles.

##  3 Action-angle variables and the Vlasov equation

Let $\gamma$ be an orbit of charged particles traveling in the Reissner-
Nordström background. The abbreviated action is expressed as

$S=-Et+l_{z}\varphi+\int_{\gamma}p_{r}dr+\int_{\gamma}p_{\theta}d\theta,$ (22)

where the integral is along the orbit $\gamma$. Since $E$ contains $Q$, the
electromagnetic interaction is already included in the abbreviated action.

Introducing the new momenta

$\displaystyle P_{0}$ $\displaystyle=m=\sqrt{-2\mathcal{H}},$ (23)
$\displaystyle P_{1}$ $\displaystyle=E=-p_{t},$ (24) $\displaystyle P_{2}$
$\displaystyle=l_{z}=p_{\varphi},$ (25) $\displaystyle P_{3}$
$\displaystyle=l=\sqrt{p_{\theta}^{2}+\frac{p_{\varphi}^{2}}{\sin^{2}\theta}},$
(26)

we can transform the canonical coordinates $(x^{\mu},p_{\mu})$ to the new
canonical coordinates $(Q^{\mu},P_{\mu})$. The new conjugate coordinates are
defined as

$\displaystyle Q^{0}$ $\displaystyle=\frac{\partial S}{\partial
m}=-\int_{\gamma}\frac{r^{2}}{\sqrt{R}}\left(1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}\right)^{-1}$
$\displaystyle\cdot\left(m\left(1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}\right)-\frac{Q}{r}\frac{q}{m}\left(E+\frac{qQ}{r}\right)\right)dr,$
(27) $\displaystyle Q^{1}$ $\displaystyle=\frac{\partial S}{\partial
E}=-t+\int_{\gamma}\frac{r^{2}}{\sqrt{R}}\left(1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}\right)^{-1}$
$\displaystyle\cdot\left(E+\frac{qQ}{r}\right)dr,$ (28) $\displaystyle Q^{2}$
$\displaystyle=\frac{\partial S}{\partial l_{z}}=\varphi-
l_{z}\int_{\gamma}\frac{\sin^{-2}\theta}{\sqrt{l^{2}-\sin^{-2}\theta
l_{z}^{2}}}d\theta,$ (29) $\displaystyle Q^{3}$ $\displaystyle=\frac{\partial
S}{\partial
l}=-l\int_{\gamma}\frac{1}{\sqrt{R}}dr+l\int_{\gamma}\frac{d\theta}{\sqrt{l^{2}-\sin^{-2}\theta
l_{z}^{2}}}.$ (30)

The new canonical coordinates $(Q^{\mu},P_{\mu})$ are also known as the
action-angle variables.

Assume that there are a large number of charged particles, and each of them
travels along the orbit described by the abbreviated action (22). At the same
time, their behavior also can be treated as an ensemble, which is described by
the distribution function $f(x^{\mu},p_{\nu})$. Furthermore, we consider that
each particle interacts with others only through the gravitational and
electromagnetic fields of the background. Thus, the distribution function
$f(x^{\mu},p_{\nu})$ satisfied the Vlasov equation

$\dot{f}\equiv\\{f,\mathcal{H}\\}=0,$ (31)

where the bracket $\\{.,.\\}$ stands for the Poisson bracket.

The transformations $(x^{\mu},p_{\mu})\rightarrow(Q^{\mu},P_{\mu})$ keep the
Poisson bracket (31) covariant

$\frac{\partial H}{\partial p_{\mu}}\frac{\partial}{\partial
x^{\mu}}-\frac{\partial H}{\partial x^{\mu}}\frac{\partial}{\partial
p_{\mu}}=\frac{\partial H}{\partial P_{\mu}}\frac{\partial}{\partial
Q^{\mu}}-\frac{\partial H}{\partial Q^{\mu}}\frac{\partial}{\partial
P_{\mu}}.$ (32)

In the action-angle variables $(Q^{\mu},P_{\mu})$, the Hamiltonian becomes
$H=\frac{P_{0}^{2}}{2}$. Therefore, the Vlasov equation (31) is simplified to

$\frac{\partial}{\partial Q^{0}}f=0,$ (33)

which is to say the distribution function $f$ does not involve $Q^{0}$. Using
the symmetries, one can constrain the distribution function
$f(x^{\mu},p_{\mu})$ to be a more simple form. This is first done by Ref.
Rioseco2017-1 , which shows for the stationary spherical symmetric flow, the
distribution function $f$ is independent of $Q^{1},Q^{2},Q^{3}$ and $P_{2}$.
That is to say, the solution of the Vlasov equation in the stationary
spherical symmetric case has a general form

$f(x^{\mu},p_{\nu})=\mathcal{F}(P_{0},P_{1},P_{3}).$ (34)

##  4 degenerate relativistic Fermi gas

Astronomical observations often show black holes gradually pulling in massive
objects, many of which are made of Fermions. To paint a more realistic picture
of how black holes accrete matter, we examine a degenerate fluid composed of a
large number of Fermions. At infinity, the Reissner-Nordström metric is
asymptotically flat and we assume that the Fermi gas is in thermal
equilibrium. Due to the black hole’s gravitational pull, the Fermi gas should
be limited to a finite range along its trajectory. Fermions adhere to Pauli’s
exclusion principle, meaning that each quantum state of an electron is allowed
to be occupied by only one particle. The statistical law that describes their
behavior at infinity is known as Fermi-Dirac statistics, which is given by the
expression Rezzolla2013

$f_{\infty}=A\frac{1}{\exp[(E-\mu)/k_{B}T]+1},$ (35)

where $A$ is the normalization factor. $k_{B}$ is the Boltzmann constant, $T$
is the temperature and $\mu$ is the relativistic chemical potential. The above
distribution function (35) is often used to describe electrically neutral
Fermi systems. If particles have a charge $q$, electromagnetic interactions
between point particles may disrupt the distribution. However, we emphasize
that the Vlasov equation (31) is only valid when ignoring the electromagnetic
interactions between particles. If the particle’s electromagnetic interaction
is taken into account, the Poisson bracket (31) would no longer be zero, and
the distribution function (35) would change. Therefore, the selection of
distribution function (35) needs to consider the conditions of Vlasov
equation. The above distribution function is valid in at least two cases. The
first case is an electrically neutral Fermion gas, such as a neutron gas. The
second case is the zero-order effect of the electron gas. Since the
interaction term is related to the charge $q$, it is sufficient to consider
the zero-order effects when $q$ is small.

Similar to Ref. Rioseco2017-1 , the distribution function throughout all
spacetime is giving by

$f=\delta(P_{0}-m_{0})f_{\infty}=A\delta(P_{0}-m_{0})\frac{1}{\exp[(P_{1}-\mu)/k_{B}T]+1},$
(36)

where $m_{0}$ is the mass of electron. The above distribution function $f$ is
a specific solution of Vlasov equation and obeys Eq.(34). It is valid
everywhere, not only at infinity. However, in a finite radius, $T$ is only a
parameter and does not stand for temperature.

##  5 Observable quantities

The observable quantities we care about are the particle current density
$J_{\mu}$ and the stress energy-momentum tensor $T_{\mu\nu}$. They are given
by a suitable integral over momenta

$\displaystyle J_{\mu}$ $\displaystyle=\int
p_{\mu}f(x,p)\mathrm{dvol}_{x}(p),$ (37) $\displaystyle T_{\mu\nu}$
$\displaystyle=\int p_{\mu}p_{\nu}f(x,p)\mathrm{dvol}_{x}(p),$ (38)

where the volume element is given by

$\mathrm{dvol}_{x}(p)=\sqrt{-\det[g^{\mu\nu}]}dp_{0}dp_{1}dp_{2}dp_{3}.$ (39)

Using Eq. (31), one can easily show that the particle current density
$J_{\mu}$ and the energy-momentum tensor $T_{\mu\nu}$ satisfy the conservation
laws

$\displaystyle\nabla_{\mu}J^{\mu}$ $\displaystyle=0,$ (40)
$\displaystyle\nabla_{\mu}T^{\mu\nu}$ $\displaystyle=0.$ (41)

Notice that the observable quantities are actually calculated in the action-
angle variables $(P_{0},P_{1},P_{2},P_{3})$. We then transform the momentum
variables $(p_{t},p_{r},p_{\theta},p_{\varphi})$ to the variables
$(m^{2},E,l_{z},l^{2})$, which are expressed as follow

$\displaystyle m^{2}$
$\displaystyle=-g^{00}(p_{t}-\frac{q}{m}\frac{Q}{r})^{2}-g^{11}p_{r}^{2}-g^{22}p_{\theta}^{2}-g^{33}p_{\varphi}^{2},$
(42) $\displaystyle E$ $\displaystyle=-p_{t},$ (43) $\displaystyle l_{z}$
$\displaystyle=p_{\varphi},$ (44) $\displaystyle l^{2}$
$\displaystyle=p_{\theta}^{2}+\sin^{-2}\theta p_{\varphi}^{2}.$ (45)

The Jacobian determinant is defined by

$J=\frac{\partial(m^{2},E,l_{z},l^{2})}{\partial(p_{t},p_{r},p_{\theta},p_{\varphi})}=\frac{4\sqrt{R}}{r^{2}}\sqrt{l^{2}-\sin^{-2}\theta
l_{z}^{2}}.$ (46)

Therefore, the volume element can be reexpressed as

$\mathrm{dvol}_{x}(p)=\frac{1}{4\sqrt{R}}\frac{1}{\sqrt{l^{2}\sin^{2}\theta-
l_{z}^{2}}}dm^{2}dEdl_{z}dl^{2}.$ (47)

In order to keep the function inside the square root non-negative, the
integral range of $l_{z}$ is $l_{z}\in[-l\sin\theta,l\sin\theta]$. We further
define

$l_{z}=l\sin\theta\sin\sigma,$ (48)

and transform the integral of $l_{z}$ to the integral of $\sigma$. It is easy
to obtain

$\mathrm{dvol}_{x}(p)=\frac{1}{4\sqrt{R}}dm^{2}dEd\sigma dl^{2},$ (49)

and the integral range of $\sigma$ is $\sigma\in[0,2\pi]$.

According to the classification of trajectories, the observable quantities
$J_{\mu}$ and $T_{\mu\nu}$ can be divided into three parts: the absorbed part,
the scattered part and the critical part, which are

$\displaystyle J_{\mu}$
$\displaystyle=J^{abs}_{\mu}+J^{scat}_{\mu}+J^{cri}_{\mu},$ (50)
$\displaystyle T_{\mu\nu}$
$\displaystyle=T^{abs}_{\mu\nu}+T^{scat}_{\mu\nu}+T^{cri}_{\mu\nu}.$ (51)

The critical part does not have a substantial contribution and attracts no
attention in most cases. The integral interval of the parameter
$(m,E,\sigma,l)$ will determine which part the quantities $J_{\mu}$ and
$T_{\mu\nu}$ belong to.

###  5.1 Absorption part

For the absorption part, there is no perihelion on the orbit. Solving the
simultaneous Eq. (20) and (21), one obtain

$r=r_{c}(E),\quad l=l_{c}(E),$ (52)

where $r_{c}(E)$ is the unstable orbit. Notice that the expression (19) of $R$
does not contained $l_{z}$, the solutions $(r_{c}(E),l_{c}(E))$ are
independent of $\sigma$. $l_{c}(E)$ determined whether there is a perihelion.
For trajectories with parameter $l>l_{c}$, the perihelion exists and the
charged particles should be scattered to infinity; otherwise, the charged
particles would be captured by the black hole. Thus, for the absorption part,
the integral interval is $l\in[0,l_{c})$.

Plugging the distribution function (36) into the observable quantities (37),
we obtain the absorption part of $J_{\mu}$ as follows

$\displaystyle J^{abs}_{t}$ $\displaystyle=-2\pi
Am_{0}\int_{m_{0}}^{\infty}dE\int_{0}^{l_{c}}dl\frac{El}{\sqrt{R_{0}}}F,$ (53)
$\displaystyle J^{abs}_{r}$ $\displaystyle=2\pi
Am_{0}\int_{m_{0}}^{\infty}dE\int_{0}^{l_{c}}dl\frac{l}{(r^{2}-2Mr+Q^{2})}F,$
(54)

where $F=\frac{1}{(\exp[(E-\mu)\beta]+1)},\beta=\frac{1}{k_{b}T}$ and $R_{0}$
stands for the function $R$ in the case of $m=m_{0}$. Since
$p_{\theta}=l\sin\sigma$ and $p_{\varphi}=l\sin\theta\sin\sigma$, integrating
$\sigma$ in a period, one can easily show $J^{abs}_{\theta}$ and
$J^{abs}_{\varphi}$ vanish.

For the same reason, one can also show that
$T_{t\theta}=T_{t\varphi}=T_{r\theta}=T_{r\varphi}=T_{\theta\varphi}=0$. The
non-zero components of $T_{\mu\nu}$ are expressed as

$\displaystyle T^{abs}_{tt}$ $\displaystyle=2\pi
Am_{0}\int_{m_{0}}^{\infty}dE\int_{0}^{l_{c}}dl\frac{E^{2}l}{\sqrt{R_{0}}}F,$
(55) $\displaystyle T^{abs}_{tr}$ $\displaystyle=-2\pi
Am_{0}\int_{m_{0}}^{\infty}dE\int_{0}^{l_{c}}dl\frac{El}{(r^{2}-2Mr+Q^{2})}F,$
(56) $\displaystyle T^{abs}_{rr}$ $\displaystyle=2\pi
Am_{0}\int_{m_{0}}^{\infty}dE\int_{0}^{l_{c}}dl\frac{\sqrt{R_{0}}l}{(r^{2}-2Mr+Q^{2})^{2}}F,$
(57) $\displaystyle T^{abs}_{\theta\theta}$
$\displaystyle=\frac{T^{abs}_{\varphi\varphi}}{\sin^{2}\theta}=\pi
Am_{0}\int_{m_{0}}^{\infty}dE\int_{0}^{l_{c}}dl\frac{l^{3}}{\sqrt{R_{0}}}F.$
(58)

###  5.2 Scattering part

For the scattering part, the condition for the existence of perihelion
requirements $l>l_{c}$. At the same time, the signature of $R_{0}$ must be
positive. Solving $R_{0}=0$, one obtains the maximum of $l_{max}(E)$, which is

$l_{max}=r\sqrt{\frac{r^{2}}{r^{2}-2Mr+Q^{2}}\left(E+\frac{qQ}{r}\right)^{2}-m_{0}^{2}}.$
(59)

Thus, the integral interval of $l$ is $l\in(l_{c},l_{max}]$.

The integral of $E$ is more complicated. The main reason is that the closer a
particle gets to the black hole, the more energy it takes to escape. We
rewrite $R_{0}$ into

$R_{0}=r^{4}\left((E+q\frac{Q}{r})^{2}-\left(1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}\right)(m_{0}^{2}+\frac{l^{2}}{r^{2}})\right),$
(60)

and define the potential

$U_{l}(r)=\left(1-\frac{2M}{r}+\frac{Q^{2}}{r^{2}}\right)(m_{0}^{2}+\frac{l^{2}}{r^{2}}).$
(61)

Since $R_{0}\geq 0$, the energy $E$ required for a particle to escape the
black hole satisfy $E\geq\sqrt{U_{l}}-\frac{qQ}{r}$. We define

$\displaystyle\Psi_{l}(r)$ $\displaystyle=\sqrt{U_{l}}-\frac{qQ}{r}.$ (62)

In the maximum of $\Psi_{l}(r)$, i.e. $\frac{d\Psi_{l}(r)}{dr}=0$, there is

$\displaystyle l^{2}_{\pm}$
$\displaystyle=\frac{1}{2Y^{2}}\bigg{(}2Y(\Delta-Y)m_{0}^{2}r^{2}+q^{2}Q^{2}r^{2}\Delta$
$\displaystyle\pm qQr^{2}\Delta\sqrt{4Ym_{0}^{2}+q^{2}Q^{2}}\bigg{)},$ (63)

where

$\displaystyle\Delta(r)$ $\displaystyle=r^{2}-2Mr+Q^{2}=(r-r_{+})(r-r_{-}),$
(64) $\displaystyle Y(r)$
$\displaystyle=r^{2}-3Mr+2Q^{2}=(r-r_{ph+})(r-r_{ph-}).$ (65)

The expressions $r_{\pm}=M\pm\sqrt{M^{2}-Q^{2}}$ are the outer or inner
horizon and $r_{ph\pm}=\frac{1}{2}(3M\pm\sqrt{9M^{2}-4Q^{2}})$ are the radii
of photo orbit. In the special case of $q=0$, we reproduce the result in ref.
Cieslik2020

$l^{2}_{+}=l^{2}_{-}=\frac{(\Delta-Y)}{Y}m_{0}^{2}r^{2}.$ (66)

Inserting $l^{2}_{+}$ into $\Psi_{l}$, one have

$\Psi_{l_{+}}=\frac{m_{0}\Delta}{r\sqrt{Y}}\sqrt{1+\frac{1}{2Y}\frac{qQ}{m_{0}}\left(\frac{qQ}{m_{0}}-\sqrt{4Y+\frac{q^{2}Q^{2}}{m_{0}^{2}}}\right)}-\frac{qQ}{r}.$
(67)

As $r$ approaching to outer phone sphere $r\rightarrow r_{ph+}$, the function
$\Psi_{l_{+}}\rightarrow\infty$. That is to say, there is no particle that can
escape the black hole. On the other hand, for particles that can incident from
infinity, their energy satisfies $E\geq m_{0}$. Or equivalent to say, the
minimal energy always satisfies $E_{min}\geq m_{0}$. Therefore, the radial
coordinate is divided into three segments

$E_{min}=\left\\{\begin{array}[]{ll}\infty,&\hbox{$r<r_{ph+}$;}\\\
\Psi_{l_{+}},&\hbox{$r_{1}>r>r_{ph+}$;}\\\
m_{0},&\hbox{$r>r_{1}$.}\end{array}\right.$ (68)

where $r_{1}$ is the root of $\Psi_{l_{+}}(r_{1})=m_{0}$. At last, we obtain
that the integral interval of $E$ is $E\in[E_{min},\infty)$.

Notice that the trajectory of the scattering part contains both
$+\sqrt{R_{0}}$ and $-\sqrt{R_{0}}$. Therefore, we expressed the non-zero
observable quantities as follows

$\displaystyle J^{scat}_{t}$ $\displaystyle=-2\pi
Am_{0}\sum_{\pm}\int_{E_{min}}^{\infty}dE\int_{l_{c}}^{l_{max}}dl\frac{El}{\sqrt{R_{0}}}F,$
(69) $\displaystyle J^{scat}_{r}$ $\displaystyle=2\pi
Am_{0}\sum_{\pm}\int_{E_{min}}^{\infty}dE\int_{l_{c}}^{l_{max}}dl\frac{l}{(r^{2}-2Mr+Q^{2})}F,$
(70) $\displaystyle T^{scat}_{tt}$ $\displaystyle=2\pi
Am_{0}\sum_{\pm}\int_{E_{min}}^{\infty}dE\int_{l_{c}}^{l_{max}}dl\frac{E^{2}l}{\sqrt{R_{0}}}F,$
(71) $\displaystyle T^{scat}_{tr}$ $\displaystyle=-2\pi
Am_{0}\sum_{\pm}\int_{E_{min}}^{\infty}dE\int_{l_{c}}^{l_{max}}dl\frac{El}{(r^{2}-2Mr+Q^{2})}F,$
(72) $\displaystyle T^{scat}_{rr}$ $\displaystyle=2\pi
Am_{0}\sum_{\pm}\int_{E_{min}}^{\infty}dE\int_{l_{c}}^{l_{max}}dl\frac{\sqrt{R_{0}}l}{(r^{2}-2Mr+Q^{2})^{2}}F,$
(73) $\displaystyle T^{scat}_{\theta\theta}$
$\displaystyle=\frac{T^{scat}_{\varphi\varphi}}{\sin^{2}\theta}=\pi
Am_{0}\sum_{\pm}\int_{E_{min}}^{\infty}dE\int_{l_{c}}^{l_{max}}dl\frac{l^{3}}{\sqrt{R_{0}}}F.$
(74)

##  6 Comparison with the perfect fluid case

In this section, we make a comparison between the observable quantities of
Vlasov gas with those of an isotropic perfect fluid. In the perfect fluid
model, the particle current density and the stress energy-momentum tensor are
given by

$J^{\mu}=nu^{\mu},\quad T^{\mu}{}_{\nu}=\rho
u^{\mu}u_{\nu}+p(u^{\mu}u_{\nu}+\delta^{\mu}_{\nu}),$ (75)

where $n$ is particle density, $u^{\mu}$ is four-velocity, $\rho$ is the
energy density and $p$ is the pressure. As is shown in Ref. Cieslik2020 ,
$u^{\mu}$ is an eigenvector of $T^{\mu}{}_{\nu}$ corresponding to the
eigenvalue $-\rho$. For any non-zero $k^{\mu}$ orthogonal to $u_{\mu}$ is also
an eigenvector corresponding to the eigenvalue $p$. In the perfect fluid
model, the pressure $p$ is a threefold degenerate eigenvalue of
$T^{\mu}{}_{\nu}$.

In the Vlasov case, things would be very different. Since
$T_{\theta\varphi}=0$ and $T_{\varphi\varphi}=\sin^{2}\theta
T_{\theta\theta}$, there existed twofold eigenvalue
$p_{\mathrm{tan}}=T^{\theta}{}_{\theta}=T^{\varphi}{}_{\varphi}$ for
$T^{\mu}{}_{\nu}$, corresponding to the eigenvector $(0,0,1,0)$ and
$(0,0,0,1)$ respectively. In order to calculate the other two eigenvectors, we
consider the matrix $\left(\begin{array}[]{cc}T^{t}{}_{t}&T^{t}{}_{r}\\\
T^{r}{}_{t}&T^{r}{}_{r}\\\ \end{array}\right)$. Calculating the secular
equation, one obtains two eigenvalues

$\lambda_{\pm}=\frac{T^{t}{}_{t}+T^{r}{}_{r}\pm\sqrt{(T^{t}{}_{t}-T^{r}{}_{r})^{2}+4T^{t}{}_{r}T^{r}{}_{t}}}{2},$
(76)

which corresponding to the energy density $-\rho$ and the radial pressure
$p_{\mathrm{rad}}$. Thus, the general form of energy-momentum tensor is
expressed as

$T^{\mu}{}_{\nu}=\left(\begin{array}[]{cccc}-\rho&0&0&0\\\
0&p_{\mathrm{rad}}&0&0\\\ 0&0&p_{\mathrm{tan}}&0\\\ 0&0&0&p_{\mathrm{tan}}\\\
\end{array}\right).$ (77)

In the case of $p_{\mathrm{rad}}\neq p_{\mathrm{tan}}$, as is pointed in Ref.
Letelier1980 , this is an anisotropic fluid with two-perfect-fluid components.
One can construct the above energy-momentum tensor (77) by combining two
perfect fluids, one perfect and one null fluid, or two null fluids. In the
case of Vlasov gas, the energy-momentum tensor can be decoupled into one
perfect fluid corresponding to the Fermi gas and another null fluid. It is the
null fluid that distinguishes the Vlasov model from the perfect fluid model.
Latter, we will calculate the eigenvalues
$\rho,p_{\mathrm{rad}},p_{\mathrm{tan}}$ as the function of $r$.

##  7 Accretion rates

For the spherical symmetric flow, there are two accretion rates: one is the
particle accretion rate $\dot{n}$, meaning number of particles that cross the
horizon per unit time; the other is the energy accretion rate $\dot{E}$,
meaning total energy crossing the horizon per unit time. The particle number
is defined by

$n=\sqrt{-g_{\mu\nu}J^{\mu}J^{\nu}}.$ (78)

Then, the particle accretion rate $\dot{n}$ is also defined by the particle
current density $J^{\mu}$. Notice that $\nabla_{\mu}J^{\mu}=0$ becomes

$\frac{\partial}{\partial r}(r^{2}\sin\theta J^{r})=0.$ (79)

Let $V$ be the region of a black hole and $S$ is the surface of the outer
horizon. The vector $\hat{n}=(1,0,0)$ is the unit normal vector of $S$
directed outside. Integrating the above equation and using the Stokes’
theorem, the particle accretion rate $\dot{n}$ is defined by

$\displaystyle\dot{n}$ $\displaystyle=-\int_{V}(r^{2}\sin\theta J^{r})_{,r}dV$
$\displaystyle=-\int_{S}J^{r}(r_{+})\hat{n}_{r}r_{+}^{2}\sin\theta d\theta
d\varphi$ $\displaystyle=-4\pi r_{+}^{2}J^{r}(r_{+}).$ (80)

When an electron falls into the black hole, its mass $m_{0}$ and charge $q$
are absorbed by the black hole, so that the black hole’s mass increases by
$m_{0}$ and charge increases by $q$. Therefore, the growth of mass and charge
for the black hole is caused by the particle number which falls into the black
hole. We then define the mass accretion rates $\dot{M}$ and the charge
accretion rates $\dot{Q}$ as

$\displaystyle\dot{M}$ $\displaystyle=m_{0}\dot{n}=-4\pi
r_{+}^{2}m_{0}J^{r}(r_{+}),$ (81) $\displaystyle\dot{Q}$
$\displaystyle=q\dot{n}=-4\pi r_{+}^{2}qJ^{r}(r_{+}).$ (82)

Similar to $J^{\mu}$, the stress energy momentum tensor $T_{\mu\nu}$ are also
divergence free. Then, one can define the energy accretion rate

$\dot{E}=-4\pi r_{+}^{2}T^{r}{}_{t}(r_{+}).$ (83)

##  8 Asymptotic behavior

In this section, we discussed the asymptotic behavior of $J^{\mu}$ and
$T_{\mu\nu}$ by using the Taylor series. At large distances, the absorption
part becomes zero, and only contributions of the scattering part remain. We
rewrite $R_{0}$ to

$R_{0}(r)=X(r)+\Delta(r)l^{2},$ (84)

where

$\displaystyle X(r)$
$\displaystyle=(E^{2}-m_{0}^{2})r^{4}+2(Mm_{0}^{2}+EqQ)r^{3}+\left(q^{2}-m_{0}^{2}\right)Q^{2}r^{2}.$
(85)

The following integration of $l$ is obtained

$\int_{l_{c}}^{l_{max}}\frac{l}{\sqrt{R_{0}}}dl=\int_{l_{c}}^{l_{max}}\frac{l}{\sqrt{X-\Delta
l^{2}}}dl=\frac{\sqrt{X-\Delta l_{c}^{2}}}{\Delta},$ (86)

where we have used the identity $X+\Delta l_{max}^{2}=0$. Notice that
$l_{max}$ can be expand as $\frac{1}{r}$

$l_{max}^{2}=r^{2}\left(E^{2}-m_{0}^{2}+2(E^{2}M+EqQ)\frac{1}{r}+o(\frac{1}{r})\right),$
(87)

But $l_{c}$ does not depend on $r$. Then, one obtains

$\int_{l_{c}}^{l_{max}}\frac{l}{\sqrt{R_{0}}}dl=\sqrt{E^{2}-m_{0}^{2}}\left(1+\frac{2E^{2}M+EqQ-
m_{0}^{2}M}{(E^{2}-m_{0}^{2})}\frac{1}{r}\right).$ (88)

Therefore, using the Taylor expansion of $\frac{1}{r}$, we obtain

$\displaystyle J_{t}(r\rightarrow\infty)$ $\displaystyle=-4\pi
Am_{0}\int_{E_{min}}^{\infty}dE\int_{l_{c}}^{l_{max}}dl\frac{El}{\sqrt{R_{0}}}F$
$\displaystyle=-4\pi
Am_{0}\int_{m_{0}}^{\infty}\frac{\sqrt{E^{2}-m_{0}^{2}}E}{\exp[(E-\mu)\beta]+1}\bigg{(}1$
$\displaystyle+\frac{2E^{2}M+EqQ-
m_{0}^{2}M}{E^{2}-m_{0}^{2}}\frac{1}{r}\bigg{)}dE,$ (89)

where we only list the first two terms.

The function $\frac{1}{\sqrt{R_{0}}}$ can also be expanded in terms of
$\frac{1}{r}$

$\frac{1}{\sqrt{R_{0}}}=\frac{1}{\sqrt{E^{2}-m_{0}^{2}}}\left(\frac{1}{r^{2}}-\frac{qQE+m_{0}^{2}M}{E^{2}-m_{0}^{2}}\frac{1}{r^{3}}+o(\frac{1}{r^{3}})\right).$
(90)

Through a similar calculation, we obtain

$\displaystyle J_{r}(r\rightarrow\infty)$ $\displaystyle=2\pi
Am_{0}\int_{m_{0}}^{\infty}\frac{E^{2}-m_{0}^{2}}{\exp[(E-\mu)\beta]+1}\bigg{(}1$
$\displaystyle+\frac{2E^{2}M+EqQ-
m_{0}^{2}M}{E^{2}-m_{0}^{2}}\frac{2}{r}\bigg{)}dE,$ (91) $\displaystyle
T_{tt}(r\rightarrow\infty)$ $\displaystyle=4\pi
Am_{0}\int_{m_{0}}^{\infty}\frac{\sqrt{E^{2}-m_{0}^{2}}E^{2}}{\exp[(E-\mu)\beta]+1}\bigg{(}1$
$\displaystyle+\frac{2E^{2}M+EqQ-
m_{0}^{2}M}{E^{2}-m_{0}^{2}}\frac{1}{r}\bigg{)}dE,$ (92) $\displaystyle
T_{tr}(r\rightarrow\infty)$ $\displaystyle=-2\pi
Am_{0}\int_{m_{0}}^{\infty}\frac{(E^{2}-m_{0}^{2})E}{\exp[(E-\mu)\beta]+1}\bigg{(}1$
$\displaystyle+\frac{2E^{2}M+EqQ-
m_{0}^{2}M}{E^{2}-m_{0}^{2}}\frac{2}{r}\bigg{)}dE,$ (93) $\displaystyle
T_{rr}(r\rightarrow\infty)$ $\displaystyle=\frac{4}{3}\pi
Am_{0}\int_{m_{0}}^{\infty}\frac{(E^{2}-m_{0}^{2})^{\frac{3}{2}}}{\exp[(E-\mu)\beta]+1}\bigg{(}1$
$\displaystyle+\frac{2E^{2}M+EqQ-
m_{0}^{2}M}{E^{2}-m_{0}^{2}}\frac{3}{r}\bigg{)}dE,$ (94) $\displaystyle
T_{\theta\theta}(r\rightarrow\infty)$
$\displaystyle=\frac{T_{\varphi\varphi}(r\rightarrow\infty)}{\sin^{2}\theta}$
$\displaystyle=\frac{4r^{2}}{3}\pi
Am_{0}\int_{m_{0}}^{\infty}\frac{(E^{2}-m_{0}^{2})^{\frac{3}{2}}}{\exp[(E-\mu)\beta]+1}\bigg{(}1$
$\displaystyle+\frac{4E^{2}M+3EqQ-
m_{0}^{2}M}{E^{2}-m_{0}^{2}}\frac{1}{r}\bigg{)}dE.$ (95)

It is not difficult to find, the stress energy-momentum tensor $T_{\mu\nu}$ at
infinity satisfy

$T_{\infty}{}^{r}{}_{r}=T_{\infty}{}^{\theta}{}_{\theta}=T_{\infty}{}^{\varphi}{}_{\varphi}.$
(96)

For the classical particles which satisfied Maxwell-Jüttner distribution, the
above quantities turn into

$\displaystyle J_{\infty t}$ $\displaystyle=-4\pi
Am_{0}^{4}\frac{K_{2}(z)}{z},$ (97) $\displaystyle J_{\infty r}$
$\displaystyle=4\pi Am_{0}^{4}\frac{1+z}{z^{3}}e^{-z},$ (98) $\displaystyle
T_{\infty}{}^{t}{}_{t}$ $\displaystyle=4\pi
Am_{0}^{5}\left(\frac{K_{1}(z)}{z}+3\frac{K_{2}(z)}{z^{2}}\right),$ (99)
$\displaystyle T_{\infty}{}^{t}{}_{r}$
$\displaystyle=-T_{\infty}{}^{r}{}_{t}=4\pi
Am_{0}^{5}\frac{z^{2}+3z+3}{z^{4}}e^{-z},$ (100) $\displaystyle
T_{\infty}{}^{r}{}_{r}$
$\displaystyle=T_{\infty}{}^{\theta}{}_{\theta}=T_{\infty}{}^{\varphi}{}_{\varphi}=4\pi
Am_{0}^{5}\frac{K_{2}(z)}{z^{2}},$ (101)

where $K_{1}(x)$ is the first Bessel function, $K_{2}(x)$ is the secondary
Bessel function and $z=m_{0}\beta$. Since
$T_{\infty}{}^{t}{}_{r}=-T_{\infty}{}^{r}{}_{t}$ is not vanish, the Vlasov gas
is actually an anisotropic fluid with two perfect fluids even at infinity. The
energy-momentum tensor for the Fermi gas is an isotropic perfect fluid, which
has a form

$T^{(fer)}_{\infty}{}^{\mu}{}_{\nu}=4\pi
Am_{0}^{5}\left(\begin{array}[]{cccc}\frac{K_{1}(z)}{z}+3\frac{K_{2}(z)}{z^{2}}&0&0&0\\\
0&\frac{K_{2}(z)}{z^{2}}&0&0\\\ 0&0&\frac{K_{2}(z)}{z^{2}}&0\\\
0&0&0&\frac{K_{2}(z)}{z^{2}}\\\ \end{array}\right).$ (102)

Another component is a null fluid and has the energy-momentum tensor

$T^{(null)}_{\infty}{}_{\mu\nu}=4\pi
Am_{0}^{5}\frac{z^{2}+3z+3}{z^{4}}e^{-z}(\delta^{0}_{\mu}\delta^{1}_{\nu}+\delta^{1}_{\mu}\delta^{0}_{\nu}).$
(103)

When the temperature approaching to zero $z\rightarrow\infty$, i.e.
$T_{\infty}{}^{t}{}_{r}=-T_{\infty}{}^{r}{}_{t}\rightarrow 0$, the
collisionless gas behaves an ideal fluid at infinity. Then, the particle
number $n$ at large distance mainly contributed by $J_{t}$

$n_{\infty}=4\pi Am_{0}^{4}\frac{K_{2}(z)}{z}.$ (104)

For the same reason, the collisionless Fermi gas also can not be treated as
the ideal fluid at infinity. The radial pressure $p_{\mathrm{rad}}$ is not
equal to the tangential pressure $p_{\mathrm{tan}}$.

The above equations (89) - (95) also indicate that the physical results are
independent of the parameters $q$ and $Q$ at the leader term of Taylor’s
expansion. However, the distribution function affects significantly. When
performing the numerical results, this behavior is verified.

At the horizon, the scattering part is omitted and only the absorption part
contributed. Eq.(20) and (21) only permit numerical solution for $l_{c}$. In
the Schwarzschild limit $Q=0$, $l_{c}$ have the analytical expression, that is

$l_{c}^{2}=\frac{M^{2}}{2}\frac{27E^{4}-36E^{2}m_{0}^{2}+8m_{0}^{4}+E(9E^{2}-8m_{0}^{2})^{\frac{3}{2}}}{E^{2}-m_{0}^{2}}.$
(105)

Thus, for the Schwarzschild black hole to accrete the collisionless Fermi gas,
there is

$\frac{\dot{n}}{n_{\infty}}\bigg{|}_{Sch}=\pi\frac{\int_{m_{0}}^{\infty}\frac{l_{c}^{2}}{\exp[(E-\mu)\beta]+1}dE}{\int_{m_{0}}^{\infty}\frac{\sqrt{E^{2}-m_{0}^{2}}E}{\exp[(E-\mu)\beta]+1}dE}.$
(106)

Figure 1: The numerical results $\frac{\dot{n}}{n_{\infty}}\big{|}_{Sch}$ for
the accretion onto a Schwarzschild black hole, where the unit is been chosen
as $\frac{\pi M^{2}}{2}$.

We plot the numerical results of $\frac{\dot{n}}{n_{\infty}}\big{|}_{Sch}$
both for the Fermi gas and the classical gas in Fig.1. In the low temperature
limit $z\rightarrow\infty$, the number
$\frac{\dot{n}}{n_{\infty}}\big{|}_{Sch}$ equal to the same value.

##  9 Numerical results

###  9.1 particle current density and particle density

Figure 2: $\varepsilon_{min}$ as the function of $r$. We choose
$Q=0.9,\mu=0,z=1$.

The numerical calculation is performed by Wolfram Mathematica. We introduce
the following dimensionless variables

$\varepsilon=\frac{E}{m_{0}},\quad\xi=\frac{l}{m_{0}},\quad
e=\frac{q}{m_{0}}.$ (107)

Without loss of generality, we choose $M=1$ and $e\leq 1,Q\leq 1$. The
function $\varepsilon_{min}(r)=\frac{E_{min}(r)}{m_{0}}$ is plotted in Fig.2,
which show that $\varepsilon_{min}$ grows as $e$ increasing. As a numerical
calculation trick, we don’t have to solve the equation $\Psi_{l_{+}}=m_{0}$ to
derive $r_{1}$, but choose the maximum of
$\max\\{\frac{\Psi_{l_{+}}}{m_{0}},1\\}$ to obtain $\varepsilon_{min}$.

Figure 3: Sample graphs of $J_{t}/(-2\pi Am_{0}^{4})$ (up) and $J_{r}/(2\pi
Am_{0}^{4})$ (down). The parameters are chosen to be $Q=0.1,e=0.1,\mu=0,z=1$.
The vertical dotted line represents the location of the photon sphere.

Figure 4: Sample graphs of $T_{tt}/(2\pi Am_{0}^{5})$ (up), $T_{tr}/(-2\pi
Am_{0}^{5})$ (middle) and $T_{rr}/(2\pi Am_{0}^{5})$ (down). The parameters
are chosen to be $Q=0.1,e=0.1,\mu=0,z=1$.

Using the function $\varepsilon_{min}$, we numerical calculate $J_{t}$ and
$J_{r}$ as the function of $r$, see Fig.3. The absorption part and the
scattering part are calculated separately. The scattering part vanishes when
$r$ enters inside the photon sphere. The absorption portion is also
discontinuous in the photon sphere. However, the total part is smooth in the
photon sphere. Since we using the Boyer-Lindquist coordinate system, there is
a singularity at the event horizon. The numerical calculation is performed
outside the horizon. When $r\rightarrow r_{+}$, Fig.3 shows that $J_{r}$
approaches infinity. However, as we will see later, this does not affect the
calculation of the accretion rate. We also plot $T_{tt}$, $T_{tr}$ and
$T_{rr}$ as the function of $r$ in Fig.4. Their behavior is similar to
$J_{\mu}$.

Figure 5: Sample graphs of $\frac{n}{n_{\infty}}$. The parameters are chosen
to be $Q=0.2,e=0.6,\mu=0,z=1$.

The particle number density $n$ is expressed as Eq.(78). Using the numerical
results of $J_{\mu}$, we plot the particle number density
$\frac{n}{n_{\infty}}$ as the function of $r$ in Fig.5. As $z$ increases, the
function $\frac{n}{n_{\infty}}$ grows near the horizon and still approaches 1
far away. A similar behavior is also shown by Ref. Rioseco2017-1 ; Cieslik2020
.

###  9.2 the cosmic censorship and accretion rates

We know that the outer horizon $r_{+}$ of a Reissner-Nordström black hole
requires that the mass of the black hole $M$ should be bigger than the charge
$Q$. Say there are $N_{0}$ Fermions fall into the black hole, therefore the
mass becomes $M+N_{0}m_{0}$ and the charge becomes $Q+N_{0}q$, where $M$ is
the original mass and $Q$ is the original charge of the black hole. Then, the
charge-mass ratio of a black hole $\frac{Q+N_{0}q}{M+N_{0}m_{0}}$ will change
with the number of particles falling into the black hole. For those particles
with $e\leq 1$ keeps falling into the black hole, the final charge-mass ratio
of a black hole would inevitably turn into $e\leq 1$. For those particles with
$e>1$, Wald’s research Wald1974  shows that the accretion process, which
destroys a black hole’s event horizon, will be stopped. As is pointed by
Ref.Wald1974  that the particles with $e>1$ will not be captured by the
extreme Reissner-Nordström black hole $Q=M$. This is because that in the case
of $e>1$, the electrostatic repulsion becomes greater than the gravitational
attraction and the particle will be rejected. One can shooting it in with high
velocity. But by doing so one increases the particle’s energy. Until the
energy is raised to $E>q_{0}$ which cannot violate the relationship $Q\leq M$,
particles are to enter the black hole. Therefore, the naked singularity can
never happened. To keep the horizon existing, we obtain the maximum number
$N_{m}$ for the black hole accretion Fermions, which is

$N_{m}=\left\\{\begin{array}[]{ll}\frac{M-Q}{q-m_{0}},&\hbox{for
$q>m_{0}$;}\\\ \infty,&\hbox{for $q\leq m_{0}$.}\end{array}\right.$ (108)

This analysis shows that the accretion of Fermions onto a black hole would
decrease as the black hole becomes more charged, which is shown in Fig.6.

Figure 6: Sample graphs of $\frac{\dot{n}}{n_{\infty}}$ as the function of $Q$
(left) and as the function of $e$ (right). The left parameters are chosen to
be $e=0.1,\mu=0$, the right parameters are chosen to be $Q=0.5,\mu=0$.

We plot the particle accretion rate $\frac{\dot{n}}{n_{\infty}}$ in Fig.6. The
left graph shows exactly $\dot{n}$ decrease as the charge $Q$ grows. The right
graph also shows the accretion rate $\dot{n}$ decreases as $e$ increases. As
$z$ increases, the accretion rate $\dot{n}$ increases. However, the accretion
rate $\dot{n}$ depends little on the parameter $\mu$. The reason is that the
parameter $\mu$ appears in the denominator of the integral function
$\frac{1}{e^{z\varepsilon-\mu\beta}+1}$. The constant term $e^{\mu\beta}$ can
be brought outside the integral and left the integral function to
$\frac{1}{e^{z\varepsilon}+e^{\mu\beta}}$. The constant term $e^{\mu\beta}$
appears both in the numerator $\dot{n}$ and the denominator $n$. Thus, the
effect of changing the parameter $\mu$ is equivalent to changing the $+1$ term
in the distribution function (36).

Figure 7: Sample graphs of $\frac{\dot{E}}{n_{\infty}}$. The parameters are
chosen to be $e=0.1,\mu=0$.

We plot the energy accretion rate $\frac{\dot{E}}{n_{\infty}}$ in Fig.7. The
energy accretion rate $\dot{E}$ is not monotonically dependent on $z$. As $z$
increases, $\dot{E}$ decreases to minimum and then increases. We also find
that the accretion rate $\dot{E}$ decreases as $Q$ increases.

###  9.3 energy density, radial and tangential pressures

Figure 8: Sample graphs of $\frac{\rho}{nm_{0}c^{2}}$. The parameters are
chosen to be $Q=0.1,e=0.1,\mu=0$.

Figure 9: Sample graphs of $\frac{p_{rad}}{nm_{0}c^{2}}$. The parameters are
chosen to be $Q=0.1,e=0.1,\mu=0$.

Figure 10: Sample graphs of $\frac{p_{tan}}{nm_{0}c^{2}}$. The parameters are
chosen to be $Q=0.1,e=0.1,\mu=0$.

Figure 11: A comparison among the radial, tangential pressure and
$\frac{T_{rr}}{nm_{0}c^{2}}$ for two different values of $z$.

In Fig.8, we show the behavior of energy density $\frac{\rho}{nm_{0}c^{2}}$ as
the function of $r$. Since we chose physical coordinates, one can more easily
see how the physical quantity changes with the radius. The energy density
$\frac{\rho}{nm_{0}c^{2}}$ changes similarly as the figure in Ref.
Rioseco2017-2 . However, that seems no true for the radial pressure
$\frac{p_{rad}}{nm_{0}c^{2}}$, see Fig.9. Fig.10 show that the tangential
pressures $\frac{p_{tan}}{nm_{0}c^{2}}$ is consistent with Ref. Rioseco2017-2
. However, the radial pressure $p_{rad}$ is no longer equal to the tangential
pressures $p_{tan}$ at infinity, see Fig.11. As is shown in the previous
section, $T^{r}{}_{r}$ equals to $p_{tan}$. And they are components of the
energy-momentum of the isotropic perfect fluid for Fermi gas at infinity. Our
numerical results show that the radial pressure $p_{rad}$ is always smaller
than tangential pressures $p_{tan}$ even at infinity. As $z$ growth, the
difference between $p_{rad}(r\rightarrow\infty)$ and
$p_{tan}(r\rightarrow\infty)$ would get smaller and smaller. Only in the low-
temperature limit $z\rightarrow\infty$, the contribution of null fluid
disappears and the model can be treated as the ideal fluid accretion model at
infinity.

##  10 conclusions

In this paper, we study the accretion of degenerate Fermi gas onto a Reissner-
Nordström black hole. The Boyer-Linquist coordinate is chosen to construct the
accretion theory. We assume that the Fermi gas is in thermal equilibrium at
infinity and described by Fermi-Dirac statistic. Each particle interacts with
others only through the gravitational and electromagnetic field of the
background and travels along the orbit described by the abbreviated action. We
obtain the expression both for the absorption part and scattering part of the
particle current density and the stress energy-momentum tensor. The expression
of particle accretion rate and energy accretion rate are also obtained. At
large distances, we analytically calculate the corresponding quantities. In
the case of Maxwell-Jüttner distribution, they are consistent with the Ref.
Rioseco2017-1 ; Cieslik2020 . At finite range, we numerically calculate three
accretion rates, energy density, radial pressure and tangential pressure.

The Vlasov gas accretion model describes an anisotropic fluid consisting of
two perfect-fluid components. One component is an isotropic perfect fluid of
Fermi gas, and the other is a null fluid. When using Boyer-Lindquist
coordinates, the contribution of the null fluid remains at infinity. At this
point, the radial pressure is smaller than the tangential pressure. This
difference decreases as the temperature decreases and eventually vanishes as
the temperature approaches 0.

The results also show that the chemical potential had minimal effect on the
results. The more charges of the black hole, the lower the accretion rates of
a black hole. However, in the extremal case $Q=M$, the black hole would not
accrete particles with $e>1$. The charge-mass ratio of a black hole would
never be larger than 1. Thus, the naked singularity cannot form.

Acknowledgement This work is partially supported by the National Natural
Science Foundation of China (NSFC U2031112) and the Scientific Research
Foundation of Hunan University of Arts and Sciences (23BSQD237). We also
acknowledge the science research grants from the China Manned Space Project
with NO. CMS-CSST-2021-A06.

## References

  * (1) P. Rioseco and O. Sarbach, Accretion of a relativistic, collisionless kinetic gas into a Schwarzschild black hole, Classical Quantum Gravity 34, 095007 (2017).
  * (2) The Event Horizon Telescope Collaboration, First M87 event horizon telescope results. I. The shadow of the supermassive black hole, Astrophys. J. Lett. 875, L1 (2019).
  * (3) The Event Horizon Telescope Collaboration, First M87 event horizon telescope results. V. Physical origin of the asymmetric ring, Astrophys. J. Lett. 875, L5 (2019).
  * (4) The Event Horizon Telescope Collaboration, First Sagittarius A* event horizon telescope results. I. The shadow of the supermassive black hole in the center of the Milky Way, Astrophys. J. Lett. 930, L12 (2022).
  * (5) The Event Horizon Telescope Collaboration, First Sagittarius A* event horizon telescope results. V. Testing astrophysical models of the galactic center black hole, Astrophys. J. Lett. 930, L16 (2022).
  * (6) P. Rioseco and O. Sarbach, Spherical steady-state accretion of a relativistic collisionless gas into a Schwarzschild black hole, J. Phys. Conf. Ser. 831, 012009 (2017).
  * (7) F. Hoyle and R. A. Lyttleton, The effect of interstellar matter on climatic variation, Math. Proc. Cambridge Philos. Soc. 35, 405 (1939).
  * (8) R. Lyttleton and F. Hoyle, The evolution of the stars, Nature (London) 146, 97 (1940).
  * (9) H. Bondi and F. Hoyle, On the mechanism of accretion by stars, Mon. Not. R. Astron. Soc. 104, 273 (1944).
  * (10) E. Chaverra and O. Sarbach, Radial accretion flows on static spherically symmetric black holes, Classical Quantum Gravity 32, 155006 (2015).
  * (11) A. Cieślik and P. Mach, Accretion of the Vlasov gas on Reissner-Nordström black holes, Phys. Rev. D 102, 024032 (2020).
  * (12) P. Mach and A. Odrzywołek, Accretion of Dark Matter onto a Moving Schwarzschild Black Hole: An Exact Solution, Phys. Rev. Lett. 126, 101104 (2021).
  * (13) P. Mach and A. Odrzywołek, Accretion of the relativistic Vlasov gas onto a moving Schwarzschild black hole: Exact solutions, Phys. Rev. D 103, 024044 (2021).
  * (14) P. Mach and A. Odrzywołek, Accretion of the relativistic vlasov gas onto a moving schwarzschild black hole: Low temperature limit and numerical aspects, Acta Phys. Pol. B Proc. Suppl. 15, 1 (2022).
  * (15) A. Gamboa, C. Gabarrete, P. Domínguez-Fernández, D. Núñez, and Olivier Sarbach, Accretion of a Vlasov gas onto a black hole from a sphere of finite radius and the role of angular momentum, Phys. Rev. D 104, 083001 (2021).
  * (16) J. Liao and D. Liu, Accretion of the relativistic Vlasov gas onto a Bardeen regular black hole, Astrophysics and Space Science 367, 109 (2022).
  * (17) Z. Cai and R. Yang, Accretion of the Vlasov gas onto a Schwarzschild-like black hole, Phys. Dark Universe 42, 101292 (2023).
  * (18) P. Rioseco and O. Sarbach, Phase space mixing in the equatorial plane of the Kerr black hole, Phys. Rev. D 98, 124024 (2018).
  * (19) L. Andersson, P. Blue, and J. Joudioux, Hidden symmetries and decay for the Vlasov equation on the Kerr spacetime, Commun. Partial Differ. Equations 43, 47 (2018).
  * (20) A. Cieślik , P. Mach , and A. Odrzywołek, accretion of the relativistic Vlasov gas in the equatorial plane of the Kerr black hole, Phys. Rev. D 106, 104056 (2022).
  * (21) P. Mach and A. Odrzywołek, equatorial accetion on the Kerr black hole, Acta Phys. Pol. B Proc. Suppl. 16, 13 (2023).
  * (22) Ping Li, Yong-qiang Liu and Xiang-hua Zhai, Accretion of the relativistic Vlasov gas onto a Kerr black hole, Phys. Rev. D 108, 124022 (2023).
  * (23) P. Mach, A. Cieślik and A. Odrzywołek, Monte Carlo methods for stationary solutions of general-relativistic Vlasov systems: Collisionless accretion onto black holes, Phys. Rev. D 108, 124057 (2023).
  * (24) A. Cieślik, P. Mach and A. Odrzywołek, Monte Carlo methods for stationary solutions of general-relativistic Vlasov systems: Planar accretion onto a moving Schwarzschild black hole, arXiv:2406.04471v1.
  * (25) Fabio Pacucci and Abraham Loeb, Accretion from Winds of Red Giant Branch Stars May Reveal the Supermassive Black Hole in Leo I, The Astrophysical Journal Letters 940, 2 (2022).
  * (26) Kareem El-Badry, Hans-Walter Rix, Yvette Cendes et. al., A red giant orbiting a black hole, Monthly Notices of the Royal Astronomical Society 521, 4323–4348 (2023).
  * (27) S. Deheuvels, J. Ballot, C. Gehan and and B. Mosser, Seismic signature of electron degeneracy in the core of red giants: hints for mass transfer between close red-giant companions, Astronomy & Astrophysics 659, A106 (2022).
  * (28) Robert Wald, Gedanken Experiments to Destroy a Black Hole, Annals of physics 83, 548-556 (1974).
  * (29) S. Chandrasekhar, The Mathematical Theory of Black Holes (Clarendon Press, Oxford, 1983), P345-P348.
  * (30) Luciano Rezzolla and Olindo Zanotti, Relativistic Hydrodynamics (Oxford University Press, Oxiford, 2013), P90-P101.
  * (31) Partricio S. Letelier, Anisotropic fluids with two-perfect-fluid components, Phys. Rev. D 22, 807-813 (1980).


# The origin of ferromagnetic interactions in NaMnCl3: how the response theory
reconciles with Goodenough-Kanamori-Anderson rules

I. V. Solovyev<EMAIL_ADDRESS>National Institute for Materials
Science, MANA, 1-1 Namiki, Tsukuba, Ibaraki 305-0044, Japan  Institute of
Metal Physics, S. Kovalevskaya str. 18, 620108 Ekaterinburg, Russia  A. V.
Ushakov  Institute of Metal Physics, S. Kovalevskaya str. 18, 620108
Ekaterinburg, Russia  S. V. Streltsov  Institute of Metal Physics, S.
Kovalevskaya str. 18, 620108 Ekaterinburg, Russia

###### Abstract

The on-site Coulomb repulsion $U$ is the key ingredient for describing the
magnetic properties of Mott insulators, leading to a popular believe that many
limitations of the density-functional theory based methods can be cured by
artificially incorporating such on-site interactions for localized electrons
in the model form. The layered antiferromagnet NaMnCl3 reveals quite a
different story: while the Coulomb $U$ on the Mn sites controls the strength
of antiferromagnetic superexchange interactions, an equally important
parameter is the Stoner coupling ${\cal I}_{\rm Cl}$ on the ligand sites. The
latter is responsible for large ferromagnetic contributions to the interatomic
exchange interactions, which in NaMnCl3 nearly cancel the effect of the
superexchange interactions. Although such behavior is anticipated from the
phenomenological Goodenough-Kamanori-Anderson rules, the quantitative
description of the ligand-related contributions remains disputable.
Considering NaMnCl3 as an example, we discuss how they can be generally taken
into account in the linear response theory to regain the dependence of the
exchange interactions on ${\cal I}_{\rm Cl}$. The problem is complicated by
the fact that, for the nearly filled Cl $3p$ shell, the parameters ${\cal
I}_{\rm Cl}$ are sensitive to the model assumptions.

Introduction. – In most cases, the exchange interactions between half-filled
$3d$ ions are expected to be antiferromagnetic due to the superexchange
mechanism [1]. Such situation is indeed realized in the canonical MnO [2], the
multiferroic BiFeO3 [3] and MnI2 [4], the lithium-ion batteries cathode
LiMnPO4 [5], the flocculant compound FeCl3 [6], and other prominent magnetic
materials composed from the Mn2+ or Fe3+ ions. However, when the ${\rm
T}$-${\rm L}$-${\rm T}$ exchange path between two transition-metal (${\rm T}$)
sites connected by a ligand (${\rm L}$), becomes close to $90^{\circ}$, the
situation may be less certain as there are several mechanisms operating in
opposite directions and partially cancelling each-other. Besides
antiferromagnetic superexchange, such mechanisms typically involve
ferromagnetic Hund’s coupling on the ligand states. In fact, by summarising
the famous Goodenough-Kamanori-Anderson (GKA) rules for transition-metal
compounds, Junjiro Kanamori has concluded that the sign of exchange
interaction in this $90^{\circ}$-case is “uncertain” [7].

NaMnCl3, which was recently reinvestigated by Devlin and Cava [8], is the good
example of such uncertainty. It crystallizes in the trigonal structure with
the space group $R\overline{3}$ (Fig. 1a), which is similar to that of the van
der Waals ferromagnet CrI3 [9], except that the [MnCl3]- layers are negatively
charged and interconnected by the Na+ layers. The Mn2+ ions form the distorted
honeycomb lattice and the Mn-Cl-Mn angle is $96.7^{\circ}$ [10]. NaMnCl3
orders antiferromagneticaly below $T_{\rm N}\approx 6.3$ K [8, 11].
Nevertheless, the Curie-Weiss constant is positive, $\theta_{\rm CW}\approx
4.2$ K, suggesting that some of the exchange interactions are ferromagnetic.
This finding is further elaborated by powder neutron diffraction measurements
indicating at the layered antiferromagnetic (LA) structure of NaMnCl3 in which
the ferromagnetic layers are coupled antiferromagnetically [12]. In the light
of these delicate magnetic properties, NaMnCl3 appears to be an interesting
system for testing abilities of first-principles methods based on the density-
functional theory (DFT).

Figure 1: (a) Crystal structure of NaMnCl3: top view (top) and side view
(bottom). The unit cell in the basal plane is shown by blue solid lines. Two
Mn sites in the primitive cell are denoted by different colors. (b) Partial
densities of Mn $3d$ and Cl $3p$ states, as obtained in the generalized
gradient approximation for the antiferromagnetic state. The zero energy is in
the middle of the gap between occupied and empty states. The positions of main
bands are indicated by symbols.

Construction of the spin models inevitably relies on some assumptions. One of
them is the use of magnetic force theorem (MFT) for interatomic exchange
interactions [13], which becomes a popular tool in the DFT community [14, 15,
16, 17]. However, apart from the fundamental limitations of MFT [18, 19, 20],
most of such calculations consider the exchange processes only between
magnetic $3d$ states, where ligands play a role of an effective medium, which
only helps to transfer electrons from one ${\rm T}$-site to another. The
contributions caused by the magnetic polarization of the ligand states,
$m^{z}_{\rm L}$, are typically ignored or treated in a phenomenological way,
by associating them with direct exchange integrals [21, 22, 23] or considering
an ad hoc energy gain caused by this $m^{z}_{\rm L}$ [24, 25]. In this
article, by employing formally exact technique [20], we explicitly show how
the ligand-related contributions to the exchange parameters can be naturally
evaluated in the framework of the linear response theory. We further
demonstrate that, if defined properly, such exchange interactions depend not
only on the Coulomb repulsion $U$ on the Mn sites, but also on the Stoner
coupling on the ligand sites, ${\cal I}_{\rm Cl}$, responsible for Hund’s
first rule.

Theory. – In order to evaluate the exchange interactions in DFT, we use the
exact expression for the total energy change caused by infinitesimal rotations
of spins, which is formulated in terms of the inverse response function [18,
19, 20]. This energy change can be mapped onto the spin model

${\cal
E}=-\frac{1}{2}\sum_{ij}J_{ij}\boldsymbol{e}_{i}\cdot\boldsymbol{e}_{j}$ (1)

with $\boldsymbol{e}_{i}$ being the unit vector along the direction of spin at
the $i$th ${\rm T}$-site. The corresponding exchange interactions in the
reciprocal space include two contributions [20]: $J_{\rm
TT}(\boldsymbol{q})=J_{\rm TT}^{0}(\boldsymbol{q})+\Delta J_{\rm
TT}(\boldsymbol{q})$, where $J_{\rm TT}^{0}(\boldsymbol{q})$ are the bare
interactions, while

$\Delta J_{\rm TT}(\boldsymbol{q})=-J_{\rm TL}^{0}(\boldsymbol{q})\left[J_{\rm
LL}^{0}(\boldsymbol{q})\right]^{-1}J_{\rm LT}^{0}(\boldsymbol{q})$ (2)

takes into account the contributions of the ligand states. In these notations,
the subscripts ${\rm TT}$, ${\rm TL}$, ${\rm LL}$, and ${\rm LT}$ mean that
all the quantities are the matrices with the indices ($\mu$ and $\nu$) running
over the ${\rm T}$\- and ${\rm L}$-sites in the unit cell. Such $\Delta J_{\rm
TT}(\boldsymbol{q})$ is obtained in the adiabatic limit, where the ${\rm
L}$-spins instantaneously follow each configuration of the ${\rm T}$-spins
[20]. The matrix elements of $J_{\mu\nu}^{0}(\boldsymbol{q})$ ($\mu,\nu\in{\rm
T}$ or ${\rm L}$) are expressed in terms of the response tensor
${\mathbb{R}}^{\uparrow\downarrow}(\boldsymbol{q})$ as

$J_{\mu\nu}^{0}(\boldsymbol{q})=\frac{1}{2}m_{\mu}^{z}\left(\left[{\mathbb{R}}^{\uparrow\downarrow}(\boldsymbol{q})\right]^{-1}_{\mu\nu}+{\cal
I}_{\mu}\delta_{\mu\nu}\right)m_{\nu}^{z},$ (3)

where $m_{\mu}$ is the magnetic moment and the Stoner parameter ${\cal
I}_{\mu}$ corresponds to the exchange-correlation energy of the form
$-\frac{1}{4}{\cal I}_{\mu}\big{(}m_{\nu}^{z}\big{)}^{2}$. Further details can
be found in Refs. [20, 26]. The real space parameters $J_{ij}$ are obtained by
the Fourier transform of $J_{\mu\nu}(\boldsymbol{q})$. The popular alternative
to this, formally exact, approach is MFT [13], which relies on additional
approximations and can be justified only in the long wavelength and strong-
coupling limits [20]. Furthermore, the contribution $\Delta J_{\rm
TT}(\boldsymbol{q})$ is ignored in most of the MFT based calculations or
replaced by a semi-empirical term simulating the direct exchange [21, 22, 23].
However, $\Delta J_{\rm TT}(\boldsymbol{q})$ can be very important, especially
when it is consideerd in the combination with the exact Eq. (3) for the
exchange interactions [20].

Results and Discussions. – The electronic structure calculations were
performed in the generalized gradient approximation (GGA) [27], for the
experimental structure reported in Ref. [10] using Vienna ab initio simulation
package (VASP) [28]. After that, several tight-binding (TB) models were
constructed in the maximally localized Wannier basis [29].

The electronic structure of NaMnCl3 is relatively simple, consisting of
several groups of isolated bands (Fig. 1b), which allows us to construct two
types of models. The simplest $d$-model includes ten Mn $3d$ bands located
near the Fermi level. In this case, the Wannier basis consists of five $3d$
orbitals (labeled by $m$ and $m^{\prime}$) at each Mn site. The corresponding
hopping parameters $\hat{t}_{ij}=[t_{ij}^{mm^{\prime}}]$, are associated with
the matrix elements of the GGA Hamiltonian in the Wannier basis [29]. Then,
$J_{ij}$ can be evaluated in the superexchange approximation for the half-
filled Mn $3d$ shell as: $J_{ij}=-\frac{{\rm
Tr}_{L}(\hat{t}_{ij})^{2}}{U+4J_{\rm H}}$ [1], where $U$ and $J_{\rm H}$ are
the intraatomic Coulomb repulsion and exchange coupling, respectively, and
${\rm Tr}_{L}$ denotes the trace over the orbital indices. In the mean-field
approximation, ${\cal I}_{\rm Mn}$ is related to $U$ and $J_{\rm H}$ as ${\cal
I}_{\rm Mn}=\frac{U+4J_{\rm H}}{5}$ [30, 31]. We treat ${\cal I}_{\rm Mn}$ as
an adjustable parameter and ask which ${\cal I}_{\rm Mn}$ would reproduce the
experimental $T_{\rm N}=6.3$ K [8] for the given hopping parameters
$\hat{t}_{ij}$. In order to evaluate $T_{\rm N}$, we use the random phase
approximation (RPA), as explained in Ref. [32]. This yields ${\cal I}_{\rm
Mn}\sim 5.5$ eV, which is unrealistically large and, assuming $J_{\rm H}=0.8$
eV [33], would correspond to nearly bare value of $U\sim 24$ eV. Furthermore,
all interactions $J_{ij}$ are antiferromagnetic, which cannot explain the sign
of $\theta_{\rm CW}$.

Thus, it is important to include the Cl $3p$ band explicitly into the model.
For these purposes, we construct the TB Hamiltonian in the basis of Mn $3d$
and Cl $3p$ Wannier orbitals, which reproduces the behavior of Mn $3d$ and Cl
$3p$ bands in GGA (Fig. 1b). In order to evaluate the exchange interactions
using Eq. (3), we focus on the antiferromagnetic (A) solution with two Mn
spins in the unit cell being antiparallel to each other. The ferromagnetic (F)
solution was also investigated and confirmed to provide a similar behavior for
the exchange interactions. First, we considered the contributions of the Mn
$3d$ states alone, by artificially cutting off all matrix elements of
${\mathbb{R}}^{\uparrow\downarrow}(\boldsymbol{q})$ associated with the ${\rm
L}$ sites, as is typically done in calculations based on MFT. In comparison
with the $d$-model, the additional admixture of the Mn $3d$ states into the Cl
$3p$ band in the $dp$-model already plays an important role by substantially
weakening the antiferromagnetic interactions. For instance, the nearest-
neighbor (nn) interaction $J_{1}$ is reduced from $-$$14.4$ till $-$$6.9$ meV
(in GGA, see Fig. 2 for the notations), and similar tendencies are found for
other interactions. The complete removal of the Cl $3p$ states (by
artificially shifting them downwards by $100$ Ry) yields $J_{1}=-$$4.5$ meV,
which can be regarded as the direct Mn-Mn interaction. Therefore, the
remaining Mn-Cl-Mn superexchange interactions operating via the Cl $3p$ sites
can be estimated in the $d$\- and $dp$-model as $-$$9.9$ meV and $-$$2.4$ meV,
respectively.

Figure 2: (a) Distance dependence of exchange interactions (in GGA) for three
values of ${\cal I}_{\rm Cl}$. The inset shows the dependence of $J_{1}$ on
${\cal I}_{\rm Cl}$. (b) Notations of main exchange interactions.

Up to this point, we obtain a good agreement between the exact theory and MFT
[13]. For instance, MFT yields $J_{1}=-$$14.1$ and $-$$6.5$ meV in the $d$\-
and $dp$-model, respectively, which are indeed very close to the “exact”
values. This is not surprising, because MFT is justified in the strong
coupling limit [20], which is realized for the Mn $3d$ states in NaMnCl3 even
at the level of GGA. However, as long as we deal only with the Mn $3d$ states,
all interactions remain antiferromagnetic, which is again inconsistent with
the type of the magnetic ground state of NaMnCl3 (A instead of LA) as well as
the sign of $\theta_{\rm CW}$.

The effect of the Cl $3p$ states on the exchange interactions is very
dramatic. The Cl $3p$ states are located rather deep in the occupied part and
their hybridization with the Mn $3d$ states is not particularly strong. In
such a situation, the sub-block
$\big{[}{\mathbb{R}}^{\uparrow\downarrow}(\boldsymbol{q})\big{]}_{\rm LL}$ of
the response matrix is small (and would completely vanish in the ionic limit,
where there is no hybridization with the Mn $3d$ states). Therefore, the
inverse matrix
$\left[{\mathbb{R}}^{\uparrow\downarrow}(\boldsymbol{q})\right]^{-1}$ will
nearly diverge and in the calculations of
$J_{ij}^{\phantom{0}}=J_{ij}^{0}+\Delta J_{ij}^{\phantom{0}}$ we have to deal
with large quantities. For instance, for the nn interaction, we obtain (within
GGA): $J_{1}^{0}=-$$100.0$ meV and $\Delta J_{1}^{\phantom{0}}=122.2$ meV,
which strongly cancel each other. Obviously, once the Cl $3p$ states are
explicitly included into the model, the strong coupling limit is not justified
and MFT cannot be used: in GGA, MFT yields $J_{1}^{0}=-$$7.6$ meV and $\Delta
J_{1}^{\phantom{0}}=-$$1.1$ meV, which differ drastically from the parameters
obtained in the exact approach. Furthermore, there is no more cancellation
between $J_{1}^{0}$ and $\Delta J_{1}^{\phantom{0}}$ within MFT.

Then, we stick to the exact approach and search for the condition, when such
cancellation would reproduce the experimental behavior of NaMnCl3, namely the
LA ground states with small $T_{\rm N}$ and small positive $\theta_{\rm CW}$.
First of all, we have realized that, although the parameters ${\cal I}_{\mu}$
on the ${\rm T}$-sites are well defined, the ones on the ${\rm L}$-sites
appear to be not. Generally, ${\cal I}_{\mu}$ is related to the exchange field
$b_{\mu}^{z}$ and magnetization $m_{\mu}^{z}$ as ${\cal
I}_{\mu}=-\frac{b_{\mu}^{z}}{m_{\mu}^{z}}$. However, there are several
possibilities for evaluating these quantities in the TB model. For instance,
$b_{\mu}^{z}$ can be associated with the site-diagonal elements of the TB
Hamiltonian in the Wannier basis. Such
$\hat{b}_{\mu}^{z}=\hat{H}_{\mu\mu}^{\uparrow}-\hat{H}_{\mu\mu}^{\downarrow}$
is the matrix in the subspace of the orbital indices. Then, in the spirit of
DFT [34], corresponding site-diagonal magnetization matrix $\hat{m}_{\mu}^{z}$
should be obtained from the occupied eigenstates of $\hat{H}$. Hence, ${\cal
I}_{\mu}$ can be defined as ${\cal I}_{\mu}=-\frac{{\rm
Tr}_{L}(\hat{b}_{\mu}^{z})}{{\rm Tr}_{L}(\hat{m}_{\mu}^{z})}$. Another
possibility (used above) is to employ the sum rule,
$b_{\mu}^{z}=\sum_{\nu}[{\mathbb{R}}^{\uparrow\downarrow}(0)]_{\mu\nu}^{-1}m_{\nu}^{z}$
[20], in order to evaluate $b_{\mu}^{z}$ for the given set of magnetic moments
and then use again the equation ${\cal
I}_{\mu}=-\frac{b_{\mu}^{z}}{m_{\mu}^{z}}$. Nevertheless, we would like to
emphasize that both definitions rely on assumptions. Particularly, the sum
rule implies that the intersite matrix elements of the TB Hamiltonian
$\hat{H}$ do not depend on spin, which is an approximation. On the other hand,
the change $b_{\rm T}^{z}$ in the TB Hamiltonian, will affect not only $m_{\rm
T}^{z}$ but also $m_{\rm L}^{z}$, even for fixed $b_{\rm L}^{z}$, thus making
the definition ${\cal I}_{\rm L}=-\frac{b_{\rm L}^{z}}{m_{\rm L}^{z}}$
ambiguous.

The results are summarized in Table 1.

Table 1: Effective Stoner parameters (in eV), as obtained in GGA for the antiferromagnetic (ferromagnetic) state by using site-diagonal elements of the tight-binding Hamiltonian and the sum rule. method | ${\cal I}_{\rm Mn}$ | ${\cal I}_{\rm Cl}$
---|---|---
${\cal I}_{\mu}=-\frac{{\rm Tr}_{L}(\hat{b}_{\mu}^{z})}{{\rm Tr}_{L}(\hat{m}_{\mu}^{z})}$ |  $0.88$ ($0.88$) |  $1.40$ ($1.51$)
sum rule |  $0.91$ ($0.90$) |  $10.95$ ($8.75$)

${\cal I}_{\rm Mn}$ only weakly depends on the method and the magnetic state.
The obtained values ${\cal I}_{\rm Mn}=0.88-0.91$ eV are reasonable and
consistent with ${\cal I}_{\rm Mn}=0.82$ eV reported for metallic fcc Mn [35].
On the other hand, depending on the definition and the magnetic state, ${\cal
I}_{\rm Cl}$ varies from $1.4$ eV till $11$ eV (and even larger if one uses
GGA$+$$U$ instead of GGA). The upper boundary for ${\cal I}_{\rm Cl}$, derived
from the sume rule, may look unrealistically large. However, it should be
understood that this ${\cal I}_{\rm Cl}$ is proportional to
$\left[{\mathbb{R}}^{\uparrow\downarrow}(\boldsymbol{q})\right]^{-1}$, which
nearly diverges in the case of NaMnCl3. Furthermore, the concept of the Stoner
splitting does not make sense for the totally filled Cl $3p$ shell in the
ionic limit (contrary to the half-filled Mn $3d$ shell).

Thus, one can expect that the bare interaction $J_{ij}^{0}$, which does not
explicitly depend on the Stoner coupling, is well defined. On the other hand,
$\Delta J_{ij}^{\phantom{0}}$ is sensitive to the definition of ${\cal I}_{\rm
Cl}$. In the following, we treat ${\cal I}_{\rm Cl}$ as an adjustable
parameter, by varying it between the values obtained using the sum rule and
site-diagonal elements of $\hat{H}$. Furthermore, we simulate the effect of
the on-site $U$ by artificially changing the splitting between the majority
($\uparrow$) and minority ($\downarrow$) Mn $3d$ spin-states in the TB
Hamiltonian. For the half-filled shell, this is equivalent to the change of
the parameter ${\cal I}_{\rm Mn}$ and a constant shift of all Mn $3d$ states,
which is controlled by the double-counting term in the GGA+$U$ approach. Since
the form of this double counting is largely empirical, we consider the
phenomenological procedure by fixing the positions of the occupied
$\uparrow$-spin states and shifting upwards the unoccupied $\downarrow$-spin
states, which accounts for the experimental situation in MnO [31].

The basic idea is explained Fig. 2, in GGA. As expected, ${\cal I}_{\rm Cl}$
tends to strengthen the ferromagnetic interactions. For small ${\cal I}_{\rm
Cl}=1.4$ eV, all interactions are antiferromagnetic. Then, for the
intermediate ${\cal I}_{\rm Cl}=8.1$ eV, the nn interaction in the $xy$ plane,
$J_{1}$, becomes ferromagnetic, while the interactions $J_{2}$, $J_{4}$,
$J_{5}$, and $J_{7}$ between the planes remain antiferromagnetic. Finally, for
large ${\cal I}_{\rm Cl}=14.8$ eV, the interactions in and between the planes
tends to be ferromagnetic. Thus, one can obtain the following phase diagram in
the plane ${\cal I}_{\rm Mn}$-${\cal I}_{\rm Cl}$ (Fig. 3).

Figure 3: (a) Phase diagram in the plane ${\cal I}_{\rm Mn}$-${\cal I}_{\rm
Cl}$. (b) Ferromagnetic (F), layered antiferromagnetic (LA), and
antiferromagnetic (A) structures. (c) The Néel temperature ($T_{\rm N}$) and
Curie-Weiss temperature ($\theta_{\rm CW}$) as the function of ${\cal I}_{\rm
Cl}$ for three characteristic values of ${\cal I}_{\rm Mn}$, which are
indicated by the dot-dashed lines in (a). #1 corresponds to GGA data.

For each ${\cal I}_{\rm Mn}$ ($U$), the system undergoes the transition from
the A to LA state and, then, to the F state with the increase of ${\cal
I}_{\rm Cl}$. The parameters ${\cal I}_{\rm Cl}$ required to stabilize the F
state are probably too large. Nevertheless, the A-LA transition takes place at
more realistic values of ${\cal I}_{\rm Cl}$ ($\sim$$3.1$ in GGA, which
further increases with the increase of $U$). Then, using the obtained
parameters $J_{ij}$, we evaluate $T_{\rm N}$ (in RPA [32]) and $\theta_{\rm
CW}=\frac{1+1/S}{3k_{\rm B}}\sum_{i}J_{i}$ (where $S=5/2$). The results are
summarized in Fig. 3c. We see that, in order to reproduce the experimental
behavior of NaMnCl3 (small $T_{\rm N}$ and small positive $\theta_{\rm CW}$),
it is necessary to stay in the LA region, but close to the A-LA border. For
instance, for ${\cal I}_{\rm Mn}=$ $0.8$ (the GGA value), $1.4$, and $1.9$ eV,
such a situation is realized around ${\cal I}_{\rm Cl}\sim$ $3.6$, $6.0$, and
$11.0$ eV, respectively, yielding $(T_{\rm N},\theta_{\rm CW})=$ $(11.8,1.5)$,
$(5.0,3.0)$, and $(3.1,3.2)$ K. Assuming $J_{\rm H}=0.8$ eV, ${\cal I}_{\rm
Mn}=$ $0.8$, $1.4$, and $1.9$ eV, would correspond to the Coulomb $U=$ $0.9$,
$3.6$, and $6.3$ eV, respectively. In all three cases, the required values of
${\cal I}_{\rm Cl}$ fall in-between two estimates based on the site-diagonal
elements of $\hat{H}$ and the sum rule. Thus, irrespectively on the value of
$U$, the experimental behavior of NaMnCl3 can be formally reproduced by tuning
${\cal I}_{\rm Cl}$.

Then, we turn to the brute force GGA$+$$U$ calculations with the parameters
$U=7.0$ eV and $J_{\rm H}=0.8$ eV, which are typically extracted from
constrained DFT calculations [33]. The main difference from the previous
analysis is the double counting term, which is now taken in the standard form
to reproduce the ionization potential and electron affinity of the Mn2+ ions
in the atomic limit [36]. This leads to the additional downward shift of the
Mn $3d$ states, so that the occupied $\uparrow$-spin Mn $3d$ states are now
located below the Cl $3p$ band (see Fig. 4a). Nevertheless, as far as the
exchange interactions are concerned, this shift plays only a secondary role as
the magnetic properties are mainly controlled by ${\cal I}_{\rm Cl}$, which
strongly depends on the definition and varies from $2.5$ eV, if one uses site-
diagonal elements of $\hat{H}$, till $45.5$ eV, if one uses the sum rule. The
A-LA transition occurs near the lower boundary and the experimental behavior
is captured reasonably well by ${\cal I}_{\rm Cl}\approx 3.6$ eV, yielding
$T_{\rm N}=12.1$ K and $\theta_{\rm CW}=4.6$ K.

Figure 4: Results of GGA$+$$U$ calculations: (b) Partial densities of Mn $3d$
and Cl $3p$ states. The zero energy is in the middle of the gap between
occupied and empty states. (b) The Néel temperature ($T_{\rm N}$) and Curie-
Weiss temperature ($\theta_{\rm CW}$) as the function of ${\cal I}_{\rm Cl}$.

Finally, taking into account the uncertainty with the choice of the parameters
${\cal I}_{\rm Cl}$ in the model, it is right to ask which is the realistic
estimate of ${\cal I}_{\rm Cl}$ and which method is more suitable for
capturing the experimental behavior of NaMnCl3? For these purposes we evaluate
effective exchange interactions by mapping the total energies obtained for
several magnetic configurations onto the Heisenberg model (1). We consider
five such configurations, which allow us extract three in-plane interactions
$J_{1}$, $J_{3}$, and $J_{6}$, and one out-of-plane interaction $J_{2}$. Then,
GGA (GGA$+$$U$) yields the following parameters (in meV): $J_{1}=-1.1$
($0.6$), $J_{2}=-1.3$ ($-0.2$), $J_{3}=-0.3$ ($0.1$), and $J_{6}=-0.6$
($0.2$). Thus, the ground state of NaMnCl3 is expected to be of the A (LA)
type with $T_{\rm N}=16.6$ ($10.8$) K and $\theta_{\rm CW}=-45.1$ ($13.3$) K
in GGA (GGA$+$$U$). Definitely, GGA$+$$U$ does a better job, while plain GGA
fails to reproduce the magnetic ground state of NaMnCl3, though only by a tiny
margin as the relatively small values of $T_{\rm N}$ and $\theta_{\rm CW}$
indicate that the situation is pretty close to the A-LA boundary. If we tried
to reproduce such $T_{\rm N}$ and $\theta_{\rm CW}$ in the linear response
theory, using Eq. (3) and the downfolding procedure for the ligand spins, we
would have to use ${\cal I}_{\rm Cl}\approx 2.6$ ($4.4$) eV in GGA
(GGA$+$$U$). These values are close to the lower boundaries of ${\cal I}_{\rm
Cl}$, evaluated via site-diagonal elements of $\hat{H}$, which still need to
be increased by about 80% in order to reproduce the results of the total
energy calculations.

Summary and Conclusions. – The response theory for interatomic exchange
interactions is the powerful modern tool as it provides not only quantitative
estimates, but also a deep microscopic understanding of the origin of these
interactions. Nevertheless, in order to describe these interactions properly
in NaMnCl3 and other $90^{\circ}$-exchange systems, it is essential to take
into account the polarization of the ligand states – the contribution, which
was typically ignored in most of the approximate techniques based on MFT. We
explain the origin of this contribution, which naturally incorporates the
dependence of $J_{ij}$ on not only the Coulomb $U$ on the Mn sites, but also
the Stoner coupling ${\cal I}_{\rm Cl}$ on the ligand sites, as was long
advocated by the phenomenological GKA rules. Yet, the main obstacle for truly
ab initio description along this line is the definition of the parameter
${\cal I}_{\rm Cl}$ for the nearly filled Cl $3p$ shell, which does not seem
to be unique.

As the next important development of the linear response technique, it would
be interesting to include the relativistic spin-orbit interaction, which is
responsible for anisotropic Dzyaloshinskii-Moriya interactions (DMI) and
exchange anisotropy. Particularly, our study clearly demonstrate the crucial
role of the ligand $p$ states in the exchange processes. However, there are
many examples when the same ligand states play a key role in developing the
large anisotropy and DMI, as demonstrated in CrI3 [37, 38] and Cr2Ge2Te6 [39].
Furthermore, the ligand assistant exchange anisotropy in honeycomb layered
cobaltites may lead to the realization of Kitaev spin liquid state, similar to
Na2IrO3 and $\alpha$-RuCl3 [40]. In this respect, NaCoCl3 seems to be
particularly promising candidate if it could be synthesized.

S.V.S would like to thank Robert Cava for stimulating discussions. The work
was supported by the project RSF 20-62-46047 of the Russian Science Foundation
in the part of conventional DFT calculations and the program 122021000038-7
(Quantum) of the Russian Ministry of Science and High Education in the part of
construction and analysis of microscopic models of interatomic exchange
interactions.

## References

  * [1] P. W. Anderson, Phys. Rev. 115, 2 (1959).
  * [2] M. E. Lines and E. D. Jones, Phys. Rev. 139, A1313 (1965).
  * [3] I. Sosnowska and A. K. Zvezdin, J. Magn. Magn. Matter. 140-144, 167 (1995); J. Jeong, E. A. Goremychkin, T. Guidi, K. Nakajima, G. S. Jeon, S.-A. Kim, S. Furukawa, Y. B. Kim, S. Lee, V. Kiryukhin, S. W. Cheong, and J.-G. Park, Phys. Rev. Lett. 108, 077202 (2012).
  * [4] J. W. Cable, M. K. Wilkinson, E. O. Wollan, and W. C. Koehler, Phys. Rev. 125, 1860 (1962).
  * [5] J. Li, W. Tian, Y. Chen, J. L. Zarestky, J. W. Lynn, and D. Vaknin, Phys. Rev. B 79, 144410 (2009).
  * [6] E. R. Jones Jr., O. B. Morton, L. Cathey, T. Auel, and E. L. Amma, J. Chem. Phys. 50, 4755 (1969); K. Ohhashi and I. Tsujikawa, J. Phys. Soc. Jpn 36, 980 (1974).
  * [7] J. Kanamori, J. Phys. Chem. Solids 10, 97 (1959).
  * [8] K. P. Devlin and R. J. Cava, Solid State Commun. 342, 114598 (2022).
  * [9] B. Huang, G. Clark, E. Navarro-Moratalla, D. R. Klein, R. Cheng, K. L. Seyler, D. Zhong, E. Schmidgall, M. A. McGuire, D. H. Cobden, W. Yao, D. Xiao, P. Jarillo-Herrero, and X. Xu, Nature (London) 546, 270 (2017).
  * [10] C. J. J. van Loon and G. C. Verschoor, Acta Cryst. B29, 1224 (1973).
  * [11] A. R. Lim, S. H. Choh, and S.-Y. Jeon, Phys. Stat. Sol. (b) 196, 425 (1996).
  * [12] N. Fedoseeva, I. Spevakova, G. Petrakovskii, V. Chuev and S. Petrov, J. Magn. Magn. Mater. 15-18, 539 (1980).
  * [13] A. I. Liechtenstein, M. I. Katsnelson, V. P. Antropov, and V. A. Gubanov, J. Magn. Magn. Mater. 67, 65 (1987).
  * [14] I. V. Solovyev, in _Recent Res. Devel. Magnetism & Magnetic Mat._, Vol. 1 (Transworld Research Network, Kerala, India, 2003), pp 253-294, ISBN: 81-7895-100-2; arXiv:cond-mat/0305668 (2003).
  * [15] Y. O. Kvashnin, O. Grånäs, I. Di Marco, M. I. Katsnelson, A. I. Lichtenstein, and O. Eriksson, Phys. Rev. B 91, 125133 (2015).
  * [16] Z. V. Pchelkina and S. V. Streltsov, Phys. Rev. B 88, 054424 (2013); Dm. M. Korotin, V. V. Mazurenko, V. I. Anisimov, and S. V. Streltsov, Phys. Rev. B 91, 224405 (2015).
  * [17] H. Yoon, T. J. Kim, J.-H. Sim, S. W. Jang, T. Ozaki, and M. J. Han, Phys. Rev. B 97, 125132 (2018).
  * [18] P. Bruno, Phys. Rev. Lett. 90, 087205 (2003).
  * [19] V. P. Antropov, M. van Schilfgaarde, S. Brink, and J. L. Xu, J. Appl. Phys. 99, 08F507 (2006).
  * [20] I. V. Solovyev, Phys. Rev. B 103, 104428 (2021).
  * [21] W. Ku, H. Rosner, W. E. Pickett, and R. T. Scalettar, Phys. Rev. Lett. 89, 167204 (2002).
  * [22] V. V. Mazurenko, S. L. Skornyakov, A. V. Kozhevnikov, F. Mila, and V. I. Anisimov, Phys. Rev. B 75, 224408 (2007).
  * [23] D. I. Badrtdinov, S. A. Nikolaev, M. I. Katsnelson, and V. V. Mazurenko, Phys. Rev. B 94, 224418 (2016).
  * [24] S. V. Streltsov and D. I. Khomskii, Phys. Rev. B, 77, 064405 (2008).
  * [25] I. Solovyev, J. Phys. Soc. Jpn. 78, 054710 (2009).
  * [26] I. V. Solovyev, S. A. Nikolaev, A. V. Ushakov, V. Yu. Irkhin, A. Tanaka, and S. V. Streltsov, Phys. Rev. B 105, 014415 (2022).
  * [27] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * [28] G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
  * [29] A. A. Mostofi, J. R. Yates, G. Pizzi, Y. S. Lee, I. Souza, D. Vanderbilt, and N. Marzari, Comput. Phys. Commun. 185, 2309 (2014).
  * [30] G. Stollhoff, A. M. Oleś, and V. Heine, Phys. Rev. B 41, 7028 (1990).
  * [31] I. V. Solovyev and K. Terakura, Phys. Rev. B 58, 15496 (1998).
  * [32] I. V. Solovyev and S. V. Streltsov, Phys. Rev. Materials 3, 114402 (2019).
  * [33] V. I. Anisimov, J. Zaanen, and O. K. Andersen, Phys. Rev. B 44, 943 (1991).
  * [34] W. Kohn and L. J. Sham, Phys. Rev. A 140, 1133 (1965).
  * [35] J. F. Janak, Phys. Rev. B 16, 255 (1977).
  * [36] I. Solovyev, N. Hamada, and K. Terakura, Phys. Rev. B 53, 7158 (1996).
  * [37] J. L. Lado and J. Fernández-Rossier, 2D Mater. 4, 035002 (2017).
  * [38] J. Liu, M. Shi, J. Lu, and M. P. Anantram, Phys. Rev. B 97, 054416 (2018).
  * [39] K. Wang, S. Nikolaev, W. Ren, and I. Solovyev, Phys. Chem. Chem. Phys. 21, 9597 (2019).
  * [40] H. Liu, J. Chaloupka, G. Khaliullin, Phys. Rev. Lett. 125, 047201 (2020).