#  Einstein–Vlasov system with equal-angular momenta in $\text{AdS}_{5}$

Hiroki
Asami111Email:mailto:asami.hiroki.b3<EMAIL_ADDRESS>Division of Particle and Astrophysical Science, Graduate School of Science,
Nagoya University, Nagoya 464-8602, Japan  Chul-Moon
Yoo222Email:mailto:yoo.chulmoon.k6<EMAIL_ADDRESS>Division of Particle and Astrophysical Science, Graduate School of Science,
Nagoya University, Nagoya 464-8602, Japan  Ryo
Kitaku333Email:mailto:kitaku.ryo.f4<EMAIL_ADDRESS>Division of Particle and Astrophysical Science, Graduate School of Science,
Nagoya University, Nagoya 464-8602, Japan  Keiya
Uemichi444Email:mailto:uemichi.keiya.j4<EMAIL_ADDRESS>Division of Particle and Astrophysical Science, Graduate School of Science,
Nagoya University, Nagoya 464-8602, Japan

###### Abstract

We investigate solutions of the $5$–dimensional rotating Einstein-Vlasov
system with an $R\times SU(2)\times U(1)$ isometry group. In a five-
dimensional spacetime, there are two independent planes of rotation, thus,
considering $U(1)$ symmetry on each rotation plane, we may impose an $R\times
U(1)\times U(1)$ isometry to a stationary spacetime. Furthermore, when the
values of the two angular momenta are equal to each other, the spatial
symmetry gets enhanced to $R\times SU(2)\times U(1)$ symmetry, and the
spacetime has a cohomogeneity-1 structure. Imposing the same symmetry to the
distribution function of the particles of which the Vlasov system consists,
the distribution function can be dependent on three mutually independent and
commutative conserved charges for particle motion (energy, total angular
momentum on $SU(2)$ and $U(1)$ angular momentum). We consider the distribution
function which exponentially depends on the $U(1)$ angular momentum and
reduces to the thermal equilibrium state in spherical symmetry. Then, in this
paper, we numerically construct solutions of the asymptotically AdS
Einstein–Vlasov system.

##  1 Introduction

General relativistic self-gravitating collisionless many-particle systems,
which are often called Einstein-Vlasov systems, have long been investigated in
astrophysics. In particular, spherically symmetric systems have been studied
in detail, and a great deal of research has been done on the existence of
solutions [1, 2, 3, 4] and the stability of the systems [5, 6, 7, 8, 9, 10]
(see also a review [11]). To realize a specific configuration of the Einstein-
Vlasov system, it is necessary to impose an appropriate ansatz on the
distribution function. For example, to give a static configuration surrounding
a black hole, we have to assume a distribution with a lower cutoff of the
angular momentum [12]. If we assume a Maxwell–Jüttner distribution, we can
realize a thermal equilibrium state of the self-gravitating system. However,
there are no thermal equilibrium states with finite mass in asymptotically
flat spacetimes because gravity is a long-range interaction. Antonov
investigated the non-relativistic many-particle systems with finite mass by
introducing an adiabatic wall confining the particles [13], and the properties
of those systems have been generalized and investigated in detail by Lynden-
Bell and Wood [14].

The thermodynamical instability of self-gravitating systems is often called
the gravothermal catastrophe which also applies to relativistic cases. It
should be noted that, however, since the thermal equilibrium states have
infinite mass with a vanishing cosmological constant, the analyses of the
gravothermal catastrophe in asymptotically flat spacetime always require an
artificial wall to confine the system. On the other hand, in a system with a
negative cosmological constant, the AdS barrier confines the particle system,
allowing it to naturally be in a thermal equilibrium state. Some of the
authors have constructed a thermal equilibrium state of such a system confined
by the AdS barrier under the assumption of static spherical symmetry and
analyzed its stability[15, 16].

Properties of the Einstein-Vlasov system in static and spherically symmetric
cases have been intensively studied, but there are still few studies on
systems with rotation[17, 18, 19, 20, 21]. The properties of a system, such as
stability, generally depend on the presence of angular momentum because the
total angular momentum may prevent the system from collapsing. In the case of
a self-gravitating many-particle system, it is expected that the instability
associated with gravitational collapse may be inhibited by the angular
momentum of the system. However, the existence of the non-zero angular
momentum inevitably reduces the spacetime symmetry, and the analyses become
much more difficult. To avoid this technical difficulty, we focus on
$5$–dimensional spacetimes because it is known that the spacetime can have a
cohomogeneity-1 structure even with non-zero angular momentum in
$5$–dimensional spacetimes. That is, we can construct a spacetime with finite
angular momentum solving a set of ordinary differential equations for unknown
variables depending only on a radial coordinate. More specifically, in a
$5$–dimensional spacetime, if the values of the angular momenta on the two
independent rotation planes are equal to each other, the spacetime symmetry
can be enhanced to $R\times SU(2)\times U(1)$. The corresponding black hole
solution is called the Myers-Perry (AdS) black hole with equal angular
momenta[22, 23, 24, 25, 26]. Due to its high symmetry, this spacetime is often
used in the analyses of gravitational perturbations and phenomena specific to
rotating systems.

In this paper, we consider rotating solutions for the Einstein-Vlasov system
with an $R\times SU(2)\times U(1)$ isometry group appropriately setting the
distribution function of the Vlasov field. A negative cosmological constant is
introduced for the realization of the stationary solutions with finite total
mass and angular momentum. In general, a resultant spacetime has an
asymptotically locally AdS structure, whose spatial geometry is given by the
foliation of squashed $S^{3}$ hyper-surfaces even at spatial infinity. The
squashing parameter at spatial infinity can be set to zero by tuning the
boundary condition at the center, and the spacetime can be asymptotically AdS
without squashing of $S^{3}$ at spatial infinity. We note that a similar
situation has been reported in the vacuum cases [27].

The purpose of introducing a negative cosmological constant is not only for
the construction of a physical solution with finite mass. Recently,
asymptotically AdS spacetimes have attracted much attention in the context of
the AdS/CFT correspondence [28, 29, 30] and the gravitational turbulent
instability [31]. The conditions for the onset of the instability have not yet
been clarified [32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45], and
there are still few clues to the final state. The system provided in this
paper may be treated as a macroscopic model of the final state of a system
complicated by turbulent phenomena, and we expect that our analyses will be
helpful to get useful insights into the final state of turbulent instability.
In addition, the instability in the Einstein-Vlasov system is also discussed
in Refs. [46, 8, 47], and the possible relation with the Hawking-Page
transition has been pointed out in Refs. [48, 49]. Another related phenomenon
is the superradiant instability of rotating black holes in asymptotically AdS
spacetime, for which the existence of finite angular momentum is essential.
The final fate of the superradiant instability has not been also clarified
yet. In order to approach the superradiant instability through the
construction of a macroscopic model with the Einstein-Vlasov system, the
introduction of finite angular momentum is a necessary step to be performed.

This paper is organized as follows. In section 2, we provide the metric ansatz
of the spacetime with an $R\times SU(2)\times U(1)$ isometry group. We also
define the conserved quantities along the geodesic of a particle and list the
conditions which we impose on the metric functions for technical and practical
reasons. In section 3, introducing a specific form of the distribution
function, we show the explicit forms for the energy-momentum tensors (detailed
calculations are shown in App. A). We write down the Einstein field equations
for our system in section 4 and numerically solve them in section 5. Section 6
is devoted to a summary and conclusion.

Throughout this paper, we use the geometrized units in which both the speed of
light and gravitational constant in $5$–dimension are unity, $c=G=1$.

##  2 Metric Ansatz and conserved quantities

###  2.1 Metric ansatz with an $R\times SU(2)\times U(1)$ isometry group

We start with the following form of the metric:

$\displaystyle\bm{g}\coloneqq-{\mathrm{e}^{2{{\mu}\pqty{r}}}}\differential{t}^{2}+{\mathrm{e}^{2{{\nu}\pqty{r}}}}\differential{r}^{2}+\frac{r^{2}}{4}\bqty{\pqty{\bm{\sigma}^{1}}^{2}+\pqty{\bm{\sigma}^{2}}^{2}+\pqty{\bm{\sigma}^{3}}^{2}}+{{h}\pqty{r}}\pqty{\differential{t}-\frac{{{a}\pqty{r}}}{2}\bm{\sigma}^{3}}^{2},$
(2.1)

where $\bm{\sigma}^{i}$ are one-forms defined as

$\displaystyle\bm{\sigma}^{1}$
$\displaystyle\coloneqq-\sin{\phi}\differential{\theta}+\sin{\theta}\cos{\phi}\differential{\psi},$
(2.2a) $\displaystyle\bm{\sigma}^{2}$
$\displaystyle\coloneqq-\cos{\phi}\differential{\theta}-\sin{\theta}\sin{\phi}\differential{\psi},$
(2.2b) $\displaystyle\bm{\sigma}^{3}$
$\displaystyle\coloneqq\differential{\phi}+\cos{\theta}\differential{\psi},$
(2.2c)

satisfying the Maurer-Cartan equation
$\differential{\bm{\sigma}^{i}}+\tfrac{1}{2}{\epsilon^{i}}_{jk}\bm{\sigma}^{j}\wedge\bm{\sigma}^{k}=0$.
The ranges of the coordinate variables are given by
$t\in\pqty{{-\infty},{\infty}}$, $r\in\left[{{0},{\infty}}\right)$,
$\theta\in\bqty{{0},{\pi}}$, $\phi\in\left[{{0},{4\pi}}\right)$ and
$\psi\in\left[{{0},{2\pi}}\right)$. The third term in the metric

$\displaystyle\bm{\gamma}\coloneqq$
$\displaystyle\frac{1}{4}\bqty{\pqty{\bm{\sigma}^{1}}^{2}+\pqty{\bm{\sigma}^{2}}^{2}+\pqty{\bm{\sigma}^{3}}^{2}}$
$\displaystyle=$
$\displaystyle\frac{1}{4}\pqty{\differential{\theta}^{2}+\differential{\phi}^{2}+\differential{\psi}^{2}+2\cos\theta\differential{\phi}\differential{\psi}}.$
(2.3)

describes the metric on $S^{3}$.

The three-sphere $S^{3}$ has two sets of $SU(2)$ generators
$\Bqty{\bm{\xi}_{i}}_{i\in\Bqty{1,2,3}}$ and
$\Bqty{\bm{\sigma}_{i}}_{i\in\Bqty{1,2,3}}$ written in the form:

$\displaystyle\bm{\xi}_{1}$
$\displaystyle=-\sin{\psi}{\partial_{\theta}}+\csc{\theta}\cos{\psi}{\partial_{\phi}}-\cot{\theta}\cos{\psi}{\partial_{\psi}},$
(2.4a) $\displaystyle\bm{\xi}_{2}$
$\displaystyle=-\cos{\psi}{\partial_{\theta}}-\csc{\theta}\sin{\psi}{\partial_{\phi}}+\cot{\theta}\sin{\psi}{\partial_{\psi}},$
(2.4b) $\displaystyle\bm{\xi}_{3}$ $\displaystyle={\partial_{\psi}},$ (2.4c)

and

$\displaystyle\bm{\sigma}_{1}$
$\displaystyle=-\sin{\phi}{\partial_{\theta}}-\cot{\theta}\cos{\phi}{\partial_{\phi}}+\csc{\theta}\cos{\phi}{\partial_{\psi}},$
(2.5a) $\displaystyle\bm{\sigma}_{2}$
$\displaystyle=-\cos{\phi}{\partial_{\theta}}+\cot{\theta}\sin{\phi}{\partial_{\phi}}-\csc{\theta}\sin{\phi}{\partial_{\psi}},$
(2.5b) $\displaystyle\bm{\sigma}_{3}$ $\displaystyle={\partial_{\phi}}.$
(2.5c)

The $SU(2)$ generators $\Bqty{\bm{\xi}_{i},\bm{\sigma}_{i}}$ satisfy the
following commutation relations:

$\displaystyle\bqty{\bm{\xi}_{i},\bm{\xi}_{j}}={\epsilon_{ij}}^{k}\bm{\xi}_{k},\quad\bqty{\bm{\sigma}_{i},\bm{\sigma}_{j}}={\epsilon_{ij}}^{k}\bm{\sigma}_{k},\quad\bqty{\bm{\xi}_{i},\bm{\sigma}_{j}}=\bqty{\bm{\sigma}_{i},\bm{\xi}_{j}}=0,$
(2.6)

or equivalently

$\displaystyle{\mathcal{L}_{\bm{\xi}_{i}}}\bm{\xi}^{j}={{\epsilon_{i}}^{j}}_{k}\bm{\xi}^{k},\quad{\mathcal{L}_{\bm{\sigma}_{i}}}\bm{\sigma}^{j}={{\epsilon_{i}}^{j}}_{k}\bm{\sigma}^{k},\quad{\mathcal{L}_{\bm{\xi}_{i}}}\bm{\sigma}^{j}=0.$
(2.7)

The Killing vectors $\Bqty{\bm{\xi}_{i}}$ and $\bm{\sigma}_{3}$ on the three-
sphere are also the Killing vectors on the spacetime due to Eq. (2.7). On the
other hand, neither $\bm{\sigma}_{1}$ nor $\bm{\sigma}_{2}$ is the Killing
vector on the spacetime due to the last term in Eq. (2.1) unless
${{a}\pqty{r}}=0$ everywhere. Since the vector $\bm{\sigma}_{3}$ generates a
$U(1)$ isometry group and the metric has a time-like Killing vector
$\bm{\eta}={\partial_{t}}$, the spacetime has the $R_{t}\times
SU(2)_{\bm{\xi}}\times U(1)_{\bm{\sigma}}$ isometry group. The black hole
solutions which have the same symmetry are known as the Myers–Perry black
holes with equal angular momenta [24, 25]. It would be worth noting that the
angular coordinates $\theta$, $\phi$ and $\psi$ are related to the Hopf
coordinates $\tilde{\theta}$, $\tilde{\phi}$ and $\tilde{\psi}$ through
$\theta=2\tilde{\theta}$, $\psi=-\tilde{\phi}+\tilde{\psi}$ and
$\phi=\tilde{\phi}+\tilde{\psi}$. Then the two equal angular momenta are the
conserved charges associated with the Killing vectors
$\partial_{\tilde{\phi}}$ and $\partial_{\tilde{\psi}}$.

In this paper, to avoid possible technical problems, we focus on the cases
satisfying the following four conditions:

  1. 1.

Non-degeneracy: $\det\bm{g}<0$,

  2. 2.

No-horizon: $\bm{r}\cdot\bm{r}>0$,

  3. 3.

Time-like Killing vector exists everywhere: $\bm{\eta}\cdot\bm{\eta}<0$,

  4. 4.

Time-like unit normal exists everywhere: $\bm{n}\cdot\bm{n}<0$,

where $\bm{r}\coloneqq\absolutevalue{g^{rr}}^{-\frac{1}{2}}\differential{r}$
and $\bm{n}\coloneqq-\absolutevalue{g^{tt}}^{-\frac{1}{2}}\differential{t}$
are the unit normal to the $r={\rm const.}$ and the $t={\rm const.}$ hyper-
surfaces, respectively. The condition 3 implies that the spacetime has no
ergo-region555In Ref. [20], the solution with an ergo-region is constructed in
a 4–dimensional asymptotically flat spacetime. . The conditions 1 \- 4 yield

$\displaystyle{\mathrm{e}^{2\mu(r)}}-h(r),\ {\mathrm{e}^{2\nu(r)}},\ F(r),\
G(r)>0,$ (2.8)

for any $r>0$, where ${{G}\pqty{r}}\coloneqq r^{2}+a^{2}h$ and
${{F}\pqty{r}}\coloneqq{\mathrm{e}^{2\mu}}G-hr^{2}$.

###  2.2 Conserved quantities for the geodesic motion

The symmetry of the spacetime indicates the existence of conserved quantities
for the geodesic motion in the form of the inner product between the Killing
vector and the momentum of a particle. Then we can consider the following
mutually independent conserved quantities for the metric (2.1):

$\displaystyle\varepsilon$ $\displaystyle=-\bm{p}\cdot\bm{\eta}=-p_{t},$
(2.9a) $\displaystyle{J_{\xi_{1}}}$
$\displaystyle=\bm{p}\cdot\bm{\xi}_{1}=-p_{\theta}\sin{\psi}+p_{\phi}\csc{\theta}\cos{\psi}-p_{\psi}\cot{\theta}\cos{\psi},$
(2.9b) $\displaystyle{J_{\xi_{2}}}$
$\displaystyle=\bm{p}\cdot\bm{\xi}_{2}=-p_{\theta}\cos{\psi}-p_{\phi}\csc{\theta}\sin{\psi}+p_{\psi}\cot{\theta}\sin{\psi},$
(2.9c) $\displaystyle{J_{\xi_{3}}}$
$\displaystyle=\bm{p}\cdot\bm{\xi}_{3}=p_{\psi},$ (2.9d)
$\displaystyle{j_{\sigma}}$
$\displaystyle=\bm{p}\cdot\bm{\sigma}_{3}=p_{\phi},$ (2.9e)

where $\bm{p}$ is the momentum of the particle and ‘$\cdot$’ denotes the inner
product concerning $\bm{g}$. The total angular momentum $J_{\xi}\geq 0$ of the
$SU(2)_{\bm{\xi}}$ sector can be defined as

$\displaystyle{J_{\xi}}^{2}\coloneqq\sum_{i=1}^{3}{J_{\xi_{i}}}^{2}=\frac{1}{4}\gamma^{\mu\nu}p_{\mu}p_{\nu},$
(2.10)

where we have defined $\gamma^{\mu\nu}$ as $\gamma^{\mu\nu}\coloneqq
4\sum_{i}{\sigma_{i}}^{\mu}{\sigma_{i}}^{\nu}$. Then $\varepsilon$,
$j_{\sigma}$ and $J_{\xi}$ are mutually commutative conserved charges and can
be used as independent coordinates in the phase space.

The momentum of the particle must satisfy the on-shell condition:
$p^{2}+m^{2}=0$, which can be rewritten as follows:

$\displaystyle\frac{{\mathrm{e}^{-2\nu}}G}{F}\pqty{\varepsilon-\frac{2ah}{G}j_{\sigma}}^{2}-{\mathrm{e}^{-2\nu}}\bqty{m^{2}+\frac{4}{r^{2}}\pqty{{J_{\xi}}^{2}-\frac{a^{2}h}{G}{j_{\sigma}}^{2}}}=\pqty{p^{r}}^{2}.$
(2.11)

We note that the left-hand side of Eq. (2.11) can be regarded as the effective
potential for the geodesic motion. Here we impose that the momentum is future
pointing $\bm{n}\cdot\bm{p}<0$. Then the positivity of the local energy of the
particle is ensured:

$\displaystyle\varepsilon-\frac{2ah}{G}j_{\sigma}>0.$ (2.12)

The allowed region in the momentum space for the particle is the subspace
satisfying the conditions (2.11) and (2.12).

##  3 Model and physical quantities

###  3.1 Distribution function and the energy-momentum tensor

Considering a collisionless many-particle system, the particles follow the
geodesic motion. Therefore the distribution function $f$ satisfies the Vlasov
(collisionless Boltzmann) equation:

$\displaystyle
p^{\mu}{\nabla_{\mu}}f=p^{\mu}\partialderivative{f}{x^{\mu}}-{\Gamma^{i}}_{\mu\nu}p^{\mu}p^{\nu}\partialderivative{f}{p^{i}}=0,$
(3.1)

which implies the conservation of the distribution function along the
geodesic. If the distribution function is written in the form

$\displaystyle f=f(\varepsilon,J_{\xi},j_{\sigma}),$ (3.2)

the Vlasov equation (3.1) is automatically satisfied because $\varepsilon$,
$J_{\xi}$ and $j_{\sigma}$ are conserved quantities along the geodesic.

In this paper, as a simple specific model, we assume the following
distribution function

$\displaystyle{{f}\pqty{\varepsilon,j_{\sigma}}}=\exp\bqty{\alpha-\beta(\varepsilon-\Omega
j_{\sigma})},$ (3.3)

with constants $\alpha\in\mathbb{R}$, $\beta>0$ and $\Omega>0$. The
distribution function (3.3) reduces to the Maxwell–Jüttner distribution, which
describes the relativistic thermal equilibrium states in static cases, with
$\Omega=0$. Thus the system with this distribution function can be regarded as
an extension of the thermal equilibrium state to the system with a finite
angular momentum. The Maxwell–Jüttner distribution is derived by extremizing
the entropy of the system fixing the total mass and the total particle number
in static cases. On the other hand, in rotating cases, it is not clear how to
determine the most probable rotating state. Therefore the functional form
(3.3) should be regarded as a working assumption for $\Omega\neq 0$.

Given a one-particle distribution function $f$, the energy-momentum tensor
$T_{\mu\nu}$ is given by

$\displaystyle
T_{\mu\nu}\coloneqq\int\differential{V_{p}}p_{\mu}p_{\nu}f(x^{\mu},p^{i}),$
(3.4)

where $\differential{V_{p}}$ is the integral measure in the momentum space:

$\displaystyle\differential{V_{p}}\coloneqq-\frac{2{{\delta}\pqty{p^{2}+m^{2}}}{{\Theta}\pqty{\varepsilon-2ahj_{\sigma}/G}}}{\sqrt{-\det\bm{g}}}\differential{p_{t}}\wedge\differential{p_{r}}\wedge\differential{p_{\theta}}\wedge\differential{p_{\phi}}\wedge\differential{p_{\psi}}.$
(3.5)

Here the delta function $\delta$ and the Heaviside’s step function $\Theta$
describe the conditions (2.11) and (2.12), respectively.

To obtain numerical solutions, we have to know the local expressions of
energy-momentum tensor by integrating over the momentum space. In spherical
cases, we can perform the integration analytically and obtain explicit
expressions because the metric is diagonal and the on-shell condition is
simple as shown in Refs. [15, 16]. On the other hand, in rotating cases, we
cannot perform analytical integrations due to the off-diagonal components of
the metric and the cross term in Eq. (2.11). As is shown in App. A, however,
we can reduce it to one-dimensional integrals of the normalized energy under
the ansatz (3.3). That is, we can obtain the values of the energy-momentum
tensor at each radial coordinate $r$ by performing one-dimensional numerical
integration. The explicit expressions are shown in App. A.

###  3.2 Komar integral and total particle number

Before solving the Einstein equations, let us define global conserved
quantities characterizing the system in terms of the Komar integrals.

We define the total mass of the system by the Komar mass:

$\displaystyle
M_{K}\coloneqq\int\differential{\Sigma}n_{\mu}{T^{\mu}}_{\nu}\eta^{\nu},$
(3.6)

where
$\differential{\Sigma}=2\pi^{2}r^{2}{\mathrm{e}^{\nu}}\sqrt{G}\differential{r}$
is the invariant volume element on a $t={\rm const.}$ hyper-surface. In our
formulation, Eq. (3.6) becomes

$\displaystyle
M_{K}=2\pi^{2}\int\differential{r}\frac{r^{2}{\mathrm{e}^{\nu}}}{\sqrt{F}}(GT_{tt}+2ahT_{t\phi}).$
(3.7)

We note that the Komar mass for the distribution (3.6) can be finite due to
the potential wall for massive particles associated with a negative
cosmological constant $\Lambda<0$.

We also define the total angular momentum of the system by

$\displaystyle
J_{\phi}\coloneqq-\int\differential{\Sigma}n_{\mu}{T^{\mu}}_{\nu}{\sigma_{3}}^{\nu}.$
(3.8)

Similarly to the Komar mass, in our formulation, we obtain the expression:

$\displaystyle
J_{\phi}=-2\pi^{2}\int\differential{r}\frac{r^{2}{\mathrm{e}^{\nu}}}{\sqrt{F}}(GT_{t\phi}+2ahT_{\phi\phi}).$
(3.9)

Since the Komar integrals associated with the Killing vectors conserve, we can
choose them as quantities characterizing the system.

##  4 Einstein equations

###  4.1 Field equations

For a 5–dimensional spacetime, the Einstein field equations
$G_{\mu\nu}+\Lambda g_{\mu\nu}=8\pi T_{\mu\nu}$ with $\Lambda<0$ can be
rewritten as

$\displaystyle
R_{\mu\nu}=8\pi\pqty{{T}_{\mu\nu}-\frac{{T}}{3}g_{\mu\nu}-\frac{g_{\mu\nu}}{2\pi
L^{2}}},$ (4.1)

where $L\coloneqq\sqrt{-6/\Lambda}$ is the AdS radius.

For numerical analyses, we rewrite the equations by dimensionless variables.
We define $\ell\coloneqq(8\pi T_{tt}(0))^{-1/2}$ and consider normalized
variables $x\coloneqq r/\ell$ and $\lambda\coloneqq L/\ell$. As alternative
functions, we introduce the following dimensionless metric functions:

$\displaystyle y_{1}(x)\coloneqq\frac{ah}{\ell},\quad
y_{2}(x)\coloneqq\frac{a^{2}h}{\ell^{2}},\quad
y_{3}(x)\coloneqq{\mathrm{e}^{2\mu}}-h,\quad
y_{4}(x)\coloneqq{\mathrm{e}^{2\nu}}.$ (4.2)

We also define $f_{1}(x)\coloneqq G/\ell^{2}=x^{2}+y_{2}$ and
$f_{2}(x)\coloneqq F/\ell^{2}={y_{1}}^{2}+f_{1}y_{3}$ for a notational
simplicity.

Defining

$\displaystyle
s_{\mu\nu}\coloneqq\frac{1}{T_{tt}(0)}\bqty{{T}_{\mu\nu}-\frac{T}{3}g_{\mu\nu}-\frac{1}{2\pi
L^{2}}g_{\mu\nu}},$ (4.3)

relevant components of the Einstein equations are as follows:

$\displaystyle{y_{1}}^{\prime\prime}$
$\displaystyle=\frac{2x({{{y_{1}}^{\prime}}}{y_{3}}-{{{y_{3}}^{\prime}}}{y_{1}})}{f_{2}}+\frac{{f_{1}}{{{y_{1}}^{\prime}}}{{{y_{3}}^{\prime}}}+{{{y_{1}}^{\prime}}}^{2}{y_{1}}+{{{y_{1}}^{\prime}}}{{{y_{2}}^{\prime}}}{y_{3}}-{{{y_{2}}^{\prime}}}{{{y_{3}}^{\prime}}}{y_{1}}}{f_{2}}$
$\displaystyle\hskip
20.0pt+4s_{t\phi}{y_{4}}+\frac{{{{y_{1}}^{\prime}}}({f_{3}}{y_{4}}-2)}{2x}+\frac{2{{{y_{1}}^{\prime}}}{y_{2}}{y_{4}}}{x^{3}}+\frac{4{f_{1}}{y_{1}}{y_{4}}}{x^{4}},$
(4.4a) $\displaystyle{y_{2}}^{\prime\prime}$
$\displaystyle=\frac{4x({{{y_{1}}^{\prime}}}{y_{1}}+{{{y_{2}}^{\prime}}}{y_{3}})}{f_{2}}+\frac{{{{y_{2}}^{\prime}}}^{2}{y_{3}}-{f_{1}}{{{y_{1}}^{\prime}}}^{2}+2{{{y_{1}}^{\prime}}}{{{y_{2}}^{\prime}}}{y_{1}}-4({y_{1}}^{2}+{y_{2}}{y_{3}})}{f_{2}}$
$\displaystyle\hskip
20.0pt+{f_{4}}{y_{4}}+\frac{{{{y_{2}}^{\prime}}}({f_{3}}{y_{4}}-2)}{2x}+\frac{2{y_{2}}{y_{4}}}{x^{2}}\pqty{6+\frac{{{y_{2}}^{\prime}}}{x}+\frac{2{y_{2}}}{x^{2}}},$
(4.4b) $\displaystyle{y_{3}}^{\prime\prime}$
$\displaystyle=\frac{{f_{1}}{{{y_{3}}^{\prime}}}^{2}-{{{y_{1}}^{\prime}}}^{2}{y_{3}}+2{{{y_{1}}^{\prime}}}{{{y_{3}}^{\prime}}}{y_{1}}}{f_{2}}$
$\displaystyle\hskip
20.0pt+2s_{tt}{y_{4}}+\frac{{{{y_{3}}^{\prime}}}({f_{3}}{y_{4}}-2)}{2x}+\frac{2{{{y_{3}}^{\prime}}}{y_{2}}{y_{4}}}{x^{3}}-\frac{4{y_{1}}^{2}{y_{4}}}{x^{4}},$
(4.4c) $\displaystyle{y_{4}}^{\prime}$
$\displaystyle={y_{4}}\bqty{\frac{2x{y_{3}}}{f_{2}}+\frac{{f_{1}}{{{y_{3}}^{\prime}}}+2{{{y_{1}}^{\prime}}}{y_{1}}+{{{y_{2}}^{\prime}}}{y_{3}}}{f_{2}}+\frac{{f_{3}}{y_{4}}+2}{x}+\frac{4{y_{2}}{y_{4}}}{x^{3}}},$
(4.4d)

where $\ {}^{\prime}\coloneqq\derivative*{x}$,
$f_{3}\coloneqq-4(1+s_{\phi\phi})+\gamma^{\mu\nu}s_{\mu\nu}$ and
$f_{4}\coloneqq-12s_{\phi\phi}+\gamma^{\mu\nu}s_{\mu\nu}$. We can show the
other components of the Einstein equations are redundant by using Eqs. (4.4)
and continuity equations.

###  4.2 Asymptotic behaviors and boundary conditions

Let us consider the boundary conditions for Eqs. (4.4) to realize
asymptotically AdS spacetimes. We can rewrite the metric as

$\displaystyle\ell^{-2}\bm{g}=-\frac{f_{2}}{f_{1}}\differential{t}^{2}+{y_{4}}\differential{x}^{2}+\frac{x^{2}}{4}\bqty{\pqty{\bm{\sigma}^{1}}^{2}+\pqty{\bm{\sigma}^{2}}^{2}+\frac{f_{1}}{x^{2}}\pqty{\bm{\sigma}^{3}-\frac{2y_{1}}{f_{1}}\differential{t}}^{2}}.$
(4.5)

Setting $T_{\mu\nu}=0$ and expanding the metric functions in the vicinity of
the boundary, we obtain the following asymptotic solutions:

$\displaystyle\frac{f_{2}}{f_{1}}$
$\displaystyle=\frac{x^{2}}{\lambda^{2}}+1-s+\order{x^{-1}},$ (4.6a)
$\displaystyle\frac{1}{y_{4}}$
$\displaystyle=\frac{x^{2}}{\lambda^{2}}+1-\frac{3}{5}s+\order{x^{-1}},$
(4.6b) $\displaystyle\frac{f_{1}}{x^{2}}$
$\displaystyle=1+s-2s\pqty{1+s}\frac{\lambda^{2}}{x^{2}}+\order{x^{-3}},$
(4.6c) $\displaystyle\frac{2y_{1}}{f_{1}}$
$\displaystyle=\omega+\order{x^{-3}},$ (4.6d)

where $s$ and $\omega$ can be regarded as free parameters for the asymptotic
solution, and the coordinates $t$ and $r$ are fixed so that the pure AdS
metric can be realized for $s=\omega=0$. The leading term $\omega$ of
$2y_{1}/f_{1}$ can be eliminated by employing the co-rotating coordinate
$\hat{\phi}\coloneqq\phi-\omega t$. Then we can introduce one-forms
$\hat{\sigma}^{i}$ at the boundary by replacing $\phi$ with $\hat{\phi}$.
Using these one-forms, we can write

$\displaystyle\ell^{-2}\bm{g}=-\frac{f_{2}}{f_{1}}\differential{t}^{2}+{y_{4}}\differential{x}^{2}+\frac{x^{2}}{4}\bqty{\pqty{\hat{\bm{\sigma}}^{1}}^{2}+\pqty{\hat{\bm{\sigma}}^{2}}^{2}+\frac{f_{1}}{x^{2}}\pqty{\hat{\bm{\sigma}}^{3}}^{2}}.$
(4.7)

As we can see from Eq. (4.6) and Eq. (4.7), the parameter $s$ should vanish so
that the spacetime has the asymptotically AdS structure. We note that the
contribution of the finite value of $s$ corresponds to the squashing of
$S^{3}$ at infinity as is reported in Ref. [27] for vacuum cases.

Next, let us investigate the asymptotic behavior around the center by
expanding the metric functions as

$\displaystyle y_{1}(x)$ $\displaystyle\simeq
x^{2}\bqty{{z_{1}}^{(0)}+{z_{1}}^{(1)}x+\frac{{z_{1}}^{(2)}}{2}x^{2}+\order{x^{3}}},$
(4.8a) $\displaystyle y_{2}(x)$ $\displaystyle\simeq
x^{2}\bqty{{z_{2}}^{(0)}+{z_{2}}^{(1)}x+\frac{{z_{2}}^{(2)}}{2}x^{2}+\order{x^{3}}},$
(4.8b) $\displaystyle y_{3}(x)$
$\displaystyle\simeq{y_{3}}^{(0)}+{y_{3}}^{(1)}x+\frac{{y_{3}}^{(2)}}{2}x^{2}+\order{x^{3}},$
(4.8c) $\displaystyle y_{4}(x)$
$\displaystyle\simeq{y_{4}}^{(0)}+{y_{4}}^{(1)}x+\frac{{y_{4}}^{(2)}}{2}x^{2}+\order{x^{3}},$
(4.8d)

with coefficients ${z_{1}}^{(i)}$, ${z_{2}}^{(i)}$, ${y_{3}}^{(i)}$ and
${y_{4}}^{(i)}$. As regularity conditions, first, we assume the metric
components for the Cartesian coordinates $(t,x_{1},x_{2},x_{3},x_{4})$666 The
spatial Cartesian coordinates can be defined by  $\displaystyle\begin{aligned}
x_{1}&=r\sin\left(\frac{\theta}{2}\right)\cos\left[\frac{1}{2}\left(\phi-\psi\right)\right],\quad
x_{2}=r\sin\left(\frac{\theta}{2}\right)\sin\left[\frac{1}{2}\left(\phi-\psi\right)\right],\\\
x_{3}&=r\cos\left(\frac{\theta}{2}\right)\cos\left[\frac{1}{2}\left(\phi+\psi\right)\right],\quad
x_{4}=r\cos\left(\frac{\theta}{2}\right)\sin\left[\frac{1}{2}\left(\phi+\psi\right)\right].\end{aligned}$
(4.9) to be $g_{tt}=-y_{3}^{(0)}$, $g_{x_{i}x_{i}}=1+\order{x^{2}}$, and
$\order{x^{2}}$ for other components. Then, we require all the Ricci tensor
components $R_{\mu\nu}$ to be $R_{\mu\nu}={\rm const.}+\order{x}$. As a
result, we obtain

$\displaystyle{z_{1}}^{(0)}={z_{1}}^{(1)}={z_{2}}^{(0)}={z_{2}}^{(1)}={y_{3}}^{(1)}=0,\quad{y_{4}}^{(0)}=1,$
(4.10)

and ${y_{3}}^{(0)}$ is kept as a free parameter. As is shown in App. B, the
system has another free parameter ${z_{2}}^{(2)}$.

When we solve the Einstein equations from the center, in general, we obtain
the solution with finite values of $s$ and $\omega$. As will be shown later,
we obtain a finite value of $s$ at infinity without fine-tuning for the
boundary condition at the center. In other words, we can set $s=0$ at infinity
by tuning ${z_{2}}^{(2)}$ appropriately as same as Ref. [27]. Therefore, we
can construct an asymptotically AdS Einstein-Vlasov system by tuning the
parameters $y_{3}^{(0)}$ and ${z_{2}}^{(2)}$. As a whole, the system has one
free parameter $y_{3}^{(0)}$ corresponding to the linear scaling of the time
coordinate, one tuning parameter ${z_{2}}^{(2)}$ and three physical parameters
$(\beta,\Omega,\lambda)$.

To obtain a solution, we solve Einstein equations (4.4) with boundary
conditions (4.10) at the center and a parameter set $(\beta,\Omega,\lambda)$.
However, for a finite value of $\Omega$, the integral (A.25) cannot be
performed analytically. Therefore we need to numerically integrate the
integral (A.25) at each step of solving the differential equations (4.4). In
our numerical simulations, we used 4th-order Runge-Kutta methods in numerical
integrations of Eq. (A.25) and Eq. (4.4).

##  5 Results

###  5.1 Static and spherically symmetric cases

First, we investigate the parameter dependence for the spherically symmetric
cases with ${z_{2}}^{(2)}=0$, which corresponds to $s=0$. Fig. 1 shows the
normalized energy density profile for the matter sector
$\rho(x)\coloneqq-{T^{t}}_{t}/{T}_{tt}(0)$ as a function of $x$ for
$\Omega=0$.

(a) $\beta$ dependence for $\lambda=1$

(b) $\lambda$ dependence for $\beta=1$

Fig. 1:  The parameter dependence of the distribution of the energy density
for the spherically symmetric cases. As the parameter $\beta$ increases, the
gravitational potential gets deeper and the system gets more compact. As the
parameter $\lambda$ decreases, the system gets more compact because the value
of $\lambda$ characterizes the radius of the potential wall associated with
the negative cosmological constant.

In the spherically symmetric cases, the system is characterized by the two
parameters $(\beta,\lambda)$. The parameter $\beta$ corresponds to the inverse
temperature of the system and it characterizes the depth of the gravitational
potential. Thus, increasing the value of $\beta$ for a fixed value of
$\lambda$, the density distribution gets sharper and the value of the total
mass increases as is shown in Fig. 1(a) and Fig. 2, respectively. Fig. 2 shows
the quasi-local mass given by the Komar integral as a function of $x$. For all
the cases, the mass increases as $M_{K}\propto x^{4}$ for $x\ll 1$ and
asymptotes to a constant value for $x\gg\lambda$. Since the AdS radius
$\lambda$ gives the characteristic radius of the potential wall which confines
the Vlasov matter, the Komar mass gets larger for large $\lambda$ if the
parameter $\beta$ is fixed.

(a) $\beta$ dependence for $\lambda=1$

(b) $\lambda$ dependence for $\beta=1$

Fig. 2:  The parameter dependence of the Komar mass for the spherically
symmetric cases. The system with relatively small $\beta$ (i.e., high-
temperature system) has a large total mass, which implies a deep gravitational
potential wall. Since the parameter $\lambda$ works like the size of the
system, the Komar mass takes a smaller value for smaller $\lambda$.

###  5.2 Stationary and rotating case

We investigate the non-spherically symmetric cases by setting finite values of
$\Omega$. Since the dependence on $\beta$ and $\lambda$ is similar to the
spherically symmetric cases, we focus on the $\Omega$ dependence and the
distribution of the angular momentum density of the system. Fig. 3 shows the
profile of the angular momentum density $j_{\phi}\coloneqq{T^{t}}_{\phi}$ and
the Komar angular momentum $J_{\phi}$ for ${z_{2}}^{(2)}=0$, $\beta=0.01$ and
$\lambda=10$. For $\Omega\neq 0$, the spherical symmetry is broken and the
system has a finite angular momentum. Since the angular momentum is carried by
the Vlasov matter, the distribution of the angular momentum density also
starts to decay around the AdS radius $x=\lambda$. Therefore, as is shown in
Fig. 3, the angular momentum density increases as $x$ increases up to a
certain value and then it decreases and vanishes asymptotically.

(a) Angular momentum density.

(b) Komar angular momentum.

Fig. 3:  The angular momentum density $j_{\phi}={T^{t}}_{\phi}$ and the total
angular momentum in the Komar expression for $\beta=0.02$ and $\lambda=10$.

Let us check the functions $f_{1}/x^{2}$ and $2y_{1}/f_{1}$ which characterize
the squashing of $S^{3}$ and the frame-dragging effect, respectively. Fig. 4
shows the functions $f_{1}/x^{2}$ and $2y_{1}/f_{1}$ for ${z_{2}}^{(2)}=0$,
$\beta=0.001$ and $\lambda=10$.

(a) The function $f_{1}/x^{2}$

(b) The function $2y_{1}/f_{1}$

Fig. 4:  The functions $f_{1}/x^{2}$ and $2y_{1}/f1$ for ${z_{2}}^{(2)}=0$,
$\beta=0.001$ and $\lambda=10$. Both of the functions asymptote to constants.

Unlike the $\Omega=0$ case, $S^{3}$ is squashed at infinity even if we set
${z_{2}}^{(2)}=0$. This is because of the contribution of the rotating Vlasov
matter. Fig. 5 shows the ${z_{2}}^{(2)}$ dependence of the functions
$f_{1}/x^{2}$ and $2y_{1}/f_{1}$ for $\Omega=0.1$, $\beta=0.001$ and
$\lambda=10$.

(a) The function $f_{1}/x^{2}$

(b) The function $2y_{1}/f_{1}$

Fig. 5:  ${z_{2}}^{(2)}$ dependence of the functions $f_{1}/x^{2}$ and
$2y_{1}/f_{1}1$ with fixed $\Omega=0.1$, $\beta=0.001$ and $\lambda=10$. We
can construct solutions with $f_{1}/x^{2}=0$ at infinity keeping
$2y_{1}/f_{1}$ finite by tuning ${z_{2}}^{(2)}$.

As we can see from Fig. 5(a), we can choose the value of $f_{1}/x^{2}$ at
infinity by tuning ${z_{2}}^{(2)}$. On the other hand, the value of
$2y_{1}/f_{1}$ at infinity is kept finite in general although it can be
eliminated at infinity by taking the co-rotating frame. As in this example, we
can always construct an asymptotically AdS solution that is not squashed at
infinity with a finite angular momentum by tuning ${z_{2}}^{(2)}$.

##  6 Conclusion

We have constructed rotating Einstein-Vlasov systems with the $R_{t}\times
SU(2)_{\bm{\xi}}\times U(1)_{\bm{\sigma}}$ isometry group imposing the
exponential form of the distribution function as
$\sim\exp\left[\beta(\varepsilon-\Omega j_{\sigma})\right]$ with $\varepsilon$
and $j_{\sigma}$ being the conserved energy and angular momentum for the
particle motion. Under this assumption, the system reduces to the thermal
equilibrium state in the limit $\Omega\to 0$. Since we assume the cosmological
constant is negative, we can construct stationary solutions with finite mass
without any artificial wall due to the confined structure of the AdS
potential. The distribution of the angular momentum density also decays
towards infinity because the angular momentum is carried by the Vlasov matter.
Therefore we can construct the system with finite mass and total angular
momentum by putting $\Omega\neq 0$ and solving the Einstein equations.

Even if we set $\Omega=0$, the system has an additional degree of freedom
describing the squashing of $S^{3}$ at infinity. We can control the squashing
parameter $s$ at infinity to vanish by tuning the free parameter
${z_{2}^{(2)}}$ for the boundary conditions at the center. For $\Omega=0$
cases, we can construct non-squashed solutions by setting ${z_{2}^{(2)}}=0$
because the spherical symmetry is trivially kept with ${z_{2}^{(2)}}=0$. On
the other hand, for $\Omega\neq 0$, each spatial hyper-surface of a given
value of radial coordinate $r$ is non-trivially squashed with
${z_{2}^{(2)}}=0$ due to the contribution of the non-spherical distribution of
the Vlasov matter. Therefore the parameter ${z_{2}}^{(2)}$ has to be adjusted
for the solution with $s=0$ and the spacetime to be asymptotically AdS with
finite total angular momentum. We have explicitly demonstrated that this fine-
tuning is possible for a specific parameter set.

In this paper, we have assumed that the distribution function of the particles
is dependent on only $(\varepsilon,j_{\sigma})$ as a simple model. The system
has, however, one more conserved quantity $J_{\xi}$ so we can construct
solutions with $J_{\xi}$ dependence keeping the system having the same
symmetry. If we assume an appropriate ansatz, we may construct rotating shell
solutions with finite total angular momentum with/without a central black
hole. Although it is hard to analyze such a system in a $4$–dimensional
spacetime without spherical symmetry, similar analyses would be relatively
easier for a $5$–dimensional spacetime with respecting the symmetry of
$R_{t}\times SU(2)_{\bm{\xi}}\times U(1)_{\bm{\sigma}}$.

## Acknowledgments

We would like to thank K. Murata for his comments and helpful discussions.
This work was supported in part by JST SPRING, Grant Number JPMJSP2125 (H.
A.), and in part by JSPS KAKENHI Grant Numbers JP19H01895 (CY), JP20H05850
(CY), JP20H05853 (CY). The author (H. A.) would like to take this opportunity
to thank the “Interdisciplinary Frontier Next-Generation Researcher Program of
the Tokai Higher Education and Research System.”

##  Appendix A Integration over the momentum space

In this appendix, we derive the expressions (A.27). In A.1, to obtain
expressions available in massless cases, we do not set the particle mass to
unity. On the other hand, in A.2, we set $m=1$ because we strict the
distribution function to Eq. (3.3).

###  A.1 Integral measure

First, let us rewrite the integral measure in the momentum space:

$\displaystyle\differential{V_{p}}=-\frac{16{{\delta}\pqty{p^{2}+m^{2}}}{{\Theta}\pqty{\varepsilon-2ahj_{\sigma}/G}}}{r^{2}\sqrt{F}{\mathrm{e}^{\nu}}\sin\theta}\differential{p_{t}}\wedge\differential{p_{r}}\wedge\differential{p_{\theta}}\wedge\differential{p_{\phi}}\wedge\differential{p_{\psi}},$
(A.1)

in terms of the conserved quantities. Since the definition of the total
angular momentum $J_{\xi}$ (2.10) can be rewritten as

$\displaystyle({J_{\xi}}^{2}-{j_{\sigma}}^{2})\sin^{2}\theta={\pqty{p_{\theta}\sin\theta}}^{2}+(p_{\psi}-p_{\phi}\cos\theta)^{2},$
(A.2)

the variables $(J_{\xi},j_{\sigma})$ satisfy
${J_{\xi}}^{2}\geq{j_{\sigma}}^{2}$. Introducing the variable $\chi$ defined
as

$\displaystyle R\cos\chi=p_{\theta}\sin\theta,\quad
R\sin\chi=p_{\psi}-p_{\phi}\cos\theta,$ (A.3)

with $R^{2}\coloneqq({J_{\xi}}^{2}-{j_{\sigma}}^{2})\sin^{2}\theta$, we obtain

$\displaystyle\differential{p_{\theta}}\wedge\differential{p_{\phi}}\wedge\differential{p_{\psi}}$
$\displaystyle=R\csc\theta\differential{R}\wedge\differential{j_{\sigma}}\wedge\differential{\chi}$
$\displaystyle=J_{\xi}\sin\theta\differential{J_{\xi}}\wedge\differential{j_{\sigma}}\wedge\differential{\chi}.$
(A.4)

Defining alternative quantities as

$\displaystyle\tilde{\varepsilon}\coloneqq\varepsilon-\frac{2ah}{G}j_{\sigma},\quad\tilde{J}^{2}\coloneqq{J_{\xi}}^{2}-\frac{a^{2}h}{G}{j_{\sigma}}^{2},$
(A.5)

we can simplify the on-shell condition (2.11) as

$\displaystyle\frac{G}{F}\tilde{\varepsilon}^{2}-\pqty{m^{2}+\frac{4\tilde{J}^{2}}{r^{2}}}=q^{2},$
(A.6)

and the positive energy condition as $\tilde{\varepsilon}\geq 0$. Using this
condition, we can eliminate the delta function. Solving Eq. (A.6) with respect
to $q\coloneqq{\mathrm{e}^{\nu}}p^{r}$, we obtain

$\displaystyle{{q}\pqty{\varepsilon,J_{\xi},j_{\sigma}}}=\bqty{\frac{G}{F}\pqty{\varepsilon-\frac{2ah}{G}j_{\sigma}}^{2}-\pqty{m^{2}+\frac{4}{r^{2}}\pqty{{J_{\xi}}^{2}-\frac{a^{2}h}{G}{j_{\sigma}}^{2}}}}^{\frac{1}{2}}$
(A.7)

for positive $q$. Using the fact that the integrand is an even function of $q$
due to the ansatz of the distribution function (3.2), we obtain

$\displaystyle-2\int{{\delta}\pqty{p^{2}+m^{2}}}\differential{p_{t}}\wedge\differential{p_{r}}=\frac{2{\mathrm{e}^{\nu}}\differential{\varepsilon}}{{{q}\pqty{\varepsilon,J_{\xi},j_{\sigma}}}}.$
(A.8)

As a result, the integral measure in the momentum space $\differential{V_{p}}$
becomes

$\displaystyle\differential{V_{p}}=\frac{16{{\Theta}\pqty{\tilde{\varepsilon}}}J_{\xi}}{r^{2}{{q}\pqty{\varepsilon,J_{\xi},j_{\sigma}}}\sqrt{F}}\differential{\varepsilon}\wedge\differential{J_{\xi}}\wedge\differential{j_{\sigma}}\wedge\differential{\chi}.$
(A.9)

By the definition of $(\tilde{\varepsilon},\tilde{J})$, it can be rewritten as

$\displaystyle\differential{V_{p}}=\frac{16{{\theta}\pqty{\tilde{\varepsilon}}}\tilde{J}}{r^{2}\sqrt{F}}\bqty{\frac{G}{F}\tilde{\varepsilon}^{2}-\pqty{m^{2}+\frac{4\tilde{J}^{2}}{r^{2}}}}^{-\frac{1}{2}}\differential{\tilde{\varepsilon}}\wedge\differential{\tilde{J}}\wedge\differential{j_{\sigma}}\wedge\differential{\chi}.$
(A.10)

We consider the domain of tilded variables based on the on-shell condition
(A.6). Since $j_{\sigma}$ satisfies
${j_{\sigma}}^{2}\leq{J_{\xi}}^{2}=\tilde{J}^{2}+\frac{a^{2}h}{G}{j_{\sigma}}^{2}$,
the region of $j_{\sigma}$ becomes:

$\displaystyle-\frac{\sqrt{G}}{r}\tilde{J}\leq
j_{\sigma}\leq\frac{\sqrt{G}}{r}\tilde{J}.$ (A.11)

With a fixed value of $\tilde{\varepsilon}$, the angular momentum $\tilde{J}$
is bounded by $J_{\mathrm{max}}$:

$\displaystyle\tilde{J}^{2}=\frac{r^{2}}{4}\bqty{\frac{G}{F}\tilde{\varepsilon}^{2}-(m^{2}+q^{2})}\leq\frac{r^{2}}{4}\bqty{\frac{G}{F}\tilde{\varepsilon}^{2}-m^{2}}\eqqcolon{J_{\mathrm{max}}}^{2},$
(A.12)

then we obtain $0\leq\tilde{J}\leq J_{\mathrm{max}}$.

The variable $\tilde{\varepsilon}$ takes the minimum value $m\sqrt{F/G}$ when
$j_{\sigma}$ and $\tilde{J}$ vanish. Since $\tilde{\varepsilon}$ does not have
an upper bound, the variable $\tilde{\varepsilon}$ takes the range

$\displaystyle m\sqrt{\frac{F}{G}}\leq\tilde{\varepsilon}<\infty.$ (A.13)

For the massive particles, we can normalize the variable $\tilde{\varepsilon}$
by the rest mass. Then we can define the variables as

$\displaystyle\bar{\varepsilon}\coloneqq\sqrt{\frac{G}{F}}\frac{\tilde{\varepsilon}}{m},\quad\bar{J}\coloneqq\frac{\tilde{J}}{J_{\mathrm{max}}},\quad\bar{\jmath}\coloneqq\frac{r}{\sqrt{G}}\frac{j_{\sigma}}{J_{\mathrm{max}}}.$
(A.14)

If we treat a massless particle system, we should replace
$\tilde{\varepsilon}/m$ with $\tilde{\varepsilon}$ so $\bar{\varepsilon}$ has
the dimension of energy. The integral regions of these variables are

$\displaystyle\bar{\varepsilon}\in\left[{{1},{\infty}}\right),\quad\bar{J}\in\bqty{{0},{1}},\mbox{\quad
and\quad}\bar{\jmath}\in\bqty{{-\bar{J}},{\bar{J}}},$ (A.15)

and the integral measure becomes:

$\displaystyle\differential{V_{p}}=\frac{2{{\Theta}\pqty{\bar{\varepsilon}}}k^{2}\bar{J}}{\sqrt{1-{\bar{J}}^{2}}}\differential{\bar{\varepsilon}}\wedge\differential{\bar{J}}\wedge\differential{\bar{\jmath}}\wedge\differential{\chi},$
(A.16)

where $k\coloneqq\sqrt{\bar{\varepsilon}^{2}-1}$ is the local kinetic energy
of the particle.

###  A.2 Performing the integration

We set the particle mass to unity in this subsection. In our system, the
relevant components of $p_{\mu}p_{\nu}$ are as follows:

$\displaystyle\pqty{p_{t}}^{2}$
$\displaystyle=\frac{F}{G}\pqty{\bar{\varepsilon}^{2}+\frac{2ah}{\sqrt{F}}\bar{\varepsilon}k\bar{\jmath}+\frac{a^{2}h^{2}}{F}k^{2}\bar{\jmath}^{2}},$
(A.17a) $\displaystyle p_{t}p_{\phi}$ $\displaystyle=-\varepsilon
j_{\sigma}=-\frac{\sqrt{F}}{2}\pqty{\bar{\varepsilon}k\bar{\jmath}+\frac{ah}{\sqrt{F}}k^{2}\bar{\jmath}^{2}},$
(A.17b) $\displaystyle\pqty{p_{\phi}}^{2}$
$\displaystyle={j_{\sigma}}^{2}=\frac{G}{4}k^{2}\bar{\jmath}^{2},$ (A.17c)
$\displaystyle\pqty{p_{r}}^{2}$
$\displaystyle={\mathrm{e}^{2\nu}}k^{2}(1-\bar{J}^{2}),$ (A.17d)
$\displaystyle\gamma^{\mu\nu}p_{\mu}p_{\nu}$
$\displaystyle=4J_{\xi}=r^{2}\pqty{k^{2}\bar{J}^{2}+\frac{a^{2}h}{r^{2}}k^{2}\bar{\jmath}^{2}},$
(A.17e) $\displaystyle g^{\mu\nu}p_{\mu}p_{\nu}$ $\displaystyle=-1$ (A.17f)

in the coordinate system defined by (A.14). Similarly, the distribution
function can be rewritten as

$\displaystyle{{f}\pqty{\varepsilon,j_{\sigma}}}=\exp\bqty{\alpha-\bar{\beta}(\bar{\varepsilon}-\bar{\Omega}k\bar{\jmath})}.$
(A.18)

The integration of (A.17) is expressed in the form:

$\displaystyle
2\int_{1}^{\infty}\differential{\bar{\varepsilon}}\int_{0}^{1}\differential{\bar{J}}\int_{-\bar{J}}^{\bar{J}}\differential{\bar{\jmath}}\
\bar{\varepsilon}^{i}k^{j}{\mathrm{e}^{-\bar{\beta}\bar{\varepsilon}}}\cdot\bar{J}^{m+1}(1-\bar{J}^{2})^{\pm\frac{1}{2}}\cdot\bar{\jmath}^{n}{\mathrm{e}^{\varsigma\bar{\jmath}}},$
(A.19)

where
$\varsigma\coloneqq\bar{\beta}\bar{\Omega}k=\bar{\beta}\bar{\Omega}\sqrt{\bar{\varepsilon}^{2}-1}$
with integers $m$ and $n$. To perform the integration over the momentum space,
we divide the integrations into the angular momentum sector and the energy
sector. The angular momentum sector in Eq. (A.19) is expressed as

$\displaystyle\mathcal{J}^{\pm}_{m,n}(\varsigma)=2\int_{0}^{1}\differential{\bar{J}}\bar{J}^{m+1}(1-\bar{J}^{2})^{\pm\frac{1}{2}}\int_{-\bar{J}}^{\bar{J}}\differential{\bar{\jmath}}\
\bar{\jmath}^{n}{\mathrm{e}^{\varsigma\bar{\jmath}}}.$ (A.20)

For the even number $m$, since the integrand of Eq. (A.20) is an even
function, it can be rewritten as

$\displaystyle\mathcal{J}^{\pm}_{m,n}(\varsigma)$
$\displaystyle=\int_{-1}^{1}\differential{\bar{J}}\bar{J}^{m+1}(1-\bar{J}^{2})^{\pm\frac{1}{2}}\mathcal{K}_{n}(\varsigma,\bar{J}),$
(A.21)

where $\mathcal{K}_{n}(\varsigma,\bar{J})$ is defined by:

$\displaystyle\mathcal{K}_{n}(\varsigma,\bar{J})$
$\displaystyle=\int_{-\bar{J}}^{\bar{J}}\differential{\bar{\jmath}}\
\bar{\jmath}^{n}{\mathrm{e}^{\varsigma\bar{\jmath}}}.$ (A.22)

Explicit forms of Eq. (A.22) for some specific values of $n$ are

$\displaystyle\mathcal{K}_{0}(\varsigma,\bar{J})$
$\displaystyle=\frac{{\mathrm{e}^{\varsigma\bar{J}}}-{\mathrm{e}^{-\varsigma\bar{J}}}}{\varsigma},$
(A.23a) $\displaystyle\mathcal{K}_{1}(\varsigma,\bar{J})$
$\displaystyle=-\frac{{\mathrm{e}^{\varsigma\bar{J}}}-{\mathrm{e}^{-\varsigma\bar{J}}}}{\varsigma^{2}}+\frac{\bar{J}({\mathrm{e}^{\varsigma\bar{J}}}+{\mathrm{e}^{-\varsigma\bar{J}}})}{\varsigma},$
(A.23b) $\displaystyle\mathcal{K}_{2}(\varsigma,\bar{J})$
$\displaystyle=\frac{(2+\varsigma^{2}\bar{J}^{2})({\mathrm{e}^{\varsigma\bar{J}}}-{\mathrm{e}^{-\varsigma\bar{J}}})}{\varsigma^{3}}-\frac{2\bar{J}({\mathrm{e}^{\varsigma\bar{J}}}+{\mathrm{e}^{-\varsigma\bar{J}}})}{\varsigma^{2}},$
(A.23c)

which appear in the expression of the energy-momentum tensor. In addition,
after integration over $\bar{J}$, we obtain

$\displaystyle{{\mathcal{J}^{-}_{0,0}}\pqty{\varsigma}}$
$\displaystyle=\frac{2\pi{{I_{1}}\pqty{\varsigma}}}{\varsigma},$ (A.24a)
$\displaystyle{{\mathcal{J}^{-}_{0,1}}\pqty{\varsigma}}$
$\displaystyle=\frac{2\pi{{I_{2}}\pqty{\varsigma}}}{\varsigma},$ (A.24b)
$\displaystyle{{\mathcal{J}^{-}_{0,2}}\pqty{\varsigma}}$
$\displaystyle=\frac{2\pi\pqty{{{I_{2}}\pqty{\varsigma}}+\varsigma{{I_{3}}\pqty{\varsigma}}}}{\varsigma^{2}},$
(A.24c) $\displaystyle{{\mathcal{J}^{-}_{2,0}}\pqty{\varsigma}}$
$\displaystyle=\frac{2\pi\pqty{3{{I_{2}}\pqty{\varsigma}}+\varsigma{{I_{3}}\pqty{\varsigma}}}}{\varsigma^{2}},$
(A.24d) $\displaystyle{{\mathcal{J}^{+}_{0,0}}\pqty{\varsigma}}$
$\displaystyle=\frac{2\pi{{I_{2}}\pqty{\varsigma}}}{\varsigma^{2}},$ (A.24e)

where $I_{\nu}(\varsigma)$ is the modified Bessel function of the first kind.

The integration for $\varepsilon$ cannot be analytically performed and we
perform it numerically when we solve the Einstein equations. For notational
simplicity, let us define the function
$\mathcal{I}_{i,j,\nu}(\bar{\beta},\bar{\Omega})$ as

$\displaystyle{{\mathcal{I}_{i,j,\nu}}\pqty{\bar{\beta},\bar{\Omega}}}=\int_{1}^{\infty}\differential{\bar{\varepsilon}}\
\bar{\varepsilon}^{i}\pqty{\sqrt{\bar{\varepsilon}^{2}-1}}^{j}I_{\nu}(\bar{\beta}\bar{\Omega}\sqrt{\bar{\varepsilon}^{2}-1})\exp(-\bar{\beta}\bar{\varepsilon}),$
(A.25)

where $i$, $j$ and $\nu$ are non-negative integers, $I_{\nu}(z)$ is the
modified Bessel function of the first kind, and

$\displaystyle\bar{\beta}(r)\coloneqq\beta\sqrt{\frac{F}{G}},\quad\bar{\Omega}(r)\coloneqq\frac{G}{2\sqrt{F}}\pqty{\Omega-\frac{2ah}{G}},\quad$
(A.26)

Then, the resulting forms are given as follows:

$\displaystyle T_{tt}$
$\displaystyle=\frac{4\pi^{2}{\mathrm{e}^{\alpha}}F}{G}\bqty{\frac{\mathcal{I}_{2,1,1}}{\bar{\beta}\bar{\Omega}}+\frac{2ah}{\sqrt{F}}\frac{\mathcal{I}_{1,2,2}}{\bar{\beta}\bar{\Omega}}+\frac{a^{2}h^{2}}{F}\pqty{\frac{\mathcal{I}_{0,2,2}}{{\bar{\beta}}^{2}{\bar{\Omega}}^{2}}+\frac{\mathcal{I}_{0,3,3}}{{\bar{\beta}\bar{\Omega}}}}},$
(A.27a) $\displaystyle T_{t\phi}$
$\displaystyle=-2\pi^{2}{\mathrm{e}^{\alpha}}\sqrt{F}\bqty{\frac{\mathcal{I}_{1,2,2}}{{\bar{\beta}\bar{\Omega}}}+\frac{ah}{\sqrt{F}}\pqty{\frac{\mathcal{I}_{0,2,2}}{{\bar{\beta}}^{2}{\bar{\Omega}}^{2}}+\frac{\mathcal{I}_{0,3,3}}{{\bar{\beta}\bar{\Omega}}}}},$
(A.27b) $\displaystyle T_{\phi\phi}$
$\displaystyle=\pi^{2}{\mathrm{e}^{\alpha}}G\pqty{\frac{\mathcal{I}_{0,2,2}}{{\bar{\beta}}^{2}{\bar{\Omega}}^{2}}+\frac{\mathcal{I}_{0,3,3}}{{\bar{\beta}\bar{\Omega}}}},$
(A.27c) $\displaystyle T_{rr}$
$\displaystyle=4\pi^{2}{\mathrm{e}^{\alpha}}{\mathrm{e}^{2\nu}}\frac{\mathcal{I}_{0,2,2}}{{\bar{\beta}}^{2}{\bar{\Omega}}^{2}},$
(A.27d) $\displaystyle\gamma^{\mu\nu}T_{\mu\nu}$
$\displaystyle=4\pi^{2}{\mathrm{e}^{\alpha}}r^{2}\bqty{\pqty{3+\frac{a^{2}h}{r^{2}}}\frac{\mathcal{I}_{0,2,2}}{{\bar{\beta}}^{2}{\bar{\Omega}}^{2}}+\pqty{1+\frac{a^{2}h}{r^{2}}}\frac{\mathcal{I}_{0,3,3}}{{\bar{\beta}\bar{\Omega}}}},$
(A.27e) $\displaystyle g^{\mu\nu}T_{\mu\nu}$
$\displaystyle=-4\pi^{2}{\mathrm{e}^{\alpha}}\frac{\mathcal{I}_{0,1,1}}{{\bar{\beta}\bar{\Omega}}}.$
(A.27f)

Performing the integral with respect to $\bar{\varepsilon}$ at each $r$, we
can obtain the local expression of the energy-momentum tensor.

##  Appendix B Asymptotic solutions around the center

In this Appendix, we investigate the asymptotic behavior around the center in
detail and obtain the asymptotic solutions with Vlasov matter. That is, we
determine the coefficients of metric functions in Eq. (4.8) recursively and
confirm that ${z_{2}}^{(2)}$ is a free parameter in our system.

Substituting Eqs. (4.8) to Eqs. (A.26), we obtain

$\displaystyle\bar{\beta}(x)\simeq\beta_{c}\bqty{1+\frac{{u}^{(2)}}{4}x^{2}+\frac{{u}^{(3)}}{12}x^{3}+\order{x^{4}}},$
(B.1a)
$\displaystyle\bar{\Omega}(x)\simeq\frac{1}{2}\Omega_{c}x\bqty{1-\frac{u^{(2)}+w^{(2)}}{4}x^{2}}+\order{x^{4}},$
(B.1b)

where $\beta_{c}\coloneqq\beta\sqrt{{y_{3}}^{(0)}}$,
$\Omega_{c}\coloneqq\Omega/\sqrt{{y_{3}}^{(0)}}$,

$\displaystyle{u}^{(i)}\coloneqq\frac{{y_{3}}^{(i)}}{{y_{3}}^{(0)}}\mbox{\quad
and\quad}w^{(i)}\coloneqq\frac{4{z_{1}}^{(i)}}{\Omega_{c}\sqrt{{y_{3}}^{(0)}}}-{z_{2}}^{(i)}$
(B.2)

are constants. Then we have

$\displaystyle
I_{\nu}(uk)\exp(-\bar{\beta}\bar{\varepsilon})\simeq\frac{{\mathrm{e}^{-\beta_{c}\bar{\varepsilon}}}}{\Gamma(1+\nu)}\pqty{\frac{\beta_{c}\Omega_{c}}{4}kx}^{\nu}\bqty{1+\frac{\zeta_{2}}{4}x^{2}+\frac{\zeta_{3}}{12}x^{3}+\order{x^{4}}}$
(B.3)

with

$\displaystyle\zeta_{2}\coloneqq\frac{{\beta_{c}}^{2}{\Omega_{c}}^{2}}{4(1+\nu)}k^{2}-\beta_{c}u^{(2)}\bar{\varepsilon}-\nu
w^{(2)},\quad\zeta_{3}\coloneqq-\beta_{c}u^{(3)}\bar{\varepsilon}-\nu
w^{(3)}.$ (B.4)

Therefore we can obtain the asymptotic expression as

$\displaystyle\mathcal{I}_{i,j,\nu}\simeq\frac{1}{\Gamma(1+\nu)}\bqty{\frac{\beta_{c}\Omega_{c}}{4}x}^{\nu}\int_{1}^{\infty}\differential{\bar{\varepsilon}}\
\bar{\varepsilon}^{i}k^{j+\nu}{\mathrm{e}^{-\beta_{c}\bar{\varepsilon}}}\bqty{1+\frac{\zeta_{2}}{4}x^{2}+\frac{\zeta_{3}}{12}x^{3}+\order{x^{4}}}.$
(B.5)

Using this asymptotic expression, the central value of $T_{tt}$ becomes

$\displaystyle
T_{tt}^{c}=\frac{4\pi^{2}{\mathrm{e}^{\alpha}}(2+\beta_{c})(6+3\beta_{c}+{\beta_{c}}^{2})}{{\beta_{c}}^{5}}{y_{3}}^{(0)}{\mathrm{e}^{-\beta_{c}}}.$
(B.6)

Thus we can expand the matter sector in the Einstein equations as

$\displaystyle
s_{\mu\nu}=\sum_{n=0}^{\infty}\frac{s_{\mu\nu}^{(n)}}{n!}x^{n},\quad\gamma^{\mu\nu}s_{\mu\nu}=\sum_{n=0}^{\infty}\frac{s_{\gamma}^{(n)}}{n!}x^{n},\quad
g^{\mu\nu}s_{\mu\nu}=\sum_{n=0}^{\infty}\frac{s^{(n)}}{n!}x^{n}.$ (B.7)

We can calculate the relevant non-vanishing components of them as

$\displaystyle s_{tt}^{(0)}$
$\displaystyle=\bqty{\frac{2\pqty{18+\beta_{c}\pqty{18+\beta_{c}\pqty{7+\beta_{c}}}}}{3\pqty{2+\beta_{c}}\pqty{6+\beta_{c}\pqty{3+\beta_{c}}}{y_{3}}^{(0)}}+\frac{4}{\lambda^{2}}}{y_{3}}^{(0)},$
(B.8a) $\displaystyle s_{rr}^{(0)}$
$\displaystyle=\frac{9+\beta_{c}\pqty{9+\beta_{c}\pqty{4+\beta_{c}}}}{3\pqty{2+\beta_{c}}\pqty{6+\beta_{c}\pqty{3+\beta_{c}}}{y_{3}}^{(0)}}-\frac{4}{\lambda^{2}},$
(B.8b) $\displaystyle s^{(0)}$
$\displaystyle=\frac{2{\beta_{c}}^{2}\pqty{1+{\beta_{c}}}}{3\pqty{2+\beta_{c}}\pqty{6+\beta_{c}\pqty{3+\beta_{c}}}{y_{3}}^{(0)}}-\frac{20}{\lambda^{2}},$
(B.8c) $\displaystyle s_{t\phi}^{(2)}$
$\displaystyle=-\frac{15+\beta_{c}(15+6\beta_{c}+{\beta_{c}}^{2})}{2(2+\beta_{c})(6+3\beta_{c}+{\beta_{c}}^{2})\sqrt{{y_{3}}^{(0)}}}\Omega_{c},$
(B.8d)

and $s_{\gamma}^{(2)}=12s_{\phi\phi}^{(2)}=3s_{rr}^{(0)}$ at the leading
order. Substituting these expressions into the Einstein equations and solving
them order by order, we obtain

$\displaystyle{z_{1}}^{(2)}$
$\displaystyle=-\frac{15+\beta_{c}(15+6\beta_{c}+{\beta_{c}}^{2})}{6(2+\beta_{c})(6+3\beta_{c}+{\beta_{c}}^{2})\sqrt{{y_{3}}^{(0)}}}\Omega_{c},$
(B.9) $\displaystyle{y_{3}}^{(2)}$
$\displaystyle=\frac{18+\beta_{c}\pqty{18+\beta_{c}\pqty{7+\beta_{c}}}}{3\pqty{2+\beta_{c}}\pqty{6+\beta_{c}\pqty{3+\beta_{c}}}}+\frac{2}{\lambda^{2}}{y_{3}}^{(0)},$
(B.10) $\displaystyle{y_{4}}^{(2)}$
$\displaystyle={z_{2}}^{(2)}+\frac{1}{3{y_{3}}^{(0)}}-\frac{2}{\lambda^{2}},$
(B.11)

at the leading order. The value of ${z_{2}}^{(2)}$ is not fixed in this
expression and is regarded as a free parameter. All other coefficients can be
determined by the equations recursively. Therefore, the boundary conditions of
the system are Eq. (4.10) and the parameters are
$({y_{3}}^{(0)},{z_{2}}^{(2)},\beta,\Omega,\lambda)$ as a whole.

## References

  * [1] G. Rein and A. D. Rendall.  Global existence of solutions of the spherically symmetric Vlasov-Einstein system with small initial data.  Communications in Mathematical Physics, 150(3):561 – 583, 1992\.
  * [2] G. Rein.  Static solutions of the spherically symmetric Vlasov-Einstein system.  Math. Proc. Cambridge Phil. Soc., 115:559, 1994.
  * [3] Håkan Andréasson, David Fajman, and Maximilian Thaller.  Static solutions to the Einstein-Vlasov system with non-vanishing cosmological constant.  SIAM J. Math. Anal., 47(4):2657–2688, 2014.
  * [4] Håkan Andréasson, David Fajman, and Maximilian Thaller.  Models for Self-Gravitating Photon Shells and Geons.  Annales Henri Poincare, 18(2):681–705, 2017.
  * [5] Hakan Andreasson and Gerhard Rein.  A Numerical investigation of the stability of steady states and critical phenomena for the spherically symmetric Einstein-Vlasov system.  Class. Quant. Grav., 23:3659–3678, 2006.
  * [6] MAHIR HADŽ IĆ and GERHARD REIN.  Stability for the spherically symmetric einstein–vlasov system–a coercivity estimate.  Mathematical Proceedings of the Cambridge Philosophical Society, 155(3):529–556, sep 2013.
  * [7] Mahir Hadzic, Zhiwu Lin, and Gerhard Rein.  Stability and instability of self-gravitating relativistic matter distributions.  10 2018.
  * [8] Sebastian Günther, Jacob Körner, Timo Lebeda, Bastian Pötzl, Gerhard Rein, Christopher Straub, and Jörg Weber.  A numerical stability analysis for the Einstein-Vlasov system.  Class. Quant. Grav., 38(3):035003, 2021.
  * [9] Sebastian Günther, Christopher Straub, and Gerhard Rein.  Collisionless Equilibria in General Relativity: Stable Configurations beyond the First Binding Energy Maximum.  Astrophys. J., 918(2):48, 2021.
  * [10] Mahir Hadžić and Zhiwu Lin.  Turning point principle for relativistic stars.  Communications in Mathematical Physics, 387(2):729–759, sep 2021\.
  * [11] Hakan Andreasson.  The Einstein-Vlasov System/Kinetic Theory.  Living Rev. Rel., 14:4, 2011.
  * [12] Håkan Andréasson.  Existence of Steady States of the Massless Einstein–Vlasov System Surrounding a Schwarzschild Black Hole.  Annales Henri Poincare, 22(12):4271–4297, 2021.
  * [13] V. A. ANTONOV.  Most probable phase distribution in spherical star systems and conditions for its existence.  Vest. Leningrad Univ., 7:135, 1962.
  * [14] D. Lynden-Bell and Roger Wood.  The gravo-thermal catastrophe in isothermal spheres and the onset of red-giant structure for stellar systems.  mnras, 138:495, Jan 1968.
  * [15] Hiroki Asami and Chul-Moon Yoo.  Thermal equilibrium states and instability of self-gravitating particles in an asymptotically AdS spacetime.  Classical and Quantum Gravity, 38(16):165014, jul 2021.
  * [16] Hiroki Asami and Chul-Moon Yoo.  Gravothermal catastrophe and critical dimension in a $D$–dimensional asymptotically AdS spacetime.  Phys. Rev. D, 106(4):044065, 2022.
  * [17] Stuart L. Shapiro and Saul A. Teukolsky.  Relativistic Stellar Systems with Spindle Singularities.  Astrophys. J., 419:622, December 1993.
  * [18] Stuart L. Shapiro and Saul A. Teukolsky.  Relativistic Stellar Systems with Rotation.  Astrophys. J., 419:636, December 1993.
  * [19] Hakan Andreasson, Markus Kunze, and Gerhard Rein.  Rotating, stationary, axially symmetric spacetimes with collisionless matter.  Commun. Math. Phys., 329:787–808, 2014.
  * [20] Ellery Ames, Håkan Andréasson, and Anders Logg.  On Axisymmetric and Stationary Solutions of the Self-Gravitating Vlasov System.  Class. Quant. Grav., 33(15):155008, 2016.
  * [21] Maximilian Thaller.  Rotating Clouds of Charged Vlasov Matter in General Relativity.  Class. Quant. Grav., 37(3):035008, 2020.
  * [22] R.C Myers and M.J Perry.  Black holes in higher dimensional space-times.  Annals of Physics, 172(2):304–347, 1986.
  * [23] S. W. Hawking, C. J. Hunter, and Marika Taylor.  Rotation and the AdS / CFT correspondence.  Phys. Rev. D, 59:064005, 1999.
  * [24] G. W. Gibbons, H. Lu, Don N. Page, and C. N. Pope.  The General Kerr-de Sitter metrics in all dimensions.  J. Geom. Phys., 53:49–73, 2005.
  * [25] G. W. Gibbons, H. Lu, Don N. Page, and C. N. Pope.  Rotating black holes in higher dimensions with a cosmological constant.  Phys. Rev. Lett., 93:171102, 2004.
  * [26] Hari K. Kunduri, James Lucietti, and Harvey S. Reall.  Gravitational perturbations of higher dimensional rotating black holes: Tensor perturbations.  Phys. Rev. D, 74:084021, 2006.
  * [27] Keiju Murata, Tatsuma Nishioka, and Norihiro Tanahashi.  Warped AdS(5) Black Holes and Dual CFTs.  Prog. Theor. Phys., 121:941–957, 2009.
  * [28] Juan Maldacena.  International Journal of Theoretical Physics, 38(4), 1999.
  * [29] S.S. Gubser, I.R. Klebanov, and A.M. Polyakov.  Gauge theory correlators from non-critical string theory.  Physics Letters B, 428(1-2), May 1998.
  * [30] Edward Witten.  Anti de sitter space and holography, 1998.
  * [31] Piotr BizoÅ and Andrzej Rostworowski.  Weakly turbulent instability of antiâde sitter spacetime.  Physical Review Letters, 107(3), Jul 2011.
  * [32] Oscar J. C. Dias, Gary T. Horowitz, Don Marolf, and Jorge E. Santos.  On the Nonlinear Stability of Asymptotically Anti-de Sitter Solutions.  Class. Quant. Grav., 29:235019, 2012.
  * [33] Maciej Maliborski and Andrzej Rostworowski.  Time-Periodic Solutions in an Einstein AdS–Massless-Scalar-Field System.  Phys. Rev. Lett., 111:051102, 2013.
  * [34] Alex Buchel, Steven L. Liebling, and Luis Lehner.  Boson stars in AdS spacetime.  Phys. Rev. D, 87(12):123006, 2013.
  * [35] Venkat Balasubramanian, Alex Buchel, Stephen R. Green, Luis Lehner, and Steven L. Liebling.  Holographic Thermalization, Stability of Anti–de Sitter Space, and the Fermi-Pasta-Ulam Paradox.  Phys. Rev. Lett., 113(7):071601, 2014.
  * [36] Venkat Balasubramanian, Alex Buchel, Stephen R. Green, Luis Lehner, and Steven L. Liebling.  Reply to Comment on “Holographic Thermalization, Stability of Anti–de Sitter Space, and the Fermi-Pasta-Ulam Paradox”.  Phys. Rev. Lett., 115(4):049102, 2015.
  * [37] Piotr Bizoń, Maciej Maliborski, and Andrzej Rostworowski.  Resonant Dynamics and the Instability of Anti–de Sitter Spacetime.  Phys. Rev. Lett., 115(8):081103, 2015.
  * [38] Fotios Dimitrakopoulos and I-Sheng Yang.  Conditionally extended validity of perturbation theory: Persistence of AdS stability islands.  Phys. Rev. D, 92(8):083013, 2015.
  * [39] Stephen R. Green, Antoine Maillard, Luis Lehner, and Steven L. Liebling.  Islands of stability and recurrence times in AdS.  Phys. Rev. D, 92(8):084001, 2015.
  * [40] David Garfinkle and Leopoldo A. Pando Zayas.  Rapid Thermalization in Field Theory from Gravitational Collapse.  Phys. Rev. D, 84:066006, 2011.
  * [41] Joanna Jalmuzna, Andrzej Rostworowski, and Piotr Bizon.  A Comment on AdS collapse of a scalar field in higher dimensions.  Phys. Rev. D, 84:085021, 2011.
  * [42] P. Bizoń and A. Rostworowski.  Gravitational Turbulent Instability of AdS5.  Acta Phys. Polon. B, 48:1375, 2017.
  * [43] Ben Craps, Oleg Evnin, and Joris Vanhoof.  Renormalization group, secular term resummation and AdS (in)stability.  JHEP, 10:048, 2014.
  * [44] Ben Craps, Oleg Evnin, and Joris Vanhoof.  Renormalization, averaging, conservation laws and AdS (in)stability.  JHEP, 01:108, 2015.
  * [45] Oleg Evnin.  Resonant Hamiltonian systems and weakly nonlinear dynamics in AdS spacetimes.  Class. Quant. Grav., 38(20):203001, 2021.
  * [46] Georgios Moschidis.  A proof of the instability of AdS for the Einstein–null dust system with an inner mirror.  Anal. Part. Diff. Eq., 13(6):1671–1754, 2020.
  * [47] Georgios Moschidis.  The Einstein–null dust system in spherical symmetry with an inner mirror: structure of the maximal development and Cauchy stability.  4 2017.
  * [48] Vladislav Vaganov.  Self-gravitating radiation in AdS(d).  7 2007.
  * [49] John Hammersley.  A Critical dimension for the stability of radiating perfect fluid stars.  Class. Quant. Grav., 25:205010, 2008.


Institute of Natural Sciences, Westlake Institute for Advanced Study,
Hangzhou, Zhejiang 310024, China RENEW (Research and Education in Energy,
Environment and Water) Institute, University at Buffalo, The State University
of New York, Buffalo, NY 14260, USA Institute of Natural Sciences, Westlake
Institute for Advanced Study, Hangzhou, Zhejiang 310024, China

# Nonadiabatic quantum dynamics of tribovoltaic effects at sliding metal-
semiconductor interfaces

Guangming Liu  School of Science, Westlake University, Hangzhou, Zhejiang
310024, China  Jun Liu<EMAIL_ADDRESS>Department of Mechanical and
Aerospace Engineering, University at Buffalo, The State University of New
York, Buffalo, NY 14260, USA  Wenjie Dou  School of Science, Westlake
University, Hangzhou, Zhejiang 310024, China<EMAIL_ADDRESS>
###### Abstract

Recent experiments observe electric current generation at a sliding metal-
semiconductor interfaces. Here, we present a detailed theoretical study on how
electric voltage is generated at such a sliding interface. Our study is based
on a two-band Anderson-Holstein model, and we solve the coupled electron-
phonon dynamics using a surface hopping method. We show that the high local
temperature induced by mechanic motion at the interfaces could lead to
electron-hole pair generation through electron-phonon couplings. We quantify
the efficiency of electron-hole generation as well as electric voltage as a
function of local temperatures and semiconductor bandgaps. We find that
increasing the local temperatures can lead to higher electron-hole generations
and electric voltage. Furthermore, we find that there is a turnover for the
electric voltage as a function of the bandgap. Such an observation is in
agreement with the experimental results. Our study offers a theoretical
framework to understand tribovoltaic effects from a quantum mechanical point
of view, and our approach can be used to quantitively simulate realistic
sliding metal-semiconductor junctions.

##  1 Introduction

When a semiconductor material is involved in a sliding contact, a high direct-
current (DC) can be generated, which has been recently observed in various
contact material system (e.g. metal/semiconductor,1, 2, 3
metal/insulator/semiconductor, 3, 4 p-n junction,5, 6liquid/semicond- uctor7,
8, 9). Such a phenomenon is referred to as ’tribovoltaic’ effect10. Moreover,
new multi-physics phenomenon originated from the tribovoltaic effect such as
tribo-photovoltaic effect11, 4, 9 and tribovoltaic-thermoelectric effect 12
have been found successively in dynamic metal-semiconductor Schottky systems.
Although the fundamental mechanism remains unsolved, it is found the DC
current density output of such systems is on the order of 10-100 $A/m^{2}$,
which is 3-4 orders’ higher than that in traditional piezoelectric or
triboelectric nanogenerator and thus show a great promise on next-generation
mechanical energy harvesting devices for self-powered electronics and Internet
of Things (IoTs) sensors13.

In the past decades, tremendous efforts have been devoted to understanding and
quantifying charge transfer during contact. To date, the fundamental physics
of triboelectricity remains unsolved14, 15, 16. The challenge of an explicit
understanding of the mechanism resides in the complex multi-scale and multi-
physics interaction at a contact interface. 17, 18, 19, 20 Frictional energy
dissipation in materials may manifest in different ways: heat generation, exo-
electron, luminescence, electron-hole generation, surface charges, etc21. In
general situation, the majority of frictional energy ends up with heat either
through phonon-phonon interaction or relaxation of other electronic
excitations. When a Schottky contact consisting of a metal and a semiconductor
material is presented however, it is proposed that the input frictional energy
may induce a non-adiabatic electronic excitation of energetic carriers, and
the generated electron-hole pairs may be non-adiabatic electronic excitation,
which ultimately contribute to the DC generation at a sliding Schottky
contact22, 10.

The physical essence of charge transfer at the interfaces is the breakdown of
Born-Oppenheimer (BO) approximation: In the BO approximation, nuclear motion
is assumed to be much slower than electronic motion, such that nuclei move on
the adiabatic potential energy surface. Under external stimuli, the
accelerated nuclear motion can introduce non-adiabatic electronic transitions
between different potential energy surface near the crossing point such that
the BO approximation breaks23, 24. Nevertheless, a basic understanding of the
electronic excitation process in an electro-mechanical coupling system is
still missing. In this work, a systematic study on the non-adiabatic
electronic excitation at a metal/semiconductor sliding interface using quantum
dynamics approach is presented, in an effort to shed light on fundamental
process and provide instructional guide for materials design and optimization
for future semiconductor-based nanogenerators.

In our previous studies, we have presented a surface hopping dynamics to
describe nonadiabatic dynamics near a metal surface. 25, 26 In such a
dynamical method, we propagate trajectories on different potential energy
surfaces describing different charge states. We further introduce stochastic
hopping between them indicating charge transfer between these states. Note
that, different Tully’s surface hopping method, which is typically used to
model nonadiabatic dynamics in solution or in gas phase, our surface hopping
is intended to study nonadiabatic dynamics at metal surfaces. Our surface
hopping dynamics have been successfully used to describe electron transfer,
energy relaxation, current-voltage characteristic near metal surfaces as well
as in nano-junctions. 27, 28

In this manuscript, we present a two-level Anderson-Holstein model based on
electron-phonon couplings to describe the semiconductor-metal interfaces. We
use a surface hopping method to describe the nonadiabatic dynamics at such
interfaces. Our dynamics indicate that local temperature induced by mechanic
motion at the interfaces could lead to electron-hole pair generation through
electron-phonon couplings. We can then quantify the generation of the electron
and hole population as well as the average voltage. We find that there is a
turnover of the voltage as the function of the bandgap. This finding is in
agreement with experimental results. We believe that this article provides
atomic insights on how electric current is generated on sliding semiconductor-
metal interfaces. Our approach can be used to study realistic sliding metal-
semiconductor junctions in the future.

##  2 Theoretical model

Here, we use a two-level/two-band Anderson-Holstein model to describe
semiconductor/metal interfaces. Our total Hamiltonian consists of three parts:
$\hat{H}_{M}$ for the metal, $\hat{H}_{S}$ for the semiconductor, and
$\hat{H}_{I}$ for the couplings between the metal and semiconductors, i.e.

$\displaystyle\hat{H}=\hat{H}_{M}+\hat{H}_{S}+\hat{H}_{I}$ (1)

The metal consists of a manifold of electronic states,

$\displaystyle\hat{H}_{M}=\sum_{k}\epsilon_{k}\hat{c}^{+}_{k}\hat{c}_{k}$ (2)

For simplicity, we model the semiconductor with two levels, conduction
level/band $E_{c}(x)$, and valence level/band $E_{v}(x)$. Here we have a local
phonon degrees of freedom ($x$ and $p$ are the position and momentum of the
phonon) coupled to the electronic levels. Such that the Hamiltonian for the
semiconductor can be written as

$\displaystyle\hat{H}_{S}=E_{c}(x)\hat{d}_{c}^{+}\hat{d}_{c}+E_{v}(x)\hat{d}_{v}^{+}\hat{d}_{v}+U_{0}(x)+\frac{p^{2}}{2m}$
(3)

$U_{0}(x)$ is the potential for the local phonon mode. Below, we will
approximate the potential as a harmonic oscillator,
$U_{0}=\frac{1}{2}m\omega^{2}x^{2}$ ($\omega$ is the frequency of the
oscillator). The couplings between the metal and levels in the semi-conductor
are assumed to be bilinear:

$\displaystyle\hat{H}_{I}=\sum_{k}V_{ck}(\hat{c}^{+}_{k}d_{c}+d_{c}^{+}\hat{c}_{k})+\sum_{k}V_{vk}(\hat{c}^{+}_{k}d_{v}+d_{v}^{+}\hat{c}_{k})$
(4)

We can define the hybridization functions to characterize the strength of the
couplings:

$\displaystyle\Gamma_{c}=\sum_{k}|V_{ck}|^{2}\delta(\epsilon-\epsilon_{k})$
(5)
$\displaystyle\Gamma_{v}=\sum_{k}|V_{vk}|^{2}\delta(\epsilon-\epsilon_{k})$
(6)

We further assume that the conductor band and the valence band couples to the
local phonon as the following:

$\displaystyle E_{c}(x)=gx\sqrt{m\omega/\hbar}+E_{g}/2$ (7) $\displaystyle
E_{v}(x)=gx\sqrt{m\omega/\hbar}-E_{g}/2$ (8)

Here, $g$ characterizes the electron-phonon coupling strength. We have assumed
that the coupling is linearly proportional to local phonon displacement, which
is the typical approximation for the electron-phonon interactions. $E_{g}$ is
the band gap between the conduction band and valance band.

The above equations conclude our model. Our goal is to simulate the electron
transfer processes between the conduction/valence band and metal surface. Such
processes could lead to electron-hole pair generation in the semiconductor as
well as electric voltage. Obviously, solving the real time dynamics will be
very difficult. Below, we introduce a classical master equation and surface
hopping (SH) method to simulate the dynamical processes.

(a)

(b)

(c)

Figure 1: Semiconductor near a metal surface. (a)Sliding semiconductor near a
metal surface create electron and hole in the semiconductor. (b)Four states
will arise after surface hopping. $U_{0}$: no electron in valance band(VB) or
conduct band(CB), $U_{0}+E_{v}$: one electron occupied in VB. $U_{0}+E_{c}$:
one electron occupied in CB. $U_{0}+E_{c}+E_{v}$: two electrons occupied in CB
and VB. (c)There will be electron current flow when connected to circus. Our
model explain why electron-hole can be created at the surfaces.

###  2.1 Classical master equation

To solve the coupled electron-phonon motion at the interfaces for the two
level Anderson-Holstein model introduced above, we explore the classical
master equation method. In the limit of weak coupling ($\Gamma<kT$), provided
phonon is classical ($\hbar\omega<kT$), we can derive a classical master
equation (CME) to describe the dynamics. The exact derivation is similar to
the one for one-level Anderson-Holstein model.25, 26 Here we outline the
derivation and show the equations of motion.

We start from the Liouville equation for the total system. Under the weak
coupling assumption, after tracing over the bath degrees of freedom, we arrive
at the quantum master equation for the system only:

$\displaystyle\partial_{t}\hat{\rho}=i[\hat{H}_{s},\hat{\rho}]-\hat{\hat{\mathcal{L}}}\hat{\rho}$
(9)

Here, $\hat{\rho}$ is density operator for the system.
$\hat{\hat{\mathcal{L}}}$ is the Lindblad superoperator coming from the
system-bath coupling. Note that, the system includes the electronic states in
the semiconductor as well as the phonon degrees of freedom. We further make
classical assumption for the phonon using Wigner transformation:

$\displaystyle\hat{\rho}(x,p)=\int d\Delta x\>\langle x-\Delta
x|\hat{\rho}|x+\Delta x\rangle\exp(-ip\Delta x)$ (10)

Here, $\hat{\rho}(x,p)$ can be interpreted as the phase space density for
different electronic states. Particularly, we have $\rho_{0}$, $\rho_{v}$,
$\rho_{c}$, $\rho_{vc}$ to denote the phase space densities with empty
electron, one electron in the valance band, one electron in the conduction
band, and two electrons, respectively. The equations of motion for these
densities are:

$\displaystyle\partial_{t}\rho_{0}=\frac{p}{m}\partial_{x}\rho_{0}-\partial_{x}U_{0}\partial_{p}\rho_{0}-(\Gamma_{c}f(E_{c})+\Gamma_{v}f(E_{v}))\rho_{0}+\Gamma_{c}f(-E_{c})\rho_{c}+\Gamma_{v}f(-E_{v})\rho_{v}$
(11)
$\displaystyle\partial_{t}\rho_{v}=\frac{p}{m}\partial_{x}\rho_{v}-\partial_{x}U_{1}\partial_{p}\rho_{v}+\Gamma_{v}f(E_{v})\rho_{0}-(\Gamma_{c}f(E_{c})+\Gamma_{v}f(-E_{v}))\rho_{v}+\Gamma_{c}f(-E_{c})\rho_{vc}$
(12)
$\displaystyle\partial_{t}\rho_{c}=\frac{p}{m}\partial_{x}\rho_{c}-\partial_{x}U_{2}\partial_{p}\rho_{c}+\Gamma_{c}f(E_{c})\rho_{0}-(\Gamma_{c}f(-E_{c})+\Gamma_{v}f(E_{v}))\rho_{c}+\Gamma_{v}f(-E_{v})\rho_{vc}$
(13)
$\displaystyle\partial_{t}\rho_{vc}=\frac{p}{m}\partial_{x}\rho_{cv}-\partial_{x}U_{3}\partial_{p}\rho_{cv}+\Gamma_{c}f(E_{c})\rho_{c}+\Gamma_{v}f(E_{v})\rho_{v}-(\Gamma_{c}f(-E_{c})+\Gamma_{v}f(-E_{v}))\rho_{vc}$
(14)

Here, $U_{1}=U_{0}+E_{v}$, $U_{2}=U_{0}+E_{c}$, and $U_{3}=U_{0}+E_{c}+E_{v}$
correspond to the potential energy surfaces for one electron in the valance
band, one electron in the conduction band, and two electrons, respectively. We
note that, the densities follow the classical motion on the corresponding
potential energy surfaces with transitions between different electronic
states. The hopping rates are proportional to the hybridization function
($\Gamma_{c}$ or $\Gamma_{v}$) as well as the Fermi functions. In Figure 1, we
sketch a semiconductor with two levels near a continuum of electronic states
(metal surface). For simplicity, we do not include electron-phonon couplings
in our plots. Figure 1(b) show the four electronic states with different
electron occupations. We refer to the state with one electron on the valence
band as our ground state, and the others as excited states. In particular, the
state with no electron is referred to as one hole state, and the state with
one electron on the conduction band is referred to as electron-hole pair
state.

In Fig. 2, we plot the potential energy surfaces for different electronic
states. Notice that there are several crossings between the PESs when the
bandgap $E_{g}$ is small. When the bandgap is large, the crossing can only
happen at very high energies. At the crossing points, the hopping happens most
frequently. At very low nuclear energy, the electrons will likely remain on
the ground state. With higher nuclear energy, there will be more transitions
on the excited states. Below, we will use a surface hopping algorithm to
propagate the equation of motion and analyze non-adiabatic electron transfer
between semiconductor and metal surfaces.

(a)

(b)

Figure 2: potential energy surfaces for different charge states: no electron
in either conduction or valance band $U_{0}$, one electron in the valence band
$U_{1}=U_{0}+E_{v}$, one electron in the conduction band $U_{2}=U_{0}+E_{c}$,
and two electrons in the conduction and valance band
$U_{3}=U_{0}+E_{v}+E_{c}$. Parameters: (a) $\hbar\omega=0.03\ eV$,
$g_{c}=g_{v}=0.09\ eV$, $E_{g}=1\ eV$ (b) $\hbar\omega=0.03\ eV$,
$g_{c}=g_{v}=0.09\ eV$,$E_{g}=2\ eV$.

##  3 Numerical Results

We now use a surface hopping algorithm to solve the CME dynamics. In the
surface hopping method, we use a swarm of trajectories to represent nuclear
densities. We evolve trajectories on the active potential energy surfaces
(PESs), and we introduce hopping between the PESs. To do so, at each time
step, we generate a random number from 0 to 1 from uniform distribution; we
then compare the random number to hopping rate times $dt$ ($dt$ is the time
interval). A hopping to a different potential energy surface occurs when the
random number is smaller than the hopping rate times $dt$; otherwise, the
trajectory remains on the original PES.

In a sliding metal-semiconductor junction, the mechanic motion could induce
phonon excitation, especially at the local contact asperities, which have a
typical size of 10-100 nm. The single asperities at contact interfaces
experience a large local pressure at GPa level and a local contact temperature
$T$ as high as 1000K.29 Such excited phonon motion could then introduce
electron-hole pairs as well as voltage generation through electron-phonon
couplings. To mimic a such process, we initialize all trajectories on the
ground state (one electron in the valance band, with PES $U_{1}=U_{0}+E_{v}$),
where the momentums and positions of these trajectories satisfy a Boltzmann
distribution at a given initial phonon temperature $T_{i}$,

$\displaystyle\rho_{v}(x,p)$ $\displaystyle=$ $\displaystyle\frac{\omega}{2\pi
kT_{i}}\exp\Big{(}-\frac{1}{kT_{i}}\big{(}\frac{p^{2}}{2m}+\frac{1}{2}m\omega^{2}(x-x_{0})^{2}\big{)}\Big{)}$
(15) $\displaystyle\rho_{0}$ $\displaystyle=$
$\displaystyle\rho_{c}=\rho_{vc}=0$ (16)

Here, $x_{0}=-\sqrt{2}g/\hbar\omega$ is the position of lowest energy on the
ground state. We then evolve our equation of motion to see how phonon energy
relaxes and how electron-hole pairs are generated as a function of time. In
our simulation, we use 4th order Runge-Kutta to propagate the dynamics. Unless
otherwise stated, we use 10000 trajectories for our SH simulations.

###  3.1 Phonon Relaxation

We first look at the phonon relaxation at the interfaces. In Figure 3(a), we
initialize our phase space density on the ground state with different phonon
temperatures $T_{i}$, and we plot average kinetic energy as a function of
time. Here, we fix the temperature of the metal (electron bath, $kT=0.05eV$),
which only appears in the equation of motion through the Fermi distribution in
the hopping rates. we note that, the kinetic energy of the phonon will reach
to the same steady state, with the average kinetic energy being half $kT$,
i.e. the temperature of the metal, regardless of the initial conditions. In
fact, we can verify that, the steady state nuclear distribution will obey a
Boltzmann distribution with the same temperature of the metal (see Ref. 25,
26). Physically, this observation indicates that when the phonon degrees of
freedom interact with a manifold of electronic states long enough, the phonon
will reach to an equilibrium distribution with the same temperature as the
electron bath. Note also that, the relaxation rate is independent of the
initial kinetic energy. As shown in Ref. 27, 30, the relaxation rate strongly
depends on the electron-phonon coupling $g$, the timescales of electronic
motion $\Gamma_{c}$ and/or $\Gamma_{v}$, and the timescale of the nuclear
motion $\hbar\omega$.

(a)

(b)

Figure 3: Phonon relaxation:$\hbar\omega=0.03\ eV,g=0.09\ eV,kT=0.05\ eV$,
$\Gamma_{c}=\Gamma_{v}=0.02\ eV$. (a)Phonon Relaxation with different initial
conditions, $T_{i}$ represents the temperature at which the oscillator is
initialized. (b)Phonon relaxation with different Bandgaps, $E_{g}(eV)$.

In Figure 3(b), we plot the phonon relaxation for different band gaps $E_{g}$.
We have fixed the initial phonon temperature, $T_{i}=10T$. Here, we observe a
turnover of the relaxation timescales for different Bandgaps. 1) For the case
where band gaps $E_{g}>1eV$, the relaxation takes longer for larger bandgaps.
In fact, for the Bandgap as large as $E_{g}=4eV$, we do not see energy
relaxation at all. This is because when the band gap is large enough, the
crossing between PESs will only occur at very large potentials (see Fig. 2),
such that hopping events are very rare, which results in little or no energy
relaxation. 2) For the cases where band gaps $E_{g}<1eV$, the relaxation
timescales could also be larger for smaller bandgaps. This is due to the fact
that when hopping occurs, the energy loss is small if the bandgaps are small.
The relaxation timescale is a result of two competing effects: the frequency
of hopping and energy loss when hopping events occur. Increasing the bandgap
will reduce the frequency of hopping, but increase energy loss when hopping
events occur.

(a)

(b)

(c)

Figure 4: Electron and hole populations at different local temperatures.
(a)Electron population as a function of time. (b)Hole population as a function
of time. (c)Steady state population of electrons and holes as a function of
local temperatures. Parameters: $\hbar\omega=0.03\ eV,g=0.09\ eV,E_{g}=1\ eV$,
$\Gamma_{c}=\Gamma_{v}=0.02\ eV$.

###  3.2 Electron and hole populations

In a sliding metal-semiconductor interfaces, the open-circuit voltage
($V_{oc}$) generation (tribovoltaic potential difference between the metal and
the semiconductor) is associated with the electron-hole (e-h) pair generation
in a mechanically excited system.22, 10 To quantify the electron-hole (e-h)
pair generation, we define the electron population on the conduction band
$N_{c}$ as well as the hole population on the valence band $N_{h}$ as
following:

$\displaystyle N_{c}=\int\big{(}\rho_{c}(x,p)+\rho_{vc}(x,p)\big{)}dxdp,$ (17)
$\displaystyle N_{h}=\int\big{(}\rho_{0}(x,p)+\rho_{c}(x,p)\big{)}dxdp$ (18)

Here again, $\rho_{c}(x,p)$ is the phase space density for the one electron
and one hole state, $\rho_{vc}(x,p)$ is the density for two electron states,
and $\rho_{0}(x,p)$ is the density for one hole state (see Figure 2).

In Figure 4(a) and Figure 4(b), we plot hole and electron populations at
different temperatures as a function of time. Here the initial phonon
temperature $T_{i}$ is set to be equal to the temperature of the electron or
the local temperature of the asperities $T$. We further vary the contact
temperature at the single asperities depending on mechanical motion parameters
(e.g. pressure and speed) at a sliding metal-semiconductor junction. We find
that the timescales for electron and hole generation are roughly the same. In
general, the steady state electron generation is larger than the hole
generation. When the local temperature is small enough($\leq 0.1eV$), there
will be nearly no hole generation. This can be seen from Figure 2(a): $U_{3}$
is the lowest excited states, whereas $U_{0}$ and $U_{2}$ are excited states
with higher energies. At small local temperatures, there will be little
distribution on potential surface $U_{0}$ and $U_{2}$, and most of the excited
states is on potential surface $U_{3}$.

To further illustrate the temperature effects, in Figure 4(c) , we further
plot electron and hole populations as a function of temperature at the steady
states. Note that both electron and hole populations increase with
temperature. This observation is in agreement with experimental results. With
faster mechanical motion and higher pressure at the interfaces, the
temperature of local asperities is larger, such that more electron-hole pairs
are generated at the interfaces. Again, we see the asymmetry of electron and
hole generations at given temperature. We hope future experimental results can
justify our simulation.

(a)

(b)

(c)

Figure 5: Electron and hole populations for different bandgaps ($eV$).
(a)Electron population as a function of time. (b) Hole population as a
function of time. (c) Steady state electron and hole populations as a function
of bandgap. Parameters: $\hbar\omega=0.03,g=0.09,kT=0.15eV$.

We now study the effects of bandgaps on the electron and hole generations. In
Figure 5(a) and Figure 5(b), we plot hole and electron populations at
different bandgaps as a function of time. Notice that the timescales for
electron and hole generations strongly depends on bandgaps. As indicated
before (see Section 3.1), the bandgap strongly affects the hopping rates,
which in turn determines the timesacle for population relaxation. In Figure
5(c), we plot the steady state electron and hole populations as a function of
bandgap. We see that the steady state electron and hole population decrease
with the bandgaps. Obviously, with larger bandgaps, the chance of occupying
excited state decreases, such that there will be less electron and hole
generations. In fact, from our calculation, we see that there is almost no
electron and hole generation for bandgap as large as 3eV. Such a fact will
show their signature in the voltage generation as depicted below.

###  3.3 Voltage generation

With the definition of the electron and hole generations in Eq. 17 and 18, it
is natural to define the open-circuit voltage ($V_{oc}$) as

$\displaystyle V_{oc}=\frac{1}{2}E_{g}(N_{c}+N_{h})$ (19)

Here again, $E_{g}$ is the bandgap. The factor $1/2$ is introduced due to the
fact that we have set Fermi level in the middle of the conduction level and
valence level. Obviously, the voltage depends on the bandgap as well as the
electron and hole generations.

In Figure 6(a), we plot voltage as a function of time for different
temperatures. Note that at longer time, the voltages reach to stationary
values. To further illustrate the effects of the temperature, we plot the
steady state voltages as a function of the temperature in Figure 6(c). Here we
see that the steady state voltage increases with the temperature. This result
is consistent with the observation in Figure 4(c), where both $N_{c}$ and
$N_{h}$ increase with temperature. In general, larger local temperature will
lead to higher voltage generation. Again, this is in agreement with the
experimental results.

(a)

(b)

(c)

(d)

Figure 6: Average voltage ($V_{oc}$) for different local temperatures and
bandgaps. (a) $V_{oc}$ as a function of time for different temperatures,
$E_{g}=1eV$. (b)$V_{oc}$ as a function of time for different bandgaps,
$kT=0.15eV$. (c) Steady state $V_{oc}$ as a function of temperature,
$E_{g}=1eV$. (d) Steady state $V_{oc}$ as a function of bandgap, $kT=0.15eV$.
Other parameters: $\hbar\omega=0.03\ eV,g=0.09\ eV$.

In Figure 6(b), we plot $V_{oc}$ as a function of time for different bandgaps
($E_{g}$). Notice that for smaller bandgaps (< 0.75 eV), the steady state
voltage increase with the bandgap; Whereas for larger bandgaps (> 0.75 eV) the
steady state voltage decreases with bandgap. This is further shown in Figure
6(d), where we clearly see a turnover for the voltage as a function of the
bandgap. The turning point is around 0.75 eV for the bandgap. This turnover is
a result of two competing factors: the electron-hole populations as well as
the bandgaps. See the definition of voltage in Eq. 19. With small bandgap,
there is little voltage generation. Increasing the bandgap will lead to less
electron and hole generations. When the bandgap is large enough, there will be
no electron and hole generations, hence no voltage generation. Such a turnover
is being observed in experiments. We hope to quantify this turnover more
precisely for realistic system in the future study.

##  4 Conclusions

We have offered a quantum mechanic perspective of the tribovaltaic effects at
a sliding metal-semiconductor interfaces. We have introduced a two-level
Anderson-Holstein model to describe nonadiabatic electron and energy transfer
at the interfaces. Furthermore, we have used a classical master equation (CME)
as well as a surface hopping algorithm to model the coupled electron-phonon
dynamics. Using such a dynamical method, we can quantify the electron-hole
pair generations as well as voltage generations. We find that, the electron
and hole populations increase with local temperature, but decrease with the
bandgaps. However, there is a turnover for the voltage generation as a
function of the bandgaps. This finding is strongly supported by the
experimental results. We believe this work present an atomic description of
tribovaltaic. Future work must apply our approach to realistic systems for
better design and control of sliding metal-semiconduction junction.

WD acknowledges the startup funding from Westlake University.

## Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Liu et al. 2018 Liu, J.; Goswami, A.; Jiang, K.; Khan, F.; Kim, S.; McGee, R.; Li, Z.; Hu, Z.; Lee, J.; Thundat, T. Direct-current triboelectricity generation by a sliding Schottky nanocontact on MoS 2 multilayers. _Nature nanotechnology_ 2018, _13_ , 112–116
  * Liu et al. 2019 Liu, J.; Liu, F.; Bao, R.; Jiang, K.; Khan, F.; Li, Z.; Peng, H.; Chen, J.; Alodhayb, A.; Thundat, T. Scaled-up direct-current generation in MoS2 multilayer-based moving heterojunctions. _ACS applied materials & interfaces_ 2019, _11_ , 35404–35409
  * Liu et al. 2018 Liu, J.; Miao, M.; Jiang, K.; Khan, F.; Goswami, A.; McGee, R.; Li, Z.; Nguyen, L.; Hu, Z.; Lee, J. et al. Sustained electron tunneling at unbiased metal-insulator-semiconductor triboelectric contacts. _Nano energy_ 2018, _48_ , 320–326
  * Hao et al. 2019 Hao, Z.; Jiang, T.; Lu, Y.; Feng, S.; Shen, R.; Yao, T.; Yan, Y.; Yang, Y. M.; Lu, Y.; Lin, S. Co-harvesting light and mechanical energy based on dynamic metal/perovskite Schottky junction. _Matter_ 2019, _1_ , 639–649
  * Xu et al. 2019 Xu, R.; Zhang, Q.; Wang, J. Y.; Liu, D.; Wang, J.; Wang, Z. L. Direct current triboelectric cell by sliding an n-type semiconductor on a p-type semiconductor. _Nano Energy_ 2019, _66_ , 104185
  * Zheng et al. 2020 Zheng, M.; Lin, S.; Xu, L.; Zhu, L.; Wang, Z. L. Scanning probing of the tribovoltaic effect at the sliding interface of two semiconductors. _Advanced Materials_ 2020, _32_ , 2000928
  * Lin et al. 2020 Lin, S.; Chen, X.; Wang, Z. L. The tribovoltaic effect and electron transfer at a liquid-semiconductor interface. _Nano Energy_ 2020, _76_ , 105070
  * Lu et al. 2020 Lu, Y.; Gao, Q.; Yu, X.; Zheng, H.; Shen, R.; Hao, Z.; Yan, Y.; Zhang, P.; Wen, Y.; Yang, G. et al. Interfacial Built-In Electric Field-Driven Direct Current Generator Based on Dynamic Silicon Homojunction. _Research_ 2020, _2020_
  * Zheng et al. 2021 Zheng, M.; Lin, S.; Tang, Z.; Feng, Y.; Wang, Z. L. Photovoltaic effect and tribovoltaic effect at liquid-semiconductor interface. _Nano Energy_ 2021, _83_ , 105810
  * Zhang et al. 2020 Zhang, Z.; Jiang, D.; Zhao, J.; Liu, G.; Bu, T.; Zhang, C.; Wang, Z. L. Tribovoltaic Effect on Metal–Semiconductor Interface for Direct-Current Low-Impedance Triboelectric Nanogenerators. _Advanced Energy Materials_ 2020, _10_ , 1903713
  * Liu et al. 2019 Liu, J.; Zhang, Y.; Chen, J.; Bao, R.; Jiang, K.; Khan, F.; Goswami, A.; Li, Z.; Liu, F.; Feng, K. et al. Separation and quantum tunneling of photo-generated carriers using a tribo-induced field. _Matter_ 2019, _1_ , 650–660
  * Zhang et al. 2021 Zhang, Z.; He, T.; Zhao, J.; Liu, G.; Wang, Z.; Zhang, C. Tribo-thermoelectric and tribovoltaic coupling effect at metal-semiconductor interface. _Materials Today Physics_ 2021, _16_ , 100295
  * Yang et al. 2021 Yang, R.; Xu, R.; Dou, W.; Benner, M.; Zhang, Q.; Liu, J. Semiconductor-based dynamic heterojunctions as an emerging strategy for high direct-current mechanical energy harvesting. _Nano Energy_ 2021, 105849
  * Lacks and Shinbrot 2019 Lacks, D. J.; Shinbrot, T. Long-standing and unresolved issues in triboelectric charging. _Nature Reviews Chemistry_ 2019, _3_ , 465–476
  * Lowell and Rose-Innes 1980 Lowell, J.; Rose-Innes, A. Contact electrification. _Advances in Physics_ 1980, _29_ , 947–1023
  * Wang and Wang 2019 Wang, Z. L.; Wang, A. C. On the origin of contact-electrification. _Materials Today_ 2019, _30_ , 34–51
  * Alicki and Jenkins 2020 Alicki, R.; Jenkins, A. Quantum theory of triboelectricity. _Physical Review Letters_ 2020, _125_ , 186101
  * Willatzen et al. 2020 Willatzen, M.; Lew Yan Voon, L. C.; Wang, Z. L. Quantum theory of contact electrification for fluids and solids. _Advanced Functional Materials_ 2020, _30_ , 1910461
  * Antony et al. 2021 Antony, A. C.; Thelen, D.; Zhelev, N.; Adib, K.; Manley, R. G. Electronic charge transfer during metal/SiO2 contact: Insight from density functional theory. _Journal of Applied Physics_ 2021, _129_ , 065304
  * Abdelaziz et al. 2018 Abdelaziz, K. M.; Chen, J.; Hieber, T. J.; Leseman, Z. C. Atomistic Field Theory for contact electrification of dielectrics. _Journal of Electrostatics_ 2018, _96_ , 10–15
  * Park and Salmeron 2014 Park, J. Y.; Salmeron, M. Fundamental aspects of energy dissipation in friction. _Chemical reviews_ 2014, _114_ , 677–711
  * Liu et al. 2019 Liu, J.; Jiang, K.; Nguyen, L.; Li, Z.; Thundat, T. Interfacial friction-induced electronic excitation mechanism for tribo-tunneling current generation. _Materials Horizons_ 2019, _6_ , 1020–1026
  * Dou and Subotnik 2020 Dou, W.; Subotnik, J. E. Nonadiabatic molecular dynamics at metal surfaces. _The Journal of Physical Chemistry A_ 2020, _124_ , 757–771
  * Nienhaus 2002 Nienhaus, H. Electronic excitations by chemical reactions on metal surfaces. _Surface Science Reports_ 2002, _45_ , 1–78
  * Dou et al. 2015 Dou, W.; Nitzan, A.; Subotnik, J. E. Surface hopping with a manifold of electronic states. II. Application to the many-body Anderson-Holstein model. _The Journal of chemical physics_ 2015, _142_ , 084110
  * Dou et al. 2015 Dou, W.; Nitzan, A.; Subotnik, J. E. Surface hopping with a manifold of electronic states. III. Transients, broadening, and the Marcus picture. _The Journal of chemical physics_ 2015, _142_ , 234106
  * Ouyang et al. 2016 Ouyang, W.; Dou, W.; Jain, A.; Subotnik, J. E. Dynamics of barrier crossings for the generalized anderson–holstein model: Beyond electronic friction and conventional surface hopping. _Journal of chemical theory and computation_ 2016, _12_ , 4178–4183
  * Dou et al. 2018 Dou, W.; Schinabeck, C.; Thoss, M.; Subotnik, J. E. A broadened classical master equation approach for treating electron-nuclear coupling in non-equilibrium transport. _The Journal of chemical physics_ 2018, _148_ , 102317
  * Kennedy 2001 Kennedy, F. E. Frictional heating and contact temperatures. _Modern tribology handbook_ 2001, _1_ , 235–272
  * Dou et al. 2015 Dou, W.; Nitzan, A.; Subotnik, J. E. Frictional effects near a metal surface. _The Journal of chemical physics_ 2015, _143_ , 054103