# Quantum walk on simplicial complexes for simplicial community detection

Euijun Song  Present: Independent Researcher, Gyeonggi, Republic of Korea;
<EMAIL_ADDRESS>0000-0002-5886-4210

###### Abstract

Quantum walks have emerged as a transformative paradigm in quantum information
processing and can be applied to various graph problems. This study explores
discrete-time quantum walks on simplicial complexes, a higher-order
generalization of graph structures. Simplicial complexes, encoding higher-
order interactions through simplices, offer a richer topological
representation of complex systems. Leveraging algebraic topology and discrete-
time quantum walk, we present a quantum walk algorithm for detecting higher-
order community structures called simplicial communities. We utilize the
Fourier coin to produce entangled translation states among adjacent simplices
in a simplicial complex. The potential of our quantum algorithm is tested on
Zachary’s karate club network. This study may contribute to understanding
complex systems at the intersection of algebraic topology and quantum
algorithms.

##  1 Introduction

Quantum walks have emerged as a powerful paradigm, offering a unique and
promising avenue for quantum information processing. Quantum walks, inspired
by classical random walks, serve as a cornerstone in the exploration of
quantum algorithms and quantum-enhanced computing [1, 2]. The versatility of
quantum walks extends far beyond algorithmic enhancements. Discrete-time
quantum walks represent a quantum counterpart to the discrete-time classical
random walks, where a quantum particle traverses lattices or graphs in
discrete steps [1]. The discrete nature of the quantum steps introduces an
interplay between quantum coherence and entanglement, offering novel avenues
for information processing beyond the classical random walks. Discrete-time
quantum walks have potential applications in various quantum algorithms for
graph problems, including quantum search [3, 4] and graph community detection
[5]. Recently, quantum walks have been studied on the generalized graph
structures called simplicial complexes [6, 7], which have much more
topological information than graphs or networks.

Network science approaches can extract an interplay between network topology
and dynamics of complex systems. While graphs or networks are based on
pairwise interactions, higher-order interactions between elements often better
capture the topology and dynamics of complex systems [8]. Mathematically,
higher-order interactions can be represented as simplicial complexes [9, 10],
consisting of polytopes called simplices, such as vertices (0-simplices),
edges (1-simplices), triangles (2-simplices), etc. Simplicial homology, an
algebraic topology tool, can extract rich information about topological
features such as high-dimensional holes in simplicial complexes. Using this
algebraic topology method, recent studies have explored higher-order dynamics
[11] and homology-based topological data analysis [12, 13].

In many real networks, the distribution of edges is globally and locally
inhomogeneous. The graph community detection [14, 15] has shed light on how
edges are locally organized within graphs or networks; that is, graphs are
clustered into groups that are densely connected internally, but sparsely
connected between distinct groups. In spectral graph theory, the graph
community structure can be captured by the eigenvectors corresponding to non-
zero eigenvalues of the graph Laplacian [16]. Quantum algorithms have also
been developed for detecting the graph community structure [17, 5]. The graph
community structures can be extended to simplicial complexes, and simplicial
communities are known to be associated with the spectrum of the higher-order
Laplacian, also known as the Hodge Laplacian [18]. This novel concept of the
higher-order community structure has been poorly investigated.

In this study, we present discrete-time quantum walk on simplicial complexes
for detecting simplicial communities, as a higher-order generalization of the
method proposed by Mukai and Hatano [5]. In Section 2, we review the
fundamentals of algebraic topology, such as simplicial complexes and higher-
order Laplacian. In Section 3, the simplicial community and modularity are
mathematically defined, and the symmetry between up-communities and down-
communities is presented. In Section 4, we present the discrete-time quantum
walk method for detecting simplicial communities, based on the Fourier quantum
walk [19, 20]. The long-time average of the transition probability is
calculated to identify simplicial communities. Finally, in Section 5, we test
our quantum algorithm in real networks.

##  2 Simplicial complex and higher-order Laplacian

In this section, we briefly review the key algebraic topology background,
including the concept of simplicial complexes [9, 10] and the Hodge Laplacian
[21].

###  2.1 Simplicial complex

A simplicial complex is a type of higher-order network encoding information
about higher-order interactions. A simplicial complex consists of a set of
polytopes called simplices. An $n$-dimensional simplex ($n$-simplex) $\sigma$
is defined as a set of $n+1$ nodes ($n\geq 0$) with an assigned orientation:

$\sigma=(v_{0},v_{1},\ldots,v_{n}).$ (1)

An $n^{\prime}$-face of an $n$-simplex $\sigma$ ($n^{\prime}<n$) is an
$n^{\prime}$-simplex generated by a subset of the nodes of $\sigma$. A
simplicial complex $\mathcal{K}$ is then defined as a collection of simplices
closed under the inclusion of the faces of simplices. A simplicial $n$-chain
$c_{n}$ is a finite linear combination of $n$-simplices in a simplicial
complex $\mathcal{K}$, as follows:

$c_{n}=\sum_{i}{b_{i}\sigma_{i}}$ (2)

where each $b_{i}$ is an integer coefficient, or generally an element of any
field. A set of $n$-chains on $\mathcal{K}$ forms the free abelian group
$C_{n}$. For $n>0$, the boundary operator $\partial_{n}:C_{n}\to C_{n-1}$ is
defined as the homomorphism that maps $n$-chains to $(n-1)$-chains [9, 10]:

$(v_{0},v_{1},\ldots,v_{n})\mapsto\sum_{i=0}^{n}{(-1)^{i}(v_{0},\ldots,v_{i-1},v_{i+1},\ldots,v_{n})}.$
(3)

One fundamental property of the boundary operator is that the boundary of the
boundary is zero, i.e., $\partial_{n}\partial_{n+1}=0$ (see Hatcher [9]). Let
$N_{n}$ be the number of $n$-simplices in a simplicial complex $\mathcal{K}$.
Using a basis for simplices of $\mathcal{K}$, the boundary operator
$\partial_{n}$ can be represented as an incidence matrix
$B_{n}\in\mathbb{R}^{N_{n-1}\times N_{n}}$, as follows [22]:

$B_{n}(\sigma_{\sim k},\sigma)=(-1)^{k}$ (4)

where $\sigma=(v_{0},v_{1},\ldots,v_{n})$ and $\sigma_{\sim
k}=(v_{0},\ldots,v_{k-1},v_{k+1},\ldots,v_{n})$. Since the boundary of the
boundary is zero, $B_{n}B_{n+1}=0$ holds for $n>0$.

###  2.2 Higher-order Laplacian and adjacency

In graphs or networks, the graph Laplacian is defined as $L=D-A$, where $D$ is
the degree matrix and $A$ is the adjacency matrix of a graph. Using an
incidence matrix $B_{1}$, the graph Laplacian can be written as
$L=B_{1}B_{1}^{T}$. The graph Laplacian is known to be central to the
diffusion process on a graph. The graph Laplacian was generalized to
simplicial complexes by Eckmann [23]. The higher-order Laplacian, also known
as the Hodge Laplacian [21], encodes the higher-order topology and can be used
to examine high-dimensional diffusion on a simplicial complex.

We first define higher-order adjacency for a simplicial complex, as previously
studied [24, 25, 26]. Two $n$-simplices $\sigma$ and $\sigma^{\prime}$ ($n\geq
0$) are upper adjacent if there exists an $(n+1)$-simplex $\gamma$ such that
$\sigma\subset\gamma$ and $\sigma^{\prime}\subset\gamma$; that is, they are
both faces of a common $(n+1)$-simplex. This upper adjacency is denoted by
$\sigma\frown\sigma^{\prime}$. Similarly, two $n$-simplices $\sigma$ and
$\sigma^{\prime}$ ($n>0$) are lower adjacent if there exists an
$(n-1)$-simplex $\gamma$ such that $\gamma\subset\sigma$ and
$\gamma\subset\sigma^{\prime}$; that is they share a common $(n-1)$-face. This
lower adjacency is denoted by $\sigma\smile\sigma^{\prime}$. We also define
the upper degree and lower degree on a simplicial complex. For an $n$-simplex
$\sigma$, the upper degree $d^{u}_{n}(\sigma)$ is the number of
$(n+1)$-simplices of which $\sigma$ is a face. The lower degree
$d^{l}_{n}(\sigma)$ is the number of $(n-1)$-faces in $\sigma$. It is easy to
check $d^{l}_{n}(\sigma)=n+1$ for $n>0$. The upper and lower adjacency
matrices $A^{u}_{n},A^{l}_{n}\in\mathbb{R}^{N_{n}\times N_{n}}$ for
$n$-simplices are defined as follows [24]:

$A^{u}_{n}(\sigma,\sigma^{\prime})=\begin{cases}1&\text{if
}\sigma\frown\sigma^{\prime}\\\ 0&\text{otherwise}\end{cases}$ (5)

and

$A^{l}_{n}(\sigma,\sigma^{\prime})=\begin{cases}1&\text{if
}\sigma\smile\sigma^{\prime}\\\ 0&\text{otherwise.}\end{cases}$ (6)

The Hodge Laplacian describes higher-dimensional diffusion on a simplicial
complex. The diffusion process on a simplicial complex can be performed
through either upper adjacent simplices or lower adjacent simplices. Formally,
the Hodge Laplacian $L_{n}\in\mathbb{R}^{N_{n}\times N_{n}}$ on a simplicial
complex $\mathcal{K}$ is defined as follows [22, 18, 21]:

$L_{n}=\begin{cases}L^{up}_{0}&\text{for }n=0\\\
L^{up}_{n}+L^{down}_{n}&\text{for }n>0\end{cases}$ (7)

where $L^{up}_{n}=B_{n+1}B_{n+1}^{T}$ and $L^{down}_{n}=B_{n}^{T}B_{n}$
represent diffusion through upper adjacent simplicies and lower adjacent
simplices, respectively. For $n>0$ and $n$-simplices
$\sigma,\sigma^{\prime}\in\mathcal{K}$, the matrix elements of $L^{up}_{n}$
and $L^{down}_{n}$ can be written as follows [26]:

$L^{up}_{n}(\sigma,\sigma^{\prime})=\begin{cases}d^{u}_{n}(\sigma)&\text{if
}\sigma=\sigma^{\prime}\\\ -1&\text{if
}\sigma\frown\sigma^{\prime},~{}\Omega(\sigma)=\Omega(\sigma^{\prime})\\\
1&\text{if
}\sigma\frown\sigma^{\prime},~{}\Omega(\sigma)\neq\Omega(\sigma^{\prime})\\\
0&\text{otherwise}\end{cases}$ (8)

and

$L^{down}_{n}(\sigma,\sigma^{\prime})=\begin{cases}d^{l}_{n}(\sigma)=n+1&\text{if
}\sigma=\sigma^{\prime}\\\ 1&\text{if
}\sigma\smile\sigma^{\prime},~{}\Omega(\sigma)=\Omega(\sigma^{\prime})\\\
-1&\text{if
}\sigma\smile\sigma^{\prime},~{}\Omega(\sigma)\neq\Omega(\sigma^{\prime})\\\
0&\text{otherwise}\end{cases}$ (9)

where $\Omega(\sigma)$ indicates the orientation of $\sigma$. If two
$n$-simplices are upper adjacent, then they are also lower adjacent.
Therefore, the off-diagonal element of the Hodge Laplacian $L_{n}$ is non-zero
if and only if $\neg(\sigma\frown\sigma^{\prime})$ and
$\sigma\smile\sigma^{\prime}$. Based on this property of the Hodge Laplacian,
we say that two $n$-simplices $\sigma,\sigma^{\prime}\in\mathcal{K}$ are
adjacent in $\mathcal{K}$ if $\neg(\sigma\frown\sigma^{\prime})$ and
$\sigma\smile\sigma^{\prime}$.

The Hodge Laplacians have useful algebraic properties [22, 18]. By the
definition of the Hodge Laplacians, we obtain $L^{up}_{n}L^{down}_{n}=0$ and
$L^{down}_{n}L^{up}_{n}=0$. This implies that

$\displaystyle\textrm{im}(L^{up}_{n})$ $\displaystyle\subseteq$
$\displaystyle\ker(L^{down}_{n})$ $\displaystyle\textrm{im}(L^{down}_{n})$
$\displaystyle\subseteq$ $\displaystyle\ker(L^{up}_{n}).$ (10)

In addition, the operator
$\mathcal{L}_{n}=\partial_{n+1}\partial_{n+1}^{*}+\partial_{n}^{*}\partial_{n}$
corresponding to the Hodge Laplacian matrix $L_{n}$ has the following property
called the Hodge decomposition [21]:

$C_{n}=\textrm{im}(\partial_{n}^{*})\oplus\ker(\mathcal{L}_{n})\oplus\textrm{im}(\partial_{n+1}).$
(11)

It is known that $\ker(\mathcal{L}_{n})$ is isomorphic to the $n$th homology
group $H_{n}$, implying that the multiplicity of the zero eigenvalues of the
Hodge Laplacian is the same as the number of $n$-dimensional holes in a
simplicial complex.

##  3 Simplicial community

###  3.1 Walks on a simplicial complex

Graph community structures can be generalized to determine higher-order
community structures on a simplicial complex. We first extend the concept of
walks to a simplicial complex [25]. Two $n$-simplices
$\sigma,\sigma^{\prime}\in\mathcal{K}$ are said to be $n$-upper-connected
($n\geq 0$) if there exists a sequence of $n$-simplices
$\sigma=s_{0},s_{1},\ldots,s_{p-1},s_{p}=\sigma^{\prime}$ such that any two
consecutive simplices are upper adjacent, i.e., $s_{i}\frown s_{i+1}$ for all
$0\leq i\leq p-1$. We call this sequence of simplices an $n$-upper-walk.
Similarly, two $n$-simplices $\sigma,\sigma^{\prime}\in\mathcal{K}$ are said
to be $n$-lower-connected ($n>0$) if there exists a sequence of $n$-simplices
$\sigma=s_{0},s_{1},\ldots,s_{p-1},s_{p}=\sigma^{\prime}$ such that any two
consecutive simplices are lower adjacent, i.e., $s_{i}\smile s_{i+1}$ for all
$0\leq i\leq p-1$. We call this sequence of simplices an $n$-lower-walk.

###  3.2 Higher-order community structure

As previously studied by Krishnagopal and Bianconi [18], a simplicial
community can be defined as a set of $n$-simplicies that are $n$-connected.
Specifically, an $n$-up community is a set of $n$-simplicies that are
$n$-upper-connected, and an $n$-down community is a set of $n$-simplicies that
are $n$-lower-connected. We present one important property of the simplicial
community—the symmetry between $n$-up communities and $(n+1)$-down
communities.

Let us consider that $(n+1)$-simplices of a simplicial complex $\mathcal{K}$
are partitioned into $(n+1)$-down communities

$\Pi_{n+1}^{d}=\left\\{\pi_{1}^{d},\ldots,\pi_{g_{n+1}}^{d}\right\\}$ (12)

where $\pi_{i}^{d}$ is the set of $(n+1)$-simplices in the $i$-th $(n+1)$-down
community. Ideally, any two $(n+1)$-simplices belonging to distinct
$(n+1)$-down communities are not $(n+1)$-lower-connected, and
$\pi_{i}^{d}\cap\pi_{j}^{d}=\emptyset$ ($i\neq j$). In a similar way, consider
$n$-up communities of a simplicial complex $\mathcal{K}$:

$\Pi_{n}^{u}=\left\\{\pi_{1}^{u},\ldots,\pi_{h_{n}}^{u}\right\\}.$ (13)

We define
$\hat{\Pi}_{n}^{u}=\left\\{\pi_{i}^{u}\mid|\pi_{i}^{u}|>1,~{}\pi_{i}^{u}\in\Pi_{n}^{u}\right\\}$
by excluding the isolated $n$-simplices. Then, $(n+1)$-down communities
$\Pi_{n+1}^{d}$ is isomorphic to the $n$-up communities (excluding the
isolated $n$-simplices) $\hat{\Pi}_{n}^{u}$. We can prove this symmetry
property by defining the isomorphism $\zeta:\Pi_{n+1}^{d}\to\hat{\Pi}_{n}^{u}$
as follows:

$\pi_{i}^{d}\mapsto\left\\{\sigma_{n}^{\prime}~{}\middle|~{}\sigma_{n}^{\prime}\textrm{
is an }n\textrm{-face of }\sigma_{n+1}\in\pi_{i}^{d}\right\\}.$ (14)

We can also derive this property from the fact that simplicial communities can
be captured by the eigenvectors associated with non-zero eigenvalues of
$L^{up}_{n}=B_{n+1}B_{n+1}^{T}$ (or $L^{down}_{n+1}=B_{n+1}^{T}B_{n+1}$) [18].

Based on the symmetry between up-communities and down-communities, we only
consider $n$-down simplicial communities ($n>0$) in the rest of the present
work. The term simplicial community refers to the $n$-down simplicial
community.

###  3.3 Simplicial modularity

Ideally, any two $n$-simplices belonging to distinct simplicial communities
are not $n$-connected. However, in many complex networks, it is often
necessary to detect community structures that are densely connected internally
but sparsely connected between distinct communities. In graphs or networks,
algorithms have been developed for maximizing modularity, which is defined as
the difference between the fraction of the edges within the given communities
and the expected fraction at random distribution [15]. We can extend the
concept of modularity to the $n$-down simplicial communities
$\Pi_{n}^{d}=\left\\{\pi_{1}^{d},\ldots,\pi_{g_{n}}^{d}\right\\}$. Define a
matrix $W_{n}\in\mathbb{R}^{N_{n}\times|\Pi_{n}^{d}|}$, by
$W_{n}(\sigma,\pi_{i}^{d})$ to be 1 if $\sigma\in\pi_{i}^{d}$ and otherwise
zero. We denote the lower neighborhood, a set of lower-adjacent simplices of
an $n$-simplex $\sigma$, by

$\mathcal{N}^{l}(\sigma)=\left\\{\sigma^{\prime}~{}\middle|~{}\sigma\smile\sigma^{\prime}\right\\}.$
(15)

We then define the simplicial modularity for the $n$-down simplicial
communities as follows:

$Q_{n}=\frac{1}{m_{n}}\textrm{Tr}\left(W^{T}_{n}M_{n}W_{n}\right)$ (16)

where $m_{n}=\sum_{\sigma}\left|\mathcal{N}^{l}(\sigma)\right|$ and
$M_{n}\in\mathbb{R}^{N_{n}\times N_{n}}$ is the modularity matrix:

$M_{n}(\sigma,\sigma^{\prime})=A^{l}_{n}(\sigma,\sigma^{\prime})-\frac{\left|\mathcal{N}^{l}(\sigma)\right|\cdot\left|\mathcal{N}^{l}(\sigma^{\prime})\right|}{m_{n}}.$
(17)

Although modularity maximization is one of the most popular approaches for
detecting community structures in graphs, it has a resolution limit leading to
merging small clusters and it has to determine the community number parameter.
In the next section, we present the quantum walk algorithm for detecting
simplicial communities.

##  4 Quantum walk for simplicial community detection

###  4.1 Discrete-time quantum walk on simplicial complexes

Here, we present the discrete-time quantum walk on simplicial complexes for
detecting $n$-down simplicial communities ($n>0$). We directly generalize the
community detection method proposed by Mukai and Hatano [5] to the simplicial
community. The method utilizes the Fourier coin (or the Fourier quantum walk),
which has been studied in high-dimensional lattices [19], regular graphs [20],
and complex networks [5].

Let $\sigma_{1},\sigma_{2},\ldots,\sigma_{N_{n}}\in\mathcal{K}$ be
$n$-simplices in a simplicial complex $\mathcal{K}$. The total Hilbert space
is given by $\mathcal{H}=\bigoplus_{i=1}^{N_{n}}\mathcal{H}_{\sigma_{i}}$,
consisting of the Hilbert space corresponding to each $\sigma_{i}$. Each
Hilbert space $\mathcal{H}_{\sigma_{i}}$ is spanned by the orthonormal basis
$\left\\{|\sigma_{i}\to\sigma_{j}\rangle\mid\sigma_{j}\in\mathcal{N}^{l}(\sigma_{i})\right\\}$,
where $|\sigma_{i}\to\sigma_{j}\rangle$ represents the translation state
jumping from $\sigma_{i}$ to its lower-adjacent simplex
$\sigma_{j}\in\mathcal{N}^{l}(\sigma_{i})$. The dimension of the total Hilbert
space $\mathcal{H}$ is equal to
$m_{n}=\sum_{\sigma_{i}}\left|\mathcal{N}^{l}(\sigma_{i})\right|$. The quantum
state $|\Psi(t)\rangle\in\mathcal{H}$ is given by

$\displaystyle|\Psi(t)\rangle$ $\displaystyle=$
$\displaystyle\sum_{\sigma_{i}\smile\sigma_{j}}\psi_{\sigma_{i},\sigma_{j}}(t)|\sigma_{i}\to\sigma_{j}\rangle$
(18) $\displaystyle=$
$\displaystyle\sum_{i=1}^{N_{n}}\sum_{\sigma_{j}\in\mathcal{N}^{l}(\sigma_{i})}\psi_{\sigma_{i},\sigma_{j}}(t)|\sigma_{i}\to\sigma_{j}\rangle.$

The sum of probabilities is equal to 1, as
$\langle\Psi(t)|\Psi(t)\rangle=\sum_{\sigma_{i}\smile\sigma_{j}}\left|\psi_{\sigma_{i},\sigma_{j}}(t)\right|^{2}=1$.

We now define unitary time evolution operator $U:\mathcal{H}\to\mathcal{H}$ to
specify how our quantum walk works. The time evolution of the quantum state is
given by $|\Psi(t)\rangle=U^{t}|\Psi(0)\rangle$. The unitary operator is
basically given by $U=SC$, where $S$ is the shift operator and $C$ is the coin
operator. We define the coin operator $C:\mathcal{H}\to\mathcal{H}$ as a
direct sum of the Fourier quantum walks [19, 20]
$F_{\sigma_{i}}:\mathcal{H}_{\sigma_{i}}\to\mathcal{H}_{\sigma_{i}}$ as
follows:

$C=\bigoplus_{i=1}^{N_{n}}F_{\sigma_{i}}$ (19)

and

$F_{\sigma_{i}}=\frac{1}{\sqrt{\left|\mathcal{N}^{l}(\sigma_{i})\right|}}\sum_{\begin{subarray}{c}\sigma_{j_{\alpha}}\in\mathcal{N}^{l}(\sigma_{i})\\\
\alpha=1,2,\ldots,|\mathcal{N}^{l}(\sigma_{i})|\end{subarray}}~{}\sum_{\begin{subarray}{c}\sigma_{j_{\beta}}\in\mathcal{N}^{l}(\sigma_{i})\\\
\beta=1,2,\ldots,|\mathcal{N}^{l}(\sigma_{i})|\end{subarray}}\omega^{\alpha\beta}|\sigma_{i}\to\sigma_{j_{\beta}}\rangle\langle\sigma_{i}\to\sigma_{j_{\alpha}}|$
(20)

where $\omega=e^{\frac{2\pi\sqrt{-1}}{|\mathcal{N}^{l}(\sigma_{i})|}}$. This
Fourier coin operator $F_{\sigma_{i}}$ produces entanglement by transforming
any translation state $|\sigma_{i}\to\sigma_{j_{\alpha}}\rangle$ into an
equally-weighted superposition of all the translation states
$|\sigma_{i}\to\sigma_{j_{\beta}}\rangle$ based on $\sigma_{i}$. We define the
shift operator $S:\mathcal{H}\to\mathcal{H}$ as follows:

$S=\sum_{\sigma_{i}\smile\sigma_{j}}|\sigma_{j}\to\sigma_{i}\rangle\langle\sigma_{i}\to\sigma_{j}|.$
(21)

For lower-adjacent simplices $\sigma_{i}$ and $\sigma_{j}$, this shift
operator $S$ maps $|\sigma_{i}\to\sigma_{j}\rangle$ to
$|\sigma_{j}\to\sigma_{i}\rangle$, i.e., moving the walker from $\sigma_{i}$
to $\sigma_{j}\in\mathcal{N}^{l}(\sigma_{i})$, and vice versa. It can be
represented as a direct sum of $2\times 2$ permutation matrices:

$S=\bigoplus_{\begin{subarray}{c}\sigma_{i}\smile\sigma_{j}\\\
i<j\end{subarray}}\begin{bmatrix}0&1\\\ 1&0\end{bmatrix}.$ (22)

###  4.2 Transition probability

We calculate the transition probability from one $n$-simplex $\sigma_{x}$ to
another $n$-simplex $\sigma_{y}$ (not necessarily lower adjacent). We set the
initial quantum state $|\Psi(0)\rangle$ such that the quantum walk starts from
$|\sigma_{x}\to\sigma_{v}\rangle$ ($\sigma_{v}\in\mathcal{N}^{l}(\sigma_{x})$)
and calculate the average probability over $\mathcal{N}^{l}(\sigma_{x})$. The
normalized transition probability that the quantum walker starting from
$\sigma_{x}$ reaches $\sigma_{y}$ at time $t$ is given by

$p(\sigma_{x}\to\sigma_{y};t)=\frac{1}{\left|\mathcal{N}^{l}(\sigma_{x})\right|\cdot\left|\mathcal{N}^{l}(\sigma_{y})\right|}\sum_{\sigma_{w}\in\mathcal{N}^{l}(\sigma_{y})}\sum_{\sigma_{v}\in\mathcal{N}^{l}(\sigma_{x})}\left|\langle\sigma_{y}\to\sigma_{w}|U^{t}|\sigma_{x}\to\sigma_{v}\rangle\right|^{2}$
(23)

where $|\sigma_{x}\to\sigma_{v}\rangle$ is the initial state and
$|\sigma_{y}\to\sigma_{w}\rangle$ is the final state at time $t$.

Let us calculate the long-time average of the normalized transition
probability:

$q(\sigma_{x}\to\sigma_{y})=\lim_{T\to\infty}\frac{1}{T}\sum_{t=1}^{T}p(\sigma_{x}\to\sigma_{y};t).$
(24)

In the practical situation, $q(\sigma_{x}\to\sigma_{y})$ can be approximated
by finite-time results for the quantum walk. Here, we theoretically calculate
$q(\sigma_{x}\to\sigma_{y})$ using the spectral decomposition of the unitary
operator $U$. Let $e^{i\theta_{k}}$ ($i=\sqrt{-1}$) and $|\Phi_{k}\rangle$ be
the eigenvalue and the corresponding eigenstate of the unitary operator $U$,
respectively ($\theta_{k}\in[0,2\pi)$, $k=1,2,\ldots,m_{n}$). Suppose the
eigenvalues of $U$ are non-degenerate. The spectral decomposition of $U$ is
given by $U=\sum_{k=1}^{m_{n}}|\Phi_{k}\rangle
e^{i\theta_{k}}\langle\Phi_{k}|$. Hence,

$\displaystyle\lim_{T\to\infty}\frac{1}{T}\sum_{t=1}^{T}\left|\langle\sigma_{y}\to\sigma_{w}|U^{t}|\sigma_{x}\to\sigma_{v}\rangle\right|^{2}$
$\displaystyle=\lim_{T\to\infty}\frac{1}{T}\sum_{t=1}^{T}\left|\sum_{k=1}^{m_{n}}e^{i\theta_{k}t}\langle\sigma_{y}\to\sigma_{w}|\Phi_{k}\rangle\langle\Phi_{k}|\sigma_{x}\to\sigma_{v}\rangle\right|^{2}$
$\displaystyle=\lim_{T\to\infty}\frac{1}{T}\sum_{t=1}^{T}\sum_{k_{1},k_{2}}e^{i(\theta_{k_{1}}-\theta_{k_{2}})t}\langle\sigma_{y}\to\sigma_{w}|\Phi_{k_{1}}\rangle\langle\Phi_{k_{1}}|\sigma_{x}\to\sigma_{v}\rangle\langle\sigma_{x}\to\sigma_{v}|\Phi_{k_{2}}\rangle\langle\Phi_{k_{2}}|\sigma_{y}\to\sigma_{w}\rangle$
$\displaystyle=\sum_{k_{1},k_{2}}\delta_{\theta_{k_{1}},\theta_{k_{2}}}\langle\sigma_{y}\to\sigma_{w}|\Phi_{k_{1}}\rangle\langle\Phi_{k_{1}}|\sigma_{x}\to\sigma_{v}\rangle\langle\sigma_{x}\to\sigma_{v}|\Phi_{k_{2}}\rangle\langle\Phi_{k_{2}}|\sigma_{y}\to\sigma_{w}\rangle$
$\displaystyle=\sum_{k=1}^{m_{n}}\left|\langle\sigma_{y}\to\sigma_{w}|\Phi_{k}\rangle\langle\Phi_{k}|\sigma_{x}\to\sigma_{v}\rangle\right|^{2}$
(25)

where $\delta_{\theta_{k_{1}},\theta_{k_{2}}}$ is the Kronecker delta, and we
used the following property (see Appendix A for proof):

$\lim_{T\to\infty}\frac{1}{T}\sum_{t=1}^{T}e^{i(\theta_{k_{1}}-\theta_{k_{2}})t}=\delta_{\theta_{k_{1}},\theta_{k_{2}}}.$
(26)

By combining Eq. 23 and Eq. 4.2, we have

$q(\sigma_{x}\to\sigma_{y})=\frac{1}{\left|\mathcal{N}^{l}(\sigma_{x})\right|\cdot\left|\mathcal{N}^{l}(\sigma_{y})\right|}\sum_{\sigma_{w}\in\mathcal{N}^{l}(\sigma_{y})}\sum_{\sigma_{v}\in\mathcal{N}^{l}(\sigma_{x})}\sum_{k=1}^{m_{n}}\left|\langle\sigma_{y}\to\sigma_{w}|\Phi_{k}\rangle\langle\Phi_{k}|\sigma_{x}\to\sigma_{v}\rangle\right|^{2}.$
(27)

We can estimate the lower bound of $q(\sigma_{x}\to\sigma_{y})$ by the average
transition amplitude, as follows:

$\displaystyle q(\sigma_{x}\to\sigma_{y})$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{m_{n}}\frac{1}{\left|\mathcal{N}^{l}(\sigma_{x})\right|\cdot\left|\mathcal{N}^{l}(\sigma_{y})\right|}\sum_{\sigma_{w}\in\mathcal{N}^{l}(\sigma_{y})}\sum_{\sigma_{v}\in\mathcal{N}^{l}(\sigma_{x})}\left|\langle\sigma_{y}\to\sigma_{w}|\Phi_{k}\rangle\langle\Phi_{k}|\sigma_{x}\to\sigma_{v}\rangle\right|^{2}$
(28) $\displaystyle\geq$
$\displaystyle\sum_{k=1}^{m_{n}}\left|\frac{1}{\left|\mathcal{N}^{l}(\sigma_{x})\right|\cdot\left|\mathcal{N}^{l}(\sigma_{y})\right|}\sum_{\sigma_{w}\in\mathcal{N}^{l}(\sigma_{y})}\sum_{\sigma_{v}\in\mathcal{N}^{l}(\sigma_{x})}\langle\sigma_{y}\to\sigma_{w}|\Phi_{k}\rangle\langle\Phi_{k}|\sigma_{x}\to\sigma_{v}\rangle\right|^{2}$
$\displaystyle=$
$\displaystyle\sum_{k=1}^{m_{n}}\left|\langle\sigma_{y}|\Phi_{k}\rangle\langle\Phi_{k}|\sigma_{x}\rangle\right|^{2}$

where $|\sigma_{x}\rangle$ and $|\sigma_{y}\rangle$ are equally-weighted
superpositions of all the translation states based on $\sigma_{x}$ and
$\sigma_{y}$, respectively:

$\displaystyle|\sigma_{x}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\left|\mathcal{N}^{l}(\sigma_{x})\right|}\sum_{\sigma_{v}\in\mathcal{N}^{l}(\sigma_{x})}|\sigma_{x}\to\sigma_{v}\rangle$
$\displaystyle|\sigma_{y}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\left|\mathcal{N}^{l}(\sigma_{y})\right|}\sum_{\sigma_{w}\in\mathcal{N}^{l}(\sigma_{y})}|\sigma_{y}\to\sigma_{w}\rangle.$
(29)

###  4.3 Simplicial community detection

Since the Fourier quantum walk on graphs is likely to be localized in a
community containing the initial node [5], we expect our generalized framework
could work for detecting simplicial communities. We present the quantum walk
algorithm for detecting simplicial communities based on the long-time average
of the normalized transition probability $q(\sigma_{x}\to\sigma_{y})$ (Eq.
24).

  1. 1.

Set a starting $n$-simplex $\sigma_{x}\in\mathcal{K}$ whose
$\left|\mathcal{N}^{l}(\sigma_{x})\right|$ is the maximum over all the
unassigned $n$-simplices.

  2. 2.

Compute $q(\sigma_{x}\to\sigma_{y})$ for each of the unassigned $n$-simplices
$\sigma_{y}\in\mathcal{K}$.

  3. 3.

Assign an $n$-simplex $\sigma_{y}$ to a member of the simplicial community of
$\sigma_{x}$, if

$q(\sigma_{x}\to\sigma_{y})>\frac{1}{m_{n}}.$ (30)

  4. 4.

Repeat until all the $n$-simplices are assigned to the simplicial communities.

The threshold (Eq. 30) is set to be the stationary probability of the
classical random walk $\frac{1}{m_{n}}$, so that the transition probability
that the quantum walker reaches the same simplicial community of the starting
simplex is higher than that expected by chance.

We can theoretically calculate the long-time average
$q(\sigma_{x}\to\sigma_{y})$ using the spectrum of the unitary operator $U$,
as shown in Eq. 27. Since the spectral decomposition requires a large
computational complexity of $O(\textrm{poly }m_{n})$, the finite-time
approximation by repeatedly applying the unitary operator $U$ is a more
efficient and natural way in quantum computation. We therefore approximate
$q(\sigma_{x}\to\sigma_{y})$ by the finite-time results for the quantum walk.
In our simulation analysis, we compute the finite-time average of the
normalized transition probability over 100 time steps from $t=1$ to $t=100$,
and denote this finite-time approximation by
$\tilde{q}(\sigma_{x}\to\sigma_{y})$.

##  5 Application to real networks

We test our quantum algorithm to Zachary’s karate club network [27],
consisting of 34 nodes and 78 edges. According to Zachary’s study [27], the
nodes are clustered into two communities: the Mr. Hi community around node 1
and the Officer community around node 34. We detect simplicial communities of
Zachary’s karate club network using our quantum algorithm. As shown in Figure
1, $N_{1}=78$ edges (1-simplices) of Zachary’s karate club network are
assigned to two 1-down communities. Most elements of the first 1-down
community (red edges in Fig. 1) consist of the edges connected within the Mr.
Hi group, whereas most elements of the second 1-down community (blue edges in
Fig. 1) consist of the edges connected within the Officer group. There are
several edges (1-simplices) connecting the Officer group and the Mr. Hi group:
(1, 32), (2, 31), (3, 10), (3, 28), and (3, 29) in the first 1-down community,
and (3, 33), (9, 31), (9, 33), (9, 34), (14, 34), and (20, 34) in the second
1-down community. The simplicial modularity (Eq. 16) for the 1-down
communities is $Q_{1}=0.434$.

Figure 1: The 1-down communities of Zachary’s karate club network.

We further detect higher-order simplicial communities of Zachary’s karate club
network. As highlighted in Table 1, $N_{2}=45$ triangles (2-simplices) are
clustered into four 2-down communities, one of which is an isolated 2-simplex,
(25, 26, 32). The simplicial modularity for the 2-down communities is
$Q_{2}=0.515$. There are also three 3-down communities, two of which are
isolated 3-simplices, (9, 31, 33, 34) and (24, 30, 33, 34). There is only one
4-down community consisting of (1, 2, 3, 4, 8) and (1, 2, 3, 4, 14) in the Mr.
Hi group. The simplicial modularity values for the 3-down communities and
4-down community are zero because all the simplicial communities except one
are isolated simplices.

Table 1: Simplicial communities of Zachary’s karate club network. $n$ | $N_{n}$ |  $n$-down simplicial communities
---|---|---
1 | 78 |  (see Fig. 1)
2 | 45 |  • (1, 2, 3), (1, 2, 4), (1, 2, 8), (1, 2, 14), (1, 2, 18), (1, 2, 20), (1, 2, 22), (1, 3, 4), (1, 3, 8), (1, 3, 9), (1, 3, 14), (1, 4, 8), (1, 4, 13), (1, 4, 14), (2, 3, 4), (2, 3, 8), (2, 3, 14), (2, 4, 8), (2, 4, 14), (3, 4, 8), (3, 4, 14) • (3, 9, 33), (9, 31, 33), (9, 31, 34), (9, 33, 34), (15, 33, 34), (16, 33, 34), (19, 33, 34), (21, 33, 34), (23, 33, 34), (24, 28, 34), (24, 30, 33), (24, 30, 34), (24, 33, 34), (27, 30, 34), (29, 32, 34), (30, 33, 34), (31, 33, 34), (32, 33, 34) • (1, 5, 7), (1, 5, 11), (1, 6, 7), (1, 6, 11), (6, 7, 17) • (25, 26, 32)
3 | 11 |  • (1, 2, 3, 4), (1, 2, 3, 8), (1, 2, 3, 14), (1, 2, 4, 8), (1, 2, 4, 14), (1, 3, 4, 8), (1, 3, 4, 14), (2, 3, 4, 8), (2, 3, 4, 14) • (9, 31, 33, 34) • (24, 30, 33, 34)
4 | 2 |  • (1, 2, 3, 4, 8), (1, 2, 3, 4, 14)

##  6 Conclusion

This study investigates discrete-time quantum walks on simplicial complexes,
which represent a higher-order extension of graphs or networks. By leveraging
algebraic topology and discrete-time quantum walk, we present a quantum walk
algorithm for the identification of higher-order community structures called
simplicial communities. The symmetry between up-communities and down-
communities is presented. Our approach utilizes the Fourier coin to generate
entangled translation states among neighboring simplices in a simplicial
complex. Our quantum algorithm is tested in Zachary’s karate club network as a
test case; however, it needs to be validated through various examples and
mathematical proof. This research has the potential to contribute to the
understanding of complex systems at the crossroads of algebraic topology and
quantum algorithms.

## Acknowledgement

This independent research received no external funding. The author has no
conflicts of interest to declare. The author would like to thank anonymous
reviewers for their valuable comments.

## Authorship contribution

Euijun Song: Conceptualization, Methodology, Formal analysis, Investigation,
Visualization, Writing – original draft.

##  Appendix A Appendix

We prove the following property of the long-time average (Eq. 26):

$\lim_{T\to\infty}\frac{1}{T}\sum_{t=1}^{T}e^{i(\theta_{k_{1}}-\theta_{k_{2}})t}=\delta_{\theta_{k_{1}},\theta_{k_{2}}}$

where $i=\sqrt{-1}$ and $\theta_{k_{1}},\theta_{k_{2}}\in[0,2\pi)$.

###### Proof.

(i) If $\theta_{k_{1}}=\theta_{k_{2}}$,

$\lim_{T\to\infty}\frac{1}{T}\sum_{t=1}^{T}e^{i(\theta_{k_{1}}-\theta_{k_{2}})t}=\lim_{T\to\infty}\frac{1}{T}\sum_{t=1}^{T}1=1.$

(ii) If $\theta_{k_{1}}\neq\theta_{k_{2}}$,

$\displaystyle\lim_{T\to\infty}\frac{1}{T}\sum_{t=1}^{T}e^{i(\theta_{k_{1}}-\theta_{k_{2}})t}$
$\displaystyle=$
$\displaystyle\lim_{T\to\infty}\frac{1}{T}\cdot\frac{e^{i(\theta_{k_{1}}-\theta_{k_{2}})}\left(1-e^{i(\theta_{k_{1}}-\theta_{k_{2}})T}\right)}{1-e^{i(\theta_{k_{1}}-\theta_{k_{2}})}}$
$\displaystyle=$
$\displaystyle\frac{e^{i(\theta_{k_{1}}-\theta_{k_{2}})}}{1-e^{i(\theta_{k_{1}}-\theta_{k_{2}})}}\lim_{T\to\infty}\frac{1-e^{i(\theta_{k_{1}}-\theta_{k_{2}})T}}{T}$
$\displaystyle=$ $\displaystyle 0.$

The last line holds because
$\left|1-e^{i(\theta_{k_{1}}-\theta_{k_{2}})T}\right|\leq 2$. ∎

## References

  * Aharonov et al. [1993] Y. Aharonov, L. Davidovich, and N. Zagury.  Quantum random walks.  _Phys. Rev. A_ , 48:1687–1690, 1993.  doi: 10.1103/PhysRevA.48.1687.
  * Travaglione and Milburn [2002] B. C. Travaglione and G. J. Milburn.  Implementing the quantum random walk.  _Phys. Rev. A_ , 65:032310, 2002.  doi: 10.1103/PhysRevA.65.032310.
  * Berry and Wang [2010] Scott D. Berry and Jingbo B. Wang.  Quantum-walk-based search and centrality.  _Phys. Rev. A_ , 82:042333, 2010.  doi: 10.1103/PhysRevA.82.042333.
  * Grover [1996] Lov K. Grover.  A fast quantum mechanical algorithm for database search.  In _Proceedings of the Twenty-Eighth Annual ACM Symposium on Theory of Computing_ , STOC ’96, pages 212–219. Association for Computing Machinery, 1996.  doi: 10.1145/237814.237866.
  * Mukai and Hatano [2020] Kanae Mukai and Naomichi Hatano.  Discrete-time quantum walk on complex networks for community detection.  _Phys. Rev. Res._ , 2:023378, 2020.  doi: 10.1103/PhysRevResearch.2.023378.
  * Matsue et al. [2016] Kaname Matsue, Osamu Ogurisu, and Etsuo Segawa.  Quantum walks on simplicial complexes.  _Quantum Information Processing_ , 15(5):1865–1896, 2016.  doi: 10.1007/s11128-016-1247-6.
  * Matsue et al. [2018] Kaname Matsue, Osamu Ogurisu, and Etsuo Segawa.  Quantum search on simplicial complexes.  _Quantum Studies: Mathematics and Foundations_ , 5(4):551–577, 2018.  doi: 10.1007/s40509-017-0144-8.
  * Giusti et al. [2016] Chad Giusti, Robert Ghrist, and Danielle S. Bassett.  Two’s company, three (or more) is a simplex.  _Journal of Computational Neuroscience_ , 41(1):1–14, 2016.  doi: 10.1007/s10827-016-0608-6.
  * Hatcher [2002] Allen Hatcher.  _Algebraic topology_.  Cambridge University Press, 2002.
  * Jonsson [2008] Jakob Jonsson.  _Simplicial complexes of graphs_ , volume 3.  Springer, 2008.
  * Millán et al. [2020] Ana P. Millán, Joaquín J. Torres, and Ginestra Bianconi.  Explosive higher-order kuramoto dynamics on simplicial complexes.  _Phys. Rev. Lett._ , 124:218301, 2020.  doi: 10.1103/PhysRevLett.124.218301.
  * Song [2023] Euijun Song.  Persistent homology analysis of type 2 diabetes genome-wide association studies in protein–protein interaction networks.  _Frontiers in Genetics_ , 14:1270185, 2023.  doi: 10.3389/fgene.2023.1270185.
  * Vipond et al. [2021] Oliver Vipond, Joshua A. Bull, Philip S. Macklin, Ulrike Tillmann, Christopher W. Pugh, Helen M. Byrne, and Heather A. Harrington.  Multiparameter persistent homology landscapes identify immune cell spatial patterns in tumors.  _Proceedings of the National Academy of Sciences_ , 118(41):e2102166118, 2021.  doi: 10.1073/pnas.2102166118.
  * Girvan and Newman [2002] M. Girvan and M. E. J. Newman.  Community structure in social and biological networks.  _Proceedings of the National Academy of Sciences_ , 99(12):7821–7826, 2002.  doi: 10.1073/pnas.122653799.
  * Newman [2006] M. E. J. Newman.  Modularity and community structure in networks.  _Proceedings of the National Academy of Sciences_ , 103(23):8577–8582, 2006.  doi: 10.1073/pnas.0601602103.
  * Newman [2013] M. E. J. Newman.  Spectral methods for community detection and graph partitioning.  _Phys. Rev. E_ , 88:042822, 2013.  doi: 10.1103/PhysRevE.88.042822.
  * Faccin et al. [2014] Mauro Faccin, Piotr Migdał, Tomi H. Johnson, Ville Bergholm, and Jacob D. Biamonte.  Community detection in quantum complex networks.  _Phys. Rev. X_ , 4:041012, 2014.  doi: 10.1103/PhysRevX.4.041012.
  * Krishnagopal and Bianconi [2021] Sanjukta Krishnagopal and Ginestra Bianconi.  Spectral detection of simplicial communities via hodge laplacians.  _Phys. Rev. E_ , 104:064303, 2021.  doi: 10.1103/PhysRevE.104.064303.
  * Mackay et al. [2002] T. D. Mackay, S. D. Bartlett, L. T. Stephenson, and B. C. Sanders.  Quantum walks in higher dimensions.  _Journal of Physics A: Mathematical and General_ , 35(12):2745, 2002.  doi: 10.1088/0305-4470/35/12/304.
  * Saito [2019] Kei Saito.  Periodicity for the fourier quantum walk on regular graphs.  _Quantum Info. Comput._ , 19(1-2):23–34, 2019\.  doi: 10.26421/QIC19.1-2-3.
  * Lim [2020] Lek-Heng Lim.  Hodge laplacians on graphs.  _SIAM Review_ , 62(3):685–715, 2020.  doi: 10.1137/18M1223101.
  * Baccini et al. [2022] Federica Baccini, Filippo Geraci, and Ginestra Bianconi.  Weighted simplicial complexes and their representation power of higher-order network data and topology.  _Phys. Rev. E_ , 106:034319, 2022.  doi: 10.1103/PhysRevE.106.034319.
  * Eckmann [1944] Beno Eckmann.  Harmonische funktionen und randwertaufgaben in einem komplex.  _Commentarii Mathematici Helvetici_ , 17(1):240–255, 1944.
  * Estrada and Ross [2018] Ernesto Estrada and Grant J. Ross.  Centralities in simplicial complexes. applications to protein interaction networks.  _Journal of Theoretical Biology_ , 438:46–60, 2018.  doi: 10.1016/j.jtbi.2017.11.003.
  * Hernández Serrano and Sánchez Gómez [2020] Daniel Hernández Serrano and Darío Sánchez Gómez.  Centrality measures in simplicial complexes: Applications of topological data analysis to network science.  _Applied Mathematics and Computation_ , 382:125331, 2020\.  doi: 10.1016/j.amc.2020.125331.
  * Muhammad and Egerstedt [2006] Abubakr Muhammad and Magnus Egerstedt.  Control using higher order laplacians in network topologies.  In _Proc. of 17th International Symposium on Mathematical Theory of Networks and Systems_ , pages 1024–1038, 2006.
  * Zachary [1977] Wayne W. Zachary.  An information flow model for conflict and fission in small groups.  _Journal of Anthropological Research_ , 33(4):452–473, 1977.  doi: 10.1086/jar.33.4.3629752.


# Existence of steady states of the massless Einstein-Vlasov system
surrounding a Schwarzschild black hole

Håkan Andréasson
Mathematical Sciences
Chalmers University of Technology
University of Gothenburg
SE-41296 Gothenburg, Sweden
email<EMAIL_ADDRESS>
###### Abstract

We show that there exist steady states of the massless Einstein-Vlasov system
which surround a Schwarzschild black hole. The steady states are (thick)
shells with finite mass and compact support. Furthermore we prove that an
arbitrary number of shells, necessarily well separated, can surround the black
hole. To our knowledge this is the first result of static self-gravitating
solutions to any massless Einstein-matter system which surround a black hole.
We also include a numerical investigation about the properties of the shells.

##  1 Introduction

The Einstein-Vlasov system typically models self-gravitating particle
ensembles such as galaxies or clusters of galaxies. The particles in the
former case are stars and in the latter case they are galaxies. Clearly, the
particles carry mass in these two situations. In this work we are instead
interested in the case of massless particles, e.g. photons, and we show that
there exist self-gravitating ensembles of massless particles with finite mass
and compact support surrounding a Schwarzschild black hole. To put our result
in context let us briefly review some related results. Existence of steady
states to the Einstein-Vlasov system in the case of massive particles was
first established in [23]. The steady states constructed in this work are
spherically symmetric with a regular centre. Several simplifications and
generalizations have since then been obtained and we refer to [20] for a
simplified and general approach, to [6] for the existence of highly
relativistic static solutions and to [11] for the existence of stationary
solutions in the axisymmetric case. There are several other existence results
and also results about the properties of the static solutions in the
literature, and we refer to [8] for a review and to [9, 2, 3] for more recent
results.

By relaxing the condition of a regular centre the case with a Schwarzschild
black hole was considered in [21], where the existence of massive static
shells of Vlasov matter surrounding a black hole was shown. A different method
leading to a similar result was more recently given by Jabiri in [18]. For a
fluid, the first result of a massive static shell surrounding a black hole was
obtained in [15]. If the matter model originates from quantum mechanics
similar results need not be true as is for instance shown in [16], where the
absence of static black hole solutions is shown for the Einstein-Dirac-
Yang/Mills equations. It is argued in [16] that a reason for the difference
between the classical and the quantum mechanical case is that classical
particles are prevented from falling into the black hole by the centrifugal
barrier, whereas quantum particles can tunnel through this barrier.

Solutions of the Einstein-Vlasov system can also model ensembles of massless
particles, e.g. photons. The first mathematical study of the massless
Einstein-Vlasov system is to our knowledge the work [24] by Rendall, where the
dynamics of cosmological solutions is investigated. Only more recently results
about static solutions have been obtained. Akbarian and Choptuik constructed
massless solutions with compact support numerically in [1]. An existence proof
was obtained in [10], where also a discussion about the relation to Wheeler’s
concept of geons is given. Gundlach studied the problem by numeric and
analytic tools in [17]. An important difference between the massive and
massless case is that the existence of massless static solutions requires that
the solutions are highly relativistic in the sense that the compactness ratio
$2M/R$ is large. Here $M$ is the ADM mass of the solution and $R$ its (areal)
radial support. It is known that $2M/R$ is always bounded by $8/9$, cf. [7].
(The classical result by Buchdahl [14] does not apply in this case although
the bound is the same.) Numerically it has been found that a necessary lower
bound is roughly $2M/R>4/5$, cf. [1, 10, 17], for the existence of massless
static solutions whereas no such lower bound is needed in the massive case.

In the present work we combine the methods from [21] and [10]. We consider the
case with a Schwarzschild black hole in the centre and we show that there
exist static massless shells of Vlasov matter with compact support and finite
mass which surround the black hole. Necessarily there is a gap between the
black hole and the shell; the inner radius of the shell has to be larger than
the radius of the photon sphere of the black hole. In our proof the shell is
placed far away from the photon sphere. This is a technical condition.
Numerically, we find that there are situations when the shell can be arbitrary
close to the photon sphere, cf. Section 5. The shell solutions are highly
relativistic in the sense that $2M/R$ is large. However, when the shell can be
placed close to the photon sphere the ratio $2M/R$ is larger than, but close
to, $2/3$. Hence, the presence of a black hole reduces the required lower
bound of $2M/R$. Clearly, since the ratio $2M/R$ of the shell is larger than
$2/3$ there is a photon sphere surrounding the shell in addition to the photon
sphere which surrounds the black hole and which is situated between the black
hole and the shell. Our result can be generalized to the case of an arbitrary
number of shells. The resulting spacetime thus contains an arbitrary number of
photon spheres. This seems to contradict the result in [13] which shows that
only one photon sphere can appear in a static spacetime. However, the result
in [13] does not apply in the case when the photon spheres are nested as in
our case.

We remark that our result holds also in the case when the black hole mass
vanishes. However, the family of solutions obtained in this work is different
from the family of solutions obtained in [10]. In the present situation we
require the inner radius of the shell, $R_{0}$, to be sufficiently large,
whereas $R_{0}$ is required to be sufficiently small in [10]. In fact, the
compactness ratio $2M/R\to 8/9$ in both the limits $R_{0}\to 0$ and
$R_{0}\to\infty$.

Let us finally mention that the linear massless Einstein-Vlasov system has
been studied on a fixed black hole spacetime in [4]. The authors show that
solutions to the linear Einstein-Vlasov system on a Kerr background satisfy a
Morawetz estimate. Our result shows that an analogous result cannot hold for
the nonlinear Einstein-Vlasov system. On the other hand, the main purpose of
[4] is to understand perturbations of black hole spacetimes. The steady states
we construct require compact configurations and the matter components cannot
be made arbitrary small. Thus they should not be relevant when studying
perturbations.

The outline of the paper is as follows. In Section 2 we introduce the massless
static Einstein-Vlasov system. In Section 3 we formulate the main results and
in Section 4 we prove our main theorem. Section 5 is devoted to a numerical
investigation of the properties of the solutions.

##  2 The static Einstein -Vlasov system

The metric of a static spherically symmetric spacetime takes the following
form in Schwarzschild coordinates

$ds^{2}=-e^{2\mu(r)}dt^{2}+e^{2\lambda(r)}dr^{2}+r^{2}(d\theta^{2}+\sin^{2}{\theta}d\varphi^{2}),$

where $r\geq 0,\,\theta\in[0,\pi],\,\varphi\in[0,2\pi].$ Asymptotic flatness
is expressed by the boundary conditions

$\lim_{r\rightarrow\infty}\lambda(r)=\lim_{r\rightarrow\infty}\mu(r)=0.$

We now formulate the static massless Einstein-Vlasov system. For an
introduction to the Einstein-Vlasov system we refer to [8], [22] and [24].
Below we use units such that $c=G=1$ where $G$ is the gravitational constant
and $c$ is the speed of light. The static massless Einstein-Vlasov system is
given by the Einstein equations

$\displaystyle\displaystyle e^{-2\lambda}(2r\lambda_{r}-1)+1=8\pi r^{2}\rho,$
(2.1) $\displaystyle\displaystyle e^{-2\lambda}(2r\mu_{r}+1)-1=8\pi r^{2}p,$
(2.2)
$\displaystyle\displaystyle\mu_{rr}+(\mu_{r}-\lambda_{r})(\mu_{r}+\frac{1}{r})=8\pi
p_{T}e^{2\lambda},$ (2.3)

together with the static Vlasov equation

$\frac{w}{\varepsilon}\partial_{r}f-(\mu_{r}\varepsilon-\frac{L}{r^{3}\varepsilon})\partial_{w}f=0,$
(2.4)

where

$\varepsilon=\varepsilon(r,w,L)=\sqrt{w^{2}+L/r^{2}}.$

The matter quantities are defined by

$\displaystyle\rho(r)$ $\displaystyle=$
$\displaystyle\frac{\pi}{r^{2}}\int_{-\infty}^{\infty}\int_{0}^{\infty}\varepsilon(r,w,L)f(r,w,L)\;dLdw,$
(2.5) $\displaystyle p(r)$ $\displaystyle=$
$\displaystyle\frac{\pi}{r^{2}}\int_{-\infty}^{\infty}\int_{0}^{\infty}\frac{w^{2}}{\varepsilon(r,w,L)}f(r,w,L)\;dLdw,$
(2.6) $\displaystyle p_{T}(r)$ $\displaystyle=$
$\displaystyle\frac{\pi}{2r^{4}}\int_{-\infty}^{\infty}\int_{0}^{\infty}\frac{L}{\varepsilon(r,w,L)}f(r,w,L)\;dLdw.$
(2.7)

The variables $w$ and $L$ can be thought of as the momentum in the radial
direction and the square of the angular momentum respectively.

The matter quantities $\rho,p$ and $p_{T}$ are the energy density, the radial
pressure and the tangential pressure respectively. The system of equations
above are not independent and we study the reduced system (2.1)-(2.2) together
with (2.4) and (2.5)-(2.6). It is straightforward to show that a solution to
the reduced system is a solution to the full system.

Define

$E=e^{\mu}\varepsilon,$

then the ansatz

$f(r,w,L)=\Phi(E,L),$ (2.8)

satisfies (2.4). By inserting this ansatz into (2.5)-(2.6) the system of
equations reduce to a system where the metric coefficients $\mu$ and $\lambda$
alone are the unknowns. This has turned out to be an efficient method to
construct static solutions and we will use this approach here. The following
form of $\Phi$ will be used

$\Phi(E,L)=(E_{0}-E)^{k}_{+}(L-L_{0})_{+}^{l},$ (2.9)

where $l\geq 1/2,\,k\geq 0,\,L_{0}>0,\;E_{0}>0,$ and $x_{+}:=\max\\{x,0\\}$.
In the Newtonian case with $l=L_{0}=0,$ this ansatz leads to steady states
with a polytropic equation of state.

The aim in this work is to show that static shells of Vlasov matter exist
which surround a Schwarzschild black hole; in fact there can be arbitrary many
shells separated by vacuum surrounding the black hole. To prove our result we
construct highly compact shells, i.e., shells for which the compactness ratio

$\Gamma:=\sup\frac{2m(r)}{r},$ (2.10)

is large; roughly $\Gamma\geq\frac{4}{5}$. Here $m$ is the Hawking mass
defined for $r\geq 2M_{0}$ by

$m(r)=M_{0}+\int_{2M_{0}}^{r}s^{2}\rho(s)\,ds,\;\;r\geq 2M_{0}.$

From [7] it always holds that $\Gamma<\frac{8}{9}$. The result in [7] concerns
steady states with a regular centre but it is straightforward to show that it
holds also in the case with a Schwarzschild black hole at the centre.

###### Remark 2.1

Numerically we are able to construct solutions where the shell is close to the
photon sphere of the black hole. For such solutions it turns out that $\Gamma$
is larger than, but close to, $2/3$, cf. Section 5. Hence, the presence of a
black hole reduces the required lower bound of $\Gamma$. Indeed, recall that
the numerical studies in the regular case indicates that the required lower
bound is larger in that case, cf. [10], [1] and [17].

If the inner radius of the shell is denoted by $R_{0}$, we show that for
highly relativistic shells there is a radius $R_{1}$ such that $f(r,w,L)=0$ in
an interval $[R_{1},R_{1}+\epsilon],\epsilon>0$. This fact makes it possible
to glue a Schwarzschild solution at $r=R_{1}$ to the shell solution with
support in $[R_{0},R_{1}]$. If a Schwarzschild solution is not attached at
$r=R_{1}$, then the ansatz (2.9) implies that Vlasov matter will occur again
and there exists a radius $R_{2}$ such that $f>0$ for all $r>R_{2}$ and the
solution is not asymptotically flat. This is a general feature of massless
static solutions of the Einstein-Vlasov system obtained from an ansatz, cf.
equation (2.8). In the massive case the situation is different and solutions
generated by the ansatz (2.9) alone gives rise to compactly supported
solutions.

In the massive case the existence of shells surrounding a Schwarzschild black
hole was settled in [22]. These shells are not highly relativistic. To
construct highly relativistic shells for which $\Gamma$ is sufficiently large
we adapt the method developed in [6], which in turn was used to show existence
of massless steady states with a regular centre in [10].

##  3 Set up and main result

Let $M_{0}>0$ be the mass of the black hole. In a vacuum region in the
exterior of the black hole it holds that

$e^{2\mu(r)}=1-\frac{2M_{0}}{r}.$

Note that the ansatz (2.9) implies that $f=0$ whenever $E>E_{0}$. Accordingly
we let $f=0$ in the interval $[2M_{0},R_{0}],$ where $R_{0}$ is the largest
root to the equation

$\big{(}1-\frac{2M_{0}}{r})\frac{L_{0}}{r^{2}}=E_{0}^{2}=1.$ (3.11)

Of course, we need a condition on $L_{0}$ which guarantees that the equation
has real roots.

###### Remark 3.1

The parameter $L_{0}$ can be removed. From above we see that by replacing
$E_{0}$ by $\tilde{E_{0}}=E_{0}/\sqrt{L_{0}}$ we could consider the case
$L_{0}=1$ and use $E_{0}$ as free parameter, cf. [25, 17]. However, below we
keep $L_{0}$ as free parameter and we fix $E_{0}=1$.

An elementary computation shows that the maximum value of the left hand side
of (3.11) is

$\frac{L_{0}}{27M_{0}^{2}},$

attained at $r=3M_{0}$. This radius corresponds to the radius of the photon
sphere of the black hole. We fix $E_{0}=1$ and impose the condition that

$L_{0}>27M_{0}^{2}=:L_{*}.$

Equation (3.11) then has three real roots and we denote the largest root by
$R_{0}$. To carry out the proof of our main result we will take $L_{0}$ large.
We have the following result.

###### Lemma 3.2

There exists a constant $C>0$, depending on $M_{0}$, such that

$1-\frac{C}{\sqrt{L_{0}}}\leq\frac{R_{0}}{\sqrt{L_{0}}}\leq 1\mbox{ as
}L_{0}\to\infty.$

Proof: The proof is a straightforward application of the solution formula for
cubic equations. Indeed, equation (3.11) is equivalent to

$r^{3}-L_{0}r+2M_{0}L_{0}=0,$

which we write as $r^{3}+pr+q=0$ where $p=-L_{0}$ and $q=2M_{0}L_{0}$. Set

$D=(\frac{p}{3})^{3}+(\frac{q}{2})^{2}=L_{0}^{2}(-\frac{L_{0}}{27}+M_{0}^{2}).$

Since we choose $L_{0}>27M_{0}^{2}$ it follows that $D<0$ and a standard
result for cubic equations implies that the equation has three real roots
given by

$r_{1}=u+v,\;\;r_{2,3}=\frac{u+v}{2}\pm\frac{u-v}{2}i\sqrt{3},$

where $u=(-q/2+i\sqrt{|D|})^{1/3}$ and $v=(-q/2-i\sqrt{|D|})^{1/3}$. Note here
that $u+v$ and $(u-v)i$ are real. The largest of these roots is

$R_{0}:=r_{1}=u+v=2(\frac{q^{2}}{4}+|D|)^{1/6}\cos{\frac{\alpha}{3}},$

where

$\alpha=\frac{\pi}{2}+\arctan{\frac{(q/2)}{\sqrt{|D|}}}.$

We have

$0\leq\frac{(q/2)}{\sqrt{|D|}}=\leq\frac{M_{0}}{\sqrt{L_{0}}\sqrt{\frac{1}{27}-\frac{M_{0}^{2}}{L_{0}}}}\leq\frac{C}{\sqrt{L_{0}}},$

where we used that $L_{0}-L_{*}>0$ since we consider large $L_{0}$. Hence
there exists a constant $C>0$ such that

$\frac{\pi}{2}\leq\alpha\leq\frac{\pi}{2}+\frac{C}{\sqrt{L_{0}}},$

for large $L_{0}$. This implies that

$\frac{\sqrt{3}}{2}-\frac{C}{\sqrt{L_{0}}}\leq\cos{\frac{\alpha}{3}}\leq\frac{\sqrt{3}}{2},$

for some positive constant $C$. Moreover, we have

$2(\frac{q^{2}}{4}+|D|)^{1/6}=2(M_{0}^{2}L_{0}^{2}+L_{0}^{2}(\frac{L_{0}}{27}-M_{0}^{2}))^{1/6}=\frac{2\sqrt{L_{0}}}{\sqrt{3}}.$

Hence we conclude that for large $L_{0}$

$1-\frac{C}{\sqrt{L_{0}}}\leq\frac{R_{0}}{\sqrt{L_{0}}}\leq 1.$

$\Box$

We now consider the Einstein-Vlasov system on the domain $r\geq R_{0}$ and we
prescribe data for $\mu$ and $\lambda$ at $r=R_{0}$ by letting

$e^{2\mu(R_{0})}=e^{-2\lambda(R_{0})}=1-\frac{2M_{0}}{R_{0}}.$

Using results from previous works we can assume that there exists a solution
to the Einstein-Vlasov system which exists on $[R_{0},\infty[$ with the
property that $\Gamma<8/9$ everywhere, cf. [23] for existence and [7] for the
bound on $\Gamma$. We will show that if $R_{0}$ is sufficiently large, i.e. we
take $L_{0}$ large, there is a radius $R_{1}>R_{0}$ such that the energy
density and the pressure components vanish at $r=R_{1}$. The shell
$[R_{0},R_{1}]$ is thin in the sense that

$\frac{R_{1}}{R_{0}}\to 1\mbox{as }R_{0}\to\infty.$ (3.12)

However, the difference $R_{1}-R_{0}$ does not need to decrease. Depending on
the parameters in equation (2.9) the difference may even become unbounded, cf.
Remark (3.4). Hence the shell is thin in the sense (3.12), which is different
from the usual notion of thin shells in general relativity.

When the distribution function $f$ has the form

$f(r,w,L)=\Phi(E,L),$ (3.13)

where $\Phi=0$ whenever $E>E_{0}$, the matter quantities $\rho$ and $p$ become
functionals of $\mu$, and we have

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\frac{2\pi}{r^{2}}\int_{\frac{\sqrt{L_{0}}}{r}}^{E_{0}e^{-\mu}}\int_{L_{0}}^{r^{2}(s^{2}-1)}\Phi(e^{\mu}s,L)\frac{s^{2}}{\sqrt{s^{2}-1-\frac{L}{r^{2}}}}dLds,$
(3.14) $\displaystyle p$ $\displaystyle=$
$\displaystyle\frac{2\pi}{r^{2}}\int_{\frac{\sqrt{L_{0}}}{r}}^{E_{0}e^{-\mu}}\int_{L_{0}}^{r^{2}(s^{2}-1)}\Phi(e^{\mu}s,L)\sqrt{s^{2}-1-\frac{L}{r^{2}}}dLds.$
(3.15)

Here we have kept the parameter $E_{0}$ but in what follows we use that
$E_{0}=1$. By taking (2.9) for $\Phi$ these integrals can be computed
explicitly in the cases when $k=0,1,2,...$ and $l=1/2,3/2,...$ Let

$\gamma=-\mu-\frac{1}{2}\log{\frac{L_{0}}{r^{2}}},$

we then have the following lemma from [10].

###### Lemma 3.3

Let $k=0,1,2,...$ and let $l=1/2,3/2,5/2,...$ then there are positive
constants $\pi^{j}_{k,l},\,j=1,2,3$ such that when $\gamma\geq 0$

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\pi^{1}_{k,l}r^{2l}(\frac{L_{0}}{r^{2}})^{l+2}(e^{\gamma}-1)^{l+k+3/2}P_{l+5/2-k}(e^{\gamma}),$
(3.17) $\displaystyle p$ $\displaystyle=$
$\displaystyle\pi^{2}_{k,l}r^{2l}(\frac{L_{0}}{r^{2}})^{l+2}(e^{\gamma}-1)^{l+k+5/2}P_{l+3/2-k}(e^{\gamma}),$
(3.18)

If $\gamma<0$ then all matter components vanish. Here $P_{n}(e^{\gamma})$ is a
polynomial of degree $n$ and $P_{n}>0.$

In order to simplify the technical details we consider only the case $k=0$ and
$l=1/2$ but we emphasize that our result holds more generally as in [10] and
[6].

###### Remark 3.4

In the case $k=0$ and $l=1/2$ the support of the shell $[R_{0},R_{1}]$ satisfy
$R_{1}-R_{0}\leq C$, independently of $R_{0}$ as shown below, whereas for
other values of the parameters we claim that

$R_{1}-R_{0}\sim R_{0}^{\frac{q-2-2l}{q}},$

where $q=k+l+5/2$, cf. Section 5. This relation is obtained by performing the
analysis below in the general case, cf. [6].

Let $t=e^{\gamma}$, we then have (with $k=0$ and $l=1/2$)

$\rho(r)=\pi^{2}r\big{(}\frac{L_{0}}{r^{2}}\big{)}^{5/2}(t-1)^{2}\big{(}\frac{3t^{3}+6t^{2}+4t+2}{15}\big{)},$
(3.19)

and

$p(r)=\pi^{2}r\big{(}\frac{L_{0}}{r^{2}}\big{)}^{5/2}(t-1)^{3}\big{(}\frac{3t^{2}+9t+8}{60}\big{)}.$
(3.20)

Let $R_{0}$ be large and define

$\delta:=\big{(}\frac{1}{20\pi^{3}}\big{)}^{1/3}.$ (3.21)

The argument below will be carried out in the interval $I:=[R_{0},R_{1}]$
where $R_{1}<R_{0}+50\delta=:R_{1}^{*}$. Since $R_{0}$ will be taken large,
$R_{1}/R_{0}$ is as close to $1$ as we wish. We now formulate the main results
in this work.

###### Theorem 3.5

Let $M_{0}\geq 0$ be the ADM mass of a Schwarzschild black hole. Then there
exist static solutions with finite ADM mass to the massless spherically
symmetric Einstein-Vlasov system surrounding the black hole. The matter
components are supported on a finite interval $[R_{0},R_{1}]$, where
$R_{0}>3M_{0}$, and spacetime is asymptotically flat.

###### Remark 3.6

The arguments below lead to a solution for which $\mu(r)$ has a finite limit
$\mu(\infty)$ as $r\to\infty$. In order to obtain an asymptotically flat
solution we rescale by letting $\tilde{E_{0}}:=e^{\mu(\infty)}$ and
$\tilde{\mu}(r):=\mu(r)-\mu(\infty)$.

###### Remark 3.7

The regularity of the solution depends on the parameters $k$ and $l$, cf.
equations (3.17) and (3.18). For the values of $k$ and $l$ that we consider in
this work the matter quantities are continuously differentiable.

###### Remark 3.8

As mentioned in the introduction, note that our result holds also in the case
when $M_{0}=0$ but that the family of solutions obtained in this work is
different from the family of solutions obtained in [10]. In the present
situation we require the inner radius $R_{0}$ to be large whereas in [10]
$R_{0}$ is required to be small. Clearly, when $M_{0}>0$ it is not possible to
take $R_{0}$ small since necessarily $R_{0}>2M_{0}$. Both families share the
property that $\Gamma\to 8/9$ in the extreme limits, i.e., when $R_{0}\to 0$
as in [10] or when $R_{0}\to\infty$ as in the present case.

###### Remark 3.9

Having a black hole with one shell surrounding it, we can start from this
solution and add another shell with the strategy in the proof. Hence, our
result implies that an arbitrary number of shells can surround the black hole.

In view of the discussion in Section 2, Theorem 3.5 is a consequence of the
following result.

###### Theorem 3.10

Consider a static solution to the massless Einstein-Vlasov system,
corresponding to the ansatz (2.9) with $k=0$ and $l=1/2$, with data given at
$r=R_{0}$. For $R_{0}$ sufficiently large, there exists $R_{1}\leq
R_{0}+50\delta$ such that the matter components are supported in the interval
$[R_{0},R_{1}]$ and vanish at $r=R_{1}$. Here $\delta$ is given by (3.21).
Furthermore, $\Gamma\to\frac{8}{9}$ as $R_{0}\to\infty$.

##  4 Proof of main result

The proof of our main result will follow from a chain of lemmas.

Proof of Theorem 3.10. First we establish convenient formulas for the matter
terms. Although $L_{0}$ is our free parameter, we will instead use $R_{0}$ as
free parameter since $R_{0}\to\infty$ as $L_{0}\to\infty$ in view of Lemma
3.2. Next we note that $R_{1}^{*}/R_{0}\leq 1+C/R_{0}$. Hence, by Lemma 3.2 we
have for any $r\in[R_{0},R_{1}^{*}]$

$1-\alpha(R_{0})\leq{\frac{L_{0}}{r^{2}}}\leq 1+\alpha(R_{0}),$ (4.22)

where $\alpha(R_{0})\geq 0$ and $\alpha(R_{0})\to 0$ as $R_{0}\to\infty$. By
the definition of $R_{0}$ we have that $\gamma(R_{0})=0$. Moreover,

$\gamma^{\prime}(r)=-\mu^{\prime}(r)+\frac{1}{r}.$ (4.23)

Now, since $\rho=p=0$ and $m(r)=M_{0}$ for $r\leq R_{0}$, we have that that
$\mu^{\prime}(R_{0})<1/R_{0}$. Here we used that $M_{0}/R_{0}<1/3$ and that

$e^{2\lambda(r)}=\frac{1}{1-\frac{2m(r)}{r}}.$

The last relation is a consequence of the Einstein equation (2.1). This
implies that $\gamma^{\prime}(R_{0})>0$. Thus $\gamma(r)>0$ in a right
neighborhood of $R_{0}$ and the aim is to show that there exists
$R_{1}<R_{1}^{*}$ such that $\gamma(R_{1})=0$. Hence $\gamma>0$ on the
interior of the closed interval $I$. An upper bound on $\gamma$ follows from
(4.23). We have since $\mu^{\prime}(r)\geq 0$,

$\gamma(r)=\gamma(r)-\gamma(R_{0})\leq\log\frac{r}{R_{0}}=\log(1+\frac{r-R_{0}}{R_{0}}),$

which implies that $0\leq\gamma\leq\alpha(R_{0})$ on $I$ if $R_{0}$ is
sufficiently large. Hence for $r\in I$

$\gamma(r)\leq e^{\gamma(r)}-1\leq\gamma(r)(1+\alpha(R_{0})).$

In particular this implies that (recall $t=e^{\gamma}$)

$1\leq\frac{3t^{3}+6t^{2}+4t+2}{15}\leq 1+\alpha(R_{0}),$

and similarly for the corresponding polynomial in $p(r)$. Putting these
estimates together we conclude that

$\pi^{2}r\gamma^{2}(r)(1-\alpha(R_{0}))\leq\rho(r)\leq\pi^{2}r\gamma^{2}(r)(1+\alpha(R_{0})),$

and similarly for $p(r)$. Since for all arguments below it is sufficient to
have a lower and an upper bound on $\rho$ and $p$ we assume for simplicity
that

$\rho(r)=\pi^{2}r\gamma^{2}(r),$ (4.24)

and

$p(r)=\frac{\pi^{2}}{3}r\gamma^{3}(r),$ (4.25)

for $r\in I$.

###### Lemma 4.1

Let $\delta$ be as above. Then

$\gamma^{\prime}(r)\geq\frac{1}{2r}\;\mbox{for
}r\in[R_{0},R_{0}+\delta]=:I_{1}.$

Proof: We have by the mean value theorem that for any $\sigma\leq\delta$

$\gamma(R_{0}+\sigma)=\gamma(R_{0}+\sigma)-\gamma(R_{0})=\sigma\gamma^{\prime}(\xi)\leq\frac{\delta}{R_{0}},$

where $\xi\in[R_{0},R_{0}+\sigma]$, since $\gamma(r)\leq 1/r$. Hence,

$\rho(r)\leq\pi^{2}r\frac{\delta^{2}}{R_{0}^{2}}\;\mbox{ on }I_{1}.$

We get for $\sigma\leq\delta$

$m(R_{0}+\sigma)\leq
M_{0}+\frac{4\pi^{3}\delta^{2}}{R_{0}^{2}}\int_{R_{0}}^{R_{0}+\sigma}\eta^{3}\,d\eta\leq
M_{0}+\frac{4\pi^{3}\delta^{2}(R_{0}+\sigma)^{3}}{R_{0}^{2}}\sigma.$

By taking $R_{0}$ large we thus obtain

$\frac{m(R_{0}+\sigma)}{R_{0}+\sigma}\leq
4\pi^{3}\delta^{2}\sigma+\alpha(R_{0})\leq
4\pi^{3}\delta^{3}+\alpha(R_{0})\leq\frac{1}{5}+\alpha(R_{0}).$

This implies that for $r\in I_{1}$

$\frac{m(r)}{r}e^{2\lambda(r)}\leq\frac{1}{3}+\alpha(R_{0}).$

Moreover, we have for $r\in I_{1}$ that

$4\pi r^{2}p(r)e^{2\lambda(r)}\leq 4\pi
r^{2}p(r)(\frac{5}{3}+\alpha(R_{0}))\leq\frac{4\pi^{3}\delta^{3}}{3}(\frac{5}{3}+\alpha(R_{0}))\leq\frac{1}{9}+\alpha(R_{0}).$

Hence, for $r\in I_{1}$ we get for sufficiently large $R_{0}$

$\mu^{\prime}(r)\leq(\frac{1}{3}+\frac{1}{9}+\alpha(R_{0}))\frac{1}{r}\leq\frac{1}{2r}.$

This completes the proof of the lemma. $\Box$

The lemma implies that for $\sigma\leq\delta$

$\gamma(R_{0}+\sigma)\geq\gamma(R_{0})+\sigma\inf_{0\leq
s\leq\sigma}\gamma^{\prime}(R_{0}+s)\geq\frac{\sigma}{2(R_{0}+\sigma)}.$

Let

$\sigma_{*}:=\delta,$

and define

$\gamma_{*}=\frac{\sigma_{*}}{2(R_{0}+\sigma_{*})}.$

###### Remark 4.2

The notation $\sigma_{*}$ is in this case superfluous but for general
parameter values it is motivated, cf. [6].

Clearly we have that

$\gamma(R_{0}+\sigma_{*})\geq\gamma_{*}.$

The result in the following lemma shows that $\gamma$ will reach the
$\gamma_{*}$ level again at a larger radius.

###### Lemma 4.3

There exists a radius $r_{2}>R_{0}+\sigma_{*}$ such that
$\gamma(r_{2})=\gamma_{*}$ and such that $\gamma(r)>\gamma_{*}$ for
$R_{0}+\sigma_{*}<r<r_{2}$. Moreover, $r_{2}\leq R_{0}+11\delta$.

Proof: Since $\gamma(R_{0}+\sigma_{*})\geq\gamma_{*},$ and since Lemma 4.1
gives that $\gamma^{\prime}(R_{0}+\sigma_{*})>0,$ the radius $r_{2}$ must be
strictly greater than $R_{0}+\sigma_{*}.$ Let
$[R_{0}+\sigma_{*},R_{0}+\sigma_{*}+\Delta],$ for some $0<\Delta<10\delta,$ be
the smallest interval such that $\gamma\geq\gamma_{*}$ on the interval. We
will show that there in fact $\gamma(R_{0}+\sigma_{*}+\Delta)=\gamma_{*}$. We
have from (4.24)

$\displaystyle m(R_{0}+\sigma_{*}+\Delta)$ $\displaystyle\geq$ $\displaystyle
4\pi^{3}\int_{R_{0}+\sigma_{*}}^{R_{0}+\sigma_{*}+\Delta}\frac{\sigma_{*}^{2}}{4(R_{0}+\sigma_{*})^{2}}r^{3}\,dr$
(4.26) $\displaystyle\geq$
$\displaystyle\pi^{3}\sigma_{*}^{2}(R_{0}+\sigma_{*})\Delta.$

Hence,

$\frac{m(R_{0}+\sigma_{*}+\Delta)}{R_{0}+\sigma_{*}+\Delta}\geq\frac{\pi^{3}\sigma_{*}^{2}(R_{0}+\sigma_{*})}{R_{0}+\sigma_{*}+\Delta}\Delta\geq\frac{9}{10}\pi^{3}\sigma_{*}^{2}\Delta,$

where we used that $R_{0}\geq 100\delta$, so that
$(R_{0}+\sigma_{*})/(R_{0}+\sigma_{*}+\Delta)\geq 9/10$, since we consider
large $R_{0}$. Since for $r\geq R_{0}$ necessarily

$\frac{m(r)}{r}<\frac{4}{9},$

we conclude that there is a $\Delta$ such that

$\Delta\leq\frac{40}{81\pi^{3}\sigma_{*}^{2}}=\frac{40}{81\pi^{3}\delta^{2}}=\frac{40\cdot
20^{2/3}}{81\pi}\leq 10\delta,$

with the property that $\gamma(R_{0}+\sigma_{*}+\Delta)=\gamma_{*}$, since
otherwise we obtain a contradiction. $\Box$

Next we show an important property of the solution at the radius $r=r_{2}$
when $R_{0}$ is sufficiently large.

###### Lemma 4.4

Let $r_{2}$ be as above. If $R_{0}$ is sufficiently large then

$\frac{m(r_{2})}{r_{2}}\geq\frac{2}{5}.$

Proof:  We consider the fundamental equation ($10$) in [5]. In our case with a
black hole at the center it takes the form

$(\frac{m(r)}{r^{2}}+4\pi
rp(r))e^{\mu(r)+\lambda(r)}-\frac{M_{0}}{R_{0}^{2}}=\frac{1}{r^{2}}\int_{R_{0}}^{r}4\pi\eta^{2}e^{\mu+\lambda}(\rho+p+2p_{T})d\eta.$
(4.27)

Here we used that $\mu(R_{0})+\lambda(R_{0})=0$ and that $p(R_{0})=0$. This
equation is a consequence of the generalized Oppenheimer-Tolman-Volkov
equation. In the present massless case we have $p+2p_{T}=\rho$, so that
$\rho+p+2p_{T}=2\rho$. If we take $r=r_{2}$ we then get

$\displaystyle\frac{m(r_{2})}{r_{2}}e^{(\mu+\lambda)(r_{2})}$
$\displaystyle\geq$
$\displaystyle\frac{8\pi}{r_{2}}\int_{R_{0}}^{r_{2}}\eta^{2}e^{\mu+\lambda}\rho\,d\eta-4\pi
r_{2}^{2}pe^{(\mu+\lambda)(r_{2})}.$ (4.28)

Here we dropped the term involving $M_{0}$ due to sign. Note that it becomes
arbitrary small for $R_{0}$ large. Next we write

$\displaystyle\int_{R_{0}}^{r_{2}}4\pi\eta\rho e^{\lambda}\,d\eta$
$\displaystyle=$
$\displaystyle\int_{R_{0}}^{r_{2}}(-\frac{d}{dr}e^{-\lambda})d\eta+\int_{R_{0}}^{r_{2}}\frac{me^{\lambda}}{\eta^{2}}d\eta$
(4.29) $\displaystyle\geq$
$\displaystyle\sqrt{1-\frac{2M_{0}}{R_{0}}}-\sqrt{1-\frac{2m(r_{2})}{r_{2}}}$
$\displaystyle=$ $\displaystyle
1-\sqrt{1-\frac{2m(r_{2})}{r_{2}}}-\big{(}1-\sqrt{1-\frac{2M_{0}}{R_{0}}}\big{)}$
$\displaystyle=$
$\displaystyle\frac{2m(r_{2})}{r_{2}(1+\sqrt{1-\frac{2m(r_{2})}{r_{2}}})}-\alpha(R_{0}).$

In the inequality above we dropped the second term due to sign. We note that
this term is as small as we wish for a relatively thin shell. The reason we
point this out is that the chain of inequalities in this paragraph is close to
a chain of equalities for very thin shells. This is essential to understand
why for a thin shell $\Gamma$ approaches the limit $8/9$.

From (4.28) we obtain using that $\mu$ is increasing

$\displaystyle Y\frac{m(r_{2})}{r_{2}}e^{(\mu+\lambda)(r_{2})}$
$\displaystyle\geq$
$\displaystyle\frac{8\pi}{r_{2}}\int_{R_{0}}^{r_{2}}\eta^{2}e^{\mu+\lambda}\rho\,d\eta-4\pi
r_{2}^{2}p(r_{2})e^{(\mu+\lambda)(r_{2})}$ $\displaystyle\geq$ $\displaystyle
8\pi e^{\mu(R_{0})}\frac{R_{0}}{r_{2}}\int_{R_{0}}^{r_{2}}\eta
e^{\lambda}\rho\,d\eta-4\pi r_{2}^{2}p(r_{2})e^{(\mu+\lambda)(r_{2})}.\ $

Let us introduce the notation $P:=4\pi r_{2}^{2}p(r_{2})$. Using (4.29) we get

$\displaystyle\frac{m(r_{2})}{r_{2}}e^{(\mu+\lambda)(r_{2})}$
$\displaystyle\geq$
$\displaystyle\frac{e^{\mu(R_{0})}R_{0}}{r_{2}}\frac{4m(r_{2})}{r_{2}(1+\sqrt{1-\frac{2m(r_{2})}{r_{2}}})}$
(4.30) $\displaystyle-Pe^{(\mu+\lambda)(r_{2})}-\alpha(R_{0}).$

Letting

$Y:=\frac{m(r_{2})}{r_{2}},$

we obtain by using that $e^{-\lambda(r_{2})}=\sqrt{1-2Y(r_{2})}$

$Y\geq
e^{\mu(R_{0})-\mu(r_{2})}\left(\frac{R_{0}}{r_{2}}\right)\frac{4Y\sqrt{1-2Y}}{1+\sqrt{1-2Y}}-P-\alpha(R_{0})e^{-(\mu+\lambda)(r_{2})}.$
(4.31)

Next we show a couple of properties of the solutions that we need to proceed
with the argument. From Lemma 4 we have that $\gamma$ approaches $\gamma_{*}$
from above and therefore $\gamma^{\prime}(r_{2})\leq 0,$ which implies that

$\mu^{\prime}(r_{2})\geq\frac{1}{r_{2}},$ (4.32)

in view of (4.23). Furthermore we have from (4.25) that

$p(r_{2})=\frac{\pi^{2}}{3}r_{2}\gamma_{*}^{3}=\frac{\pi^{2}}{3}r_{2}\frac{\delta^{3}}{8(R_{0}+\delta)^{3}}.$

Now recall that

$\mu^{\prime}=(\frac{m}{r^{2}}+4\pi rp)e^{2\lambda}.$

We will show that $Y=m(r_{2})/r_{2}\geq 1/4$. Assume the contrary, i.e.,
assume

$Y<\frac{1}{4},$ (4.33)

then

$e^{2\lambda(r_{2})}=\frac{1}{1-2Y}<2,$

and thus

$4\pi
r_{2}^{2}p(r_{2})e^{2\lambda(r_{2})}\leq\frac{\pi^{3}}{3}r_{2}^{3}\frac{\delta^{3}}{(R_{0}+\delta)^{3}}=\frac{1}{60}\frac{r_{2}^{3}}{(R_{0}+\delta)^{3}}\leq\frac{1}{50}.$
(4.34)

Here we used that $\delta^{3}=1/20\pi^{3}$ and that $r_{2}\leq R_{0}+11\delta$
so that

$\frac{r_{2}^{3}}{(R_{0}+\delta)^{3}}\leq\frac{6}{5}\mbox{ for }R_{0}\mbox{
large}.$

Hence if (4.33) holds then

$\mu^{\prime}(r_{2})=(\frac{m(r_{2})}{r_{2}^{2}}+4\pi
r_{2}p(r_{2}))e^{2\lambda(r_{2})}\leq(\frac{1}{2}+\frac{1}{50})\frac{1}{r_{2}},$

which is a contradiction to (4.32) and we obtain that

$Y\geq\frac{1}{4}.$

Let us next consider the difference $\mu(R_{0})-\mu(r_{2}).$ We have

$\mu(r_{2})-\mu(R_{0})=\int_{R_{0}}^{r_{2}}(\frac{m}{\eta^{2}}+4\pi\eta
p)e^{2\lambda}\,d\eta.$

Since $\gamma^{\prime}(r)\leq 1/r$ we have for $R_{0}\leq r\leq r_{2}\leq
R_{0}+11\delta$

$\gamma(r)\leq\log{\frac{r}{R_{0}}}\leq\frac{11\delta}{R_{0}},$

and we get by using (4.25), and that $e^{2\lambda}\leq 9$

$4\pi r^{2}p(r)e^{2\lambda(r)}\leq
12\pi^{3}r^{3}\gamma^{3}(r)\leq\frac{3}{5}\big{(}\frac{r}{R_{0}}\big{)}^{3}11^{3}\leq
C.$

Hence, since $m/r\leq 4/9$ we get

$\int_{R_{0}}^{r_{2}}(\frac{m}{\eta^{2}}+4\pi\eta
p)e^{2\lambda}\,d\eta\leq\int_{R_{0}}^{r_{2}}\frac{C}{\eta}\,d\eta\leq
C\log{\frac{r_{2}}{R_{0}}}\leq\frac{C}{R_{0}}.$

As a consequence we obtain that

$1-\alpha(R_{0})\leq
e^{\mu(R_{0})-\mu(r_{2})}\left(\frac{R_{0}}{r_{2}}\right)\leq
1+\alpha(R_{0}).$

From this estimate we also have that $e^{-(\mu+\lambda)(r_{2})}$ is bounded
independently of $R_{0}$ and using that $Y\geq 1/4$ we can therefore write
(4.31) as

$1\geq\frac{4\sqrt{1-2Y}}{1+\sqrt{1-2Y}}-4P-\alpha(R_{0}).$ (4.35)

Let us introduce the notation $B=4P+\alpha(R_{0})$. We can assume that $B\geq
0$ since otherwise this term can be dropped. Multiplying (4.35) by

$1+\sqrt{1-2Y},$

and then squaring both sides we obtain

$Y\geq\frac{4}{9}-\frac{B\sqrt{1-2Y}(1+\sqrt{1-2Y})}{3},$ (4.36)

where we dropped the term involving $B^{2}$ due to sign. Using that $Y\geq
1/4$ we estimate

$\frac{\sqrt{1-2Y}(1+\sqrt{1-2Y})}{3}\leq\frac{1}{3\sqrt{2}}(1+\frac{1}{\sqrt{2}})\leq\frac{9}{20}.$

Taking $R_{0}$ large we derive from (4.36) the inequality

$Y\geq\frac{4}{9}-2P.$

The estimate (4.34) can now be used noting that $e^{2\lambda(r_{2})}\geq 2$
and we obtain

$Y\geq\frac{4}{9}-\frac{1}{50}\geq\frac{2}{5}.$

$\Box$

###### Remark 4.5

The proof above is very similar to the corresponding proof in [6]. However,
for shells for which the inner radius $R_{0}\to 0$, as in [6], the pressure
term $P\to 0$ and one can conclude from the argument in the lemma above that
as $R_{0}\to 0$ the compactness ratio $\Gamma\to 8/9$. In the present case the
situation is slightly different. However, as soon as we know that there is a
radius $R_{1}$ such that $\rho(R_{1})=p(R_{1})=0$, and such that
$R_{1}/R_{0}\to 1$ as $R_{0}\to\infty$, then we can use the argument above,
with $r_{2}$ replaced by $R_{1}$, to show that the compactness ratio
$\Gamma\to\frac{8}{9}$ as $R_{0}\to\infty$. This shows the last claim in
Theorem 3.10.

Inspired by an idea of T. Makino introduced in [19], we show that $\gamma$
necessarily must vanish close to the point $r_{2}$ if $R_{0}$ is sufficiently
large. Let

$x:=\frac{m(r)}{r\gamma(r)}.$

Using that $m^{\prime}(r)=4\pi r^{2}\rho,$ and that
$\mu^{\prime}(r)=(m/r^{2}+4\pi r)e^{2\lambda}$ it follows that

$rx^{\prime}=\frac{4\pi r^{2}\rho}{\gamma}-x+\frac{x^{2}}{1-2\gamma
x}-\frac{x}{\gamma}+\frac{4\pi xr^{2}p(r)}{\gamma(1-2\gamma x)}.$

In our case $r>R_{0}$ and $\gamma>0$ and we will show that $\gamma(r)=0$ for
some $R_{1}\in I$. Since $\gamma>0$ and $\rho,p\geq 0$ the first term and the
last term can be dropped and we have

$rx^{\prime}\geq-x+\frac{x^{2}}{1-2\gamma
x}-\frac{x}{\gamma}=\frac{x^{2}}{3(1-2\gamma x)}-x+\frac{2x^{2}}{3(1-2\gamma
x)}-\frac{x}{\gamma}.$ (4.37)

Take $R_{0}$ sufficiently large so that $m(r_{2})/r_{2}\geq 2/5$ by Lemma 4.4.
Let $r\in[r_{2},16r_{2}/15],$ then since $m$ is increasing in $r$ we get

$\frac{m(r)}{r}\geq\frac{m(r_{2})}{r}=\frac{r_{2}}{r}\frac{m(r_{2})}{r_{2}}\geq\frac{15}{16}\cdot\frac{2}{5}=\frac{3}{8}.$

Now by the definition of $x$ it follows that

$\frac{x}{1-2\gamma x}=\frac{m}{\gamma r}e^{2\lambda}=\frac{m}{\gamma
r(1-2m/r)}\geq\frac{3}{2\gamma},\mbox{ when }\frac{m}{r}\geq\frac{3}{8}.$

Thus on $[r_{2},16r_{2}/15],$

$\frac{2x^{2}}{3(1-2\gamma x)}-\frac{x}{\gamma}\geq 0,$

so that on this interval

$rx^{\prime}\geq\frac{x^{2}}{3(1-2\gamma x)}-x\geq\frac{4}{3}x^{2}-x,$ (4.38)

where we used that

$\frac{1}{1-2\gamma x}=\frac{1}{1-2m/r}\geq 4\mbox{ when
}\frac{m}{r}\geq\frac{3}{8}.$

The upper bound

$\gamma\leq\log(1+\frac{r-R_{0}}{R_{0}})$

implies that

$x(r_{2})=\frac{m(r_{2})}{r_{2}\gamma(r_{2})}\geq\frac{2}{5}\frac{R_{0}}{r_{2}-R_{0}}.$
(4.39)

In particular

$\frac{x(r_{2})}{x(r_{2})-3/4}\leq\frac{16}{15},$

for $R_{0}$ large. Solving (4.38) yields

$x(r)\geq\frac{3}{4}\left(1-\frac{r(4x(r_{2})/3-1)}{4r_{2}x(r_{2})/3}\right)^{-1},\mbox{
on }r\in[r_{2},16r_{2}/15),$

and we get that $x(r)\to\infty$ as $r\to R_{1},$ where

$R_{1}\leq r_{2}\frac{x(r_{2})}{x(r_{2})-3/4}.$ (4.40)

In view of (4.39) we estimate

$R_{1}\leq
r_{2}+\frac{\frac{15}{8}(r_{2}-R_{0})r_{2}}{R_{0}-\frac{15}{8}(r_{2}-R_{0})}.$

For sufficiently large $R_{0}$ we have
$R_{0}/10\geq\frac{15}{8}(r_{2}-R_{0})\geq 0$, recall that $r_{2}-R_{0}\leq
11\delta$, and we obtain

$R_{1}\leq r_{2}+\frac{\frac{25}{12}(r_{2}-R_{0})r_{2}}{R_{0}}\leq
r_{2}+3(r_{2}-R_{0})\leq R_{0}+50\delta.$

This completes the proof of the theorem. $\Box$

##  5 Numerical results

In this section we present some numerical results. In the analytic proof we
required $L_{0}$, or equivalently $R_{0}$, to be large. One aim is to
investigate for which $L_{0}$ solutions can be obtained numerically. Recall
that necessarily $L_{0}>L_{*}$, and one question is if solutions can be
constructed for $L_{0}$ values arbitrary close to $L_{*}$, i.e., if the shell
can occur arbitrary close to the photon sphere of the black hole. It turns out
that this depends on the parameters and on the mass of the black hole, but for
a large class of solutions it is possible. Hence, the condition in the proof
that $L_{0}$ is large is mainly a technical condition. Moreover, based on our
numerical findings together with the arguments in the proof, we conjecture
that the compactness ratio for any shell surrounding a black hole satisfies
$\Gamma>2/3$. Note that $2M_{0}/r=2/3$ at the radius of the photon sphere.
This study is presented in Subsection 5.1.

Another aim is to numerically confirm the claim in Remark (3.4) that the
support $[R_{1},R_{0}]$ of the solution satisfies

$R_{1}-R_{0}\sim R_{0}^{\frac{q-2-2l}{q}},$

where $q=k+l+5/2$. Hence, the thickness of the shells depend on the parameter
values $k$ and $l$. This investigation is presented in Subsection 5.2.

In Subsection 5.3 we construct solutions with several shells. We do this by
using two different strategies. From the set up in Section 2 we know that
having a compactly supported solution with ADM mass $M_{0}$, where $M_{0}$ in
this case is the total mass of the black hole and the shell(s) surrounding the
black hole, we can choose $L_{0}>L_{*}$ (sufficiently large) and find the
largest root $R_{0}$ to equation (3.11). We can then pose data at $r=R_{0}$
and numerically construct an additional shell. This strategy thus follows the
analytic proof in the previous section.

The second strategy is different in the sense that it does not follow the
strategy of the proof. We start with a black hole and then we use one ansatz
function for the shells, in particular we do not change $L_{0}$ for the
different shells. In this way a solution consisting of several shells,
separated by vacuum, is generated by solving equation (3.11) only once. After
a number of shells there is no longer a vacuum region separating the
neighboring shells and a Schwarzschild solution has to be glued before this
happens. The solution obtained in this way is very similar to the multi-peak
solutions obtained in the massive case [12]. In the latter case there is
however no need to glue a Schwarzschild solution since the massive solutions
have compact support.

###  5.1 A black hole with one shell

First we compute solutions for the parameter values $k=0$ and $l=1/2$ used in
the proof. We have chosen the mass of the black hole to be $M_{0}=0.25$ so
that $L_{*}=27M_{0}^{2}=27/16\approx 1.69$. First we try to put the shell
close to the photon sphere of the black hole by choosing $L_{0}=1.7$. However,
as Figure 1 shows there is no vacuum region after the first peak (and thus
there is no proper shell) and it is not possible to obtain an asymptotically
flat solution. (In all figures the energy density $\rho$ is displayed on the
vertical axis and the area radius $r$ on the horizontal axis.) The choice
$L_{0}=1.86$ gives on the other hand rise to a vacuum region after the first
peak and an asymptotically flat solution is obtained. The photon sphere is
located at $r=3M_{0}=0.75$, and the inner radius of the shell is located at
$R_{0}=0.92$. The compactness ratio $\Gamma$ of the shell is $0.74$. The shell
is depicted in Figure 2. The radius of the photon sphere preceding the shell
is denoted by $r_{*}$.

If we increase the parameter value $k$, to $k=1$ instead of $k=0$, we have to
increase $L_{0}$ to obtain a proper shell, i.e., the shell must be placed
further away from the photon sphere of the black hole. Indeed, in Figure 3 we
have used the same parameter values as in Figure 2 with the only change that
$k=1$ instead of $k=0$. We see that in this case a proper shell is not
obtained. We have to increase $L_{0}$ to $L_{0}=5.5$ in order to get a proper
shell. This case as shown in Figure 4. The inner radius is located at
$R_{0}=2.04$ and the compactness ratio is $\Gamma=0.80$. Hence, the parameter
$k$ has a considerable influence on how close to the photon sphere the shell
can be placed. We point out that the parameter $l$ has a similar impact but we
have not made a systematic study of the dependence on $k$ and $l$.

If we increase the mass of the black hole the situation changes and the shell
can be placed closer to the photon sphere. Again this depends on the parameter
values of $k$ and $l$, but generally the shells can be placed closer to the
photon sphere when $M_{0}$ is larger. If $M_{0}$ is taken sufficiently large
the shell can be placed arbitrary close to the photon sphere. This results in
a low value of the compactness ratio $\Gamma$. Since $2M_{0}/r=2/3$ at the
location of the photon sphere, we conjecture that $\Gamma>2/3$ for any shell
surrounding the black hole, where low values are attained for shells which are
located very close to the photon sphere. Recall from the proof that as the
inner radius $R_{0}\to\infty$, the compactness ratio $\Gamma\to 8/9$ (a fact
that is confirmed numerically in Figures 6(c), 7(c) and 8(c)), so that morally
it is advantageous to keep $R_{0}$ small in order to get a low value of
$\Gamma$.

We choose $M_{0}=2.0$, which implies that $L_{*}=27M_{0}^{2}=108$. The photon
sphere is located at $r=3M_{0}=6.0$. We find that the choice $L_{0}=108.1$
gives rise to a proper shell which is located very close to the photon sphere,
at $R_{0}=6.1$. The shell is depicted in Figure 5, where $k=0$ and $l=1/2$.
The compactness ratio $\Gamma$ for the shell is merely $0.68$. If $k$ is
increased to $k=1$, then we need to increase $L_{0}$, to $L_{0}=113$, in order
to obtain a proper shell, which results in a shell with an inner radius of
$R_{0}=6.9$, and a compactness ratio $\Gamma=0.73$.

Figure 1: Not a proper shell ($k=0,l=1/2$). Figure 2: A proper shell
($k=0,l=1/2$). $L_{0}=1.86,\Gamma=0.74$ and $r_{*}=0.75$. Figure 3: Not a
proper shell ($k=1,l=1/2$). Figure 4: A proper shell ($k=1,l=1/2$).
$L_{0}=5.5,\Gamma=0.8$ and $r_{*}=0.75$. Figure 5: A proper shell
($k=0,l=1/2$). $L_{0}=108.1,\Gamma=0.68$ and $r_{*}=6.0$.

###  5.2 The thickness of the shells

We claim in Remark (3.4) that the thickness of the shells, $R_{1}-R_{0}$,
depend on the parameters $k$ and $l$ as

$R_{1}-R_{0}\sim R_{0}^{\frac{q-2-2l}{q}},$ (5.41)

where $q=k+l+5/2$. In Figure (6), (7) and (8) we have computed shells for
$L_{0}=15$, $L_{0}=75$ and $L_{0}=375$ in the three cases $k=0,l=1/2$;
$k=1,l=1/2$ and $k=0,l=3/2$. The quantity

$\tau:=\frac{q-2-2l}{q},$

takes the following values in these cases, $0,1/4$ and $-1/4$, respectively.
When $\tau=0$, the width of the shell is independent on $R_{0}$ and this is
confirmed in Figure (6) where the thickness $R_{1}-R_{0}\approx 0.29$ for any
of these shells. This is the case we considered in the proof. In the second
case $\tau=1/4$, and the width of the shell thus increases as $R_{0}$ grows.
We compute the ratio

$\frac{R_{1}-R_{0}}{R_{0}^{\tau}},$ (5.42)

for the three shells in Figure (7) and we find that it is approximately
constant, roughly $0.67$. In the last case $\tau=-1/4$, and the width
decreases as $R_{0}$ increases. This can be seen in Figure (8) where the ratio
(5.42) is roughly $0.46$ for each shell. Hence, we find numerical support for
the claim (5.41). Moreover, $\Gamma$ is computed in the three cases where
$L_{0}=375$, to confirm that $\Gamma\to\frac{8}{9}$ as $R_{0}\to\infty$, cf.
Figures 6(c), 7(c) and 8(c).

(a) $L_{0}=15$

(b) $L_{0}=75$

(c) $L_{0}=375,\Gamma=0.883$

Figure 6: Three shells with constant thickness ($\tau=0$)

(a) $L_{0}=15$

(b) $L_{0}=75$

(c) $L_{0}=375,\Gamma=0.875$

Figure 7: Three shells with growing thickness ($\tau=1/4$)

(a) $L_{0}=15$

(b) $L_{0}=75$

(c) $L_{0}=375,\Gamma=0.884$

Figure 8: Three shells with decreasing thickness ($\tau=-1/4$)

###  5.3 A black hole with several nested shells

Here we construct solutions where the black hole is surrounded by several
shells. We choose the black hole mass $M_{0}=1.0$ so that $L_{*}=27.0$. We put
a shell close to the photon sphere by the choice $L_{0}=28.0$. We obtain a
solution with one shell surrounding the black hole with ADM mass of $1.3$.
This implies that $L_{*}=45.6$ and we choose $L_{0}=48$ for the second shell.
The solution we obtain is depicted in Figure 9 where we have taken $k=0$ and
$l=1/2$. The radius of the photon sphere surrounding the black hole is
$r_{*}=3.0$ and the radius of the photon sphere surrounding the first shell is
located at $r_{*}=4.2$. Similar as above, if we increase the parameter values
$k$, then we have to take larger values of $L_{0}$ to obtain proper shells.
Clearly, the procedure can be continued and an arbitrary number of shells can
be constructed which surround the black hole.

Next, we again start with a black hole but we only use one ansatz function for
the shells, in particular we only solve equation (3.11) once and $L_{0}$ is
fixed. We then obtain a solution consisting of several nested shells,
separated by vacuum. In Figure 10 such a solution is depicted where
$M_{0}=1.3$, $L_{0}=48$, $k=0$ and $l=1/2$. After the fifth peak there is no
longer a vacuum region separating the neighboring shells (which is difficult
to see in the picture since $\rho$ is very small) and a Schwarzschild solution
has to be glued after the fourth peak in order to obtain an asymptotically
flat solution. The solution obtained in this way is very similar to the multi-
peak solutions obtained in the massive case [12]. In the latter case there is
however no need to glue a Schwarzschild solution since the massive solutions
have compact support. We also include a few cases with other choices of the
parameters $k$ and $l$. In Figure 11 we have used the same parameters as in
Figure 10 with the only change that $k=1$. In this case there are no proper
shells and an asymptotically flat solution is not possible to obtain. By
increasing $L_{0}$ proper shells are obtained as depicted in Figure 12.
Finally we consider the case $l=3/2$ but otherwise the same parameters as in
Figure 10. In this case three proper shells separated by vacuum are obtained
as is depicted in Figure 13.

Acknowledgement: The author would like to thank Gerhard Rein for useful
discussions.

Figure 9: $L_{0}=28.0$ (shell one) and $L_{0}=48$ (shell two). Figure 10:
Multi-shells with one ansatz. $L_{0}=48$, $M_{0}=1.3$, $k=0$ and $l=1/2$.
Figure 11: No proper shells. $L_{0}=48$, $M_{0}=1.3$, $k=0$ and $l=1/2$.
Figure 12: Multi-shells with one ansatz. $L_{0}=80$, $M_{0}=1.3$, $k=1$ and
$l=1/2$. Figure 13: Multi-shells with one ansatz. $L_{0}=48$, $M_{0}=1.3$,
$k=0$ and $l=3/2$.

## References

  * [1] A. Akbarian, M. Choptuik Critical collapse in the spherically symmetric Einstein-Vlasov model. Phys. Rev. D 90, 104023 (2014).
  * [2] E. Ames, H. Andréasson and A. Logg On axisymmetric and stationary solutions of the self-gravitating Vlasov system. Class. Quant. Grav. 33, 155008 (2016).
  * [3] E. Ames, H. Andréasson and A. Logg Cosmic string and black hole limits of toroidal Vlasov bodies in general relativity. Phys. Rev. D 99, 024012 (2019).
  * [4] L. Andersson, P. Blue and J. Joudioux Hidden symmetries and decay for the Vlasov equation on the Kerr spacetime. Comm. Partial Differential Equations 43, 47–65 (2018).
  * [5] H. Andréasson On the Buchdahl inequality for spherically symmetric static shells. Commun. Math. Phys. 274, 399–408 (2007).
  * [6] H. Andréasson On static shells and the Buchdahl inequality for the spherically symmetric Einstein-Vlasov system. Commun. Math. Phys. 274, 409–425 (2007).
  * [7] H. Andréasson Sharp bounds on $2m/r$ of general spherically symmetric static objects. J. Differential Equations 245, 2243–2266 (2008).
  * [8] H. Andréasson The Einstein-Vlasov system/Kinetic theory. Liv. Rev. Rel. 14 (2011).
  * [9] H. Andréasson, D. Fajman, M. Thaller Static Solutions to the Einstein–Vlasov System with a Nonvanishing Cosmological Constant. SIAM J. Math. Anal. 47, 2657–2688 (2015).
  * [10] H. Andréasson, D. Fajman and M. Thaller Models for self-gravitating photon shells and geons. Ann. Henri Poincaré 18, 681–705 (2017).
  * [11] H. Andréasson, M. Kunze, G. Rein Rotating, stationary, axially symmetric spacetimes with collisionless matter. Commun. Math. Phys. 329, 787–808 (2014).
  * [12] H. Andréasson, G. Rein On the steady states of the spherically symmetric Einstein-Vlasov system. Classical Quantum Gravity 24, 1809–1832 (2007).
  * [13] C. Cederbaum and G.j. Galloway Uniqueness of photon spheres via positive mass rigidity. Comm. Anal. Geom. 25, 303–32 (2017).
  * [14] H.A. Buchdahl General relativistic fluid spheres. Phys. Rev.  116, 1027–1034 (1959).
  * [15] J. Frauendiener, C. Hoenselaers and W. Konrad A shell around a black hole. Classical Quantum Gravity 7, 585–587 (1990).
  * [16] F. Finster, J. Smoller and S.T. Yau Absence of static, spherically symmetric black hole solutions for Einstein-Dirac-Yang/Mills equations with complete fermion shells. Adv. Theor. Math. Phys. 4, 1231–1257 (2000).
  * [17] C. Gundlach Einstein-Vlasov system in spherical symmetry: reduction of the equations of motion and classification of single-shell static solutions in the limit of massless particles. Phys. Rev. D 94, 124046 (2016).
  * [18] F.E. Jabiri Static spherically symmetric Einstein-Vlasov bifurcations of the Schwarzschild spacetime. arXiv:2001.08645.
  * [19] T. Makino On spherically symmetric stellar models in general relativity. J. Math. Kyoto Univ. 38, 55–69 (1998).
  * [20] T. Ramming, G. Rein Spherically symmetric equilibria for self-gravitating kinetic or fluid models in the non-relativistic and relativistic case – a simple proof for finite extension. SIAM J. Math. Anal. 45, 900–914 (2013).
  * [21] G. Rein Static solutions of the spherically symmetric Vlasov-Einstein system. Math. Proc. Camb. Phil. Soc. 115, 559–570 (1994).
  * [22] G. Rein, The Vlasov-Einstein system with surface symmetry, Habilitationsschrift, Munich (1995).
  * [23] G. Rein, A. Rendall Smooth static solutions of the spherically symmetric Vlasov-Einstein system. Ann. Inst. H. Poincaré Phys. Théor. 59, 383–397 (1993).
  * [24] A.D. Rendall The initial singularity in solutions of the Einstein–Vlasov system of Bianchi type I. J. Math. Phys. 37, 438–451 (1996).
  * [25] A.D. Rendall and J.J.L. Velazquez Veiled singularities for the spherically symmetric massless Einstein-Vlasov system. Ann. Henri Poincaré 18, 3565–3631 (2017).