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 [25]. 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 [22] 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 [2, 3, 9] for more recent results.

By relaxing the condition of a regular centre the case with a Schwarzschild black hole was considered in [23], where the existence of massive static shells of Vlasov matter surrounding a black hole was shown. A different method leading to a similar result has been more recently given by Jabiri [20]. For a fluid, the first result of a massive static shell surrounding a black hole was obtained in [18]. If the matter model originates from quantum mechanics, similar results need not be true as is for instance shown in [17], where the absence of static black hole solutions is shown for the Einstein–Dirac–Yang/Mills equations. It is argued in [17] 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 [26] 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 [19]. 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 [13] 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, 19], 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 [10, 23]. 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. Sect. 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 [14] which shows that only one photon sphere can appear in a static spacetime. However, the result in [14] 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\rightarrow 8/9\) in both the limits \(R_0\rightarrow 0\) and \(R_0\rightarrow \infty \). (It should be pointed out that in the former case \(L_0\) is fixed and \(R_0\) becomes small by increasing \(E_0\).)

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 Sect. 2 we introduce the massless static Einstein–Vlasov system. In Sect. 3 we formulate the main results and in Sect. 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

$$\begin{aligned} 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}), \end{aligned}$$

where \(r\ge 0,\,\theta \in [0,\pi ],\,\varphi \in [0,2\pi ]\) and \(t\in {\mathbb {R}}\). Asymptotic flatness is expressed by the boundary conditions

$$\begin{aligned} \lim _{r\rightarrow \infty }\lambda (r)=\lim _{r\rightarrow \infty }\mu (r) =0. \end{aligned}$$

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

$$\begin{aligned} \displaystyle e^{-2\lambda }(2r\lambda _{r}-1)+1= & {} 8\pi r^2\rho , \end{aligned}$$
(2.1)
$$\begin{aligned} \displaystyle e^{-2\lambda }(2r\mu _{r}+1)-1= & {} 8\pi r^2 p, \end{aligned}$$
(2.2)
$$\begin{aligned} \displaystyle \mu _{rr}+(\mu _{r}-\lambda _{r})\left( \mu _{r}+\frac{1}{r}\right)= & {} 8\pi p_T e^{2\lambda }, \end{aligned}$$
(2.3)

together with the static Vlasov equation

$$\begin{aligned} \frac{w}{\varepsilon }\partial _{r}f -\left( \mu _{r} \varepsilon - \frac{L}{r^3\varepsilon }\right) \partial _{w}f=0, \end{aligned}$$
(2.4)

where

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

The matter quantities are defined by

$$\begin{aligned} \rho (r)= & {} \frac{\pi }{r^{2}} \int _{-\infty }^{\infty }\int _{0}^{\infty }\varepsilon (r,w,L) f(r,w,L)\;\mathrm{d}L\mathrm{d}w, \end{aligned}$$
(2.5)
$$\begin{aligned} p(r)= & {} \frac{\pi }{r^{2}}\int _{-\infty }^{\infty }\int _{0}^{\infty } \frac{w^{2}}{\varepsilon (r,w,L)}f(r,w,L)\;\mathrm{d}L\mathrm{d}w, \end{aligned}$$
(2.6)
$$\begin{aligned} p_T(r)= & {} \frac{\pi }{2r^{4}}\int _{-\infty }^{\infty }\int _{0}^{\infty }\frac{L}{\varepsilon (r,w,L)}f(r,w,L)\; \mathrm{d}L\mathrm{d}w. \end{aligned}$$
(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

$$\begin{aligned} E=e^{\mu }\varepsilon , \end{aligned}$$

then the ansatz

$$\begin{aligned} f(r,w,L)=\Phi (E,L), \end{aligned}$$
(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

$$\begin{aligned} \Phi (E,L)=(E_0-E)^k_{+}(L-L_0)_{+}^l, \end{aligned}$$
(2.9)

where \(l\ge 1/2,\,k\ge 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

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

is large; roughly \(\Gamma \ge \frac{4}{5}\). Here m is the Hawking mass defined for \(r\ge 2M_0\) by

$$\begin{aligned} m(r)=M_0+\int _{2M_0}^r s^2\rho (s)\, \mathrm{d}s, \;\; r\ge 2M_0. \end{aligned}$$

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. Sect. 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 indicate that the required lower bound is larger in that case, cf. [1, 10, 19].

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. Eq. (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 [23]. 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

$$\begin{aligned} e^{2\mu (r)}=1-\frac{2M_0}{r}. \end{aligned}$$

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

$$\begin{aligned} \left( 1-\frac{2M_0}{r}\right) \frac{L_0}{r^2}= E_0^2=1. \end{aligned}$$
(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. [19, 27]. 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

$$\begin{aligned} \frac{L_0}{27M_0^2}, \end{aligned}$$

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

$$\begin{aligned} L_0>27M_0^2=:L_*. \end{aligned}$$

Equation (3.11) then has three real roots and we denote the largest root by \(R_0\) where \(R_0>3M_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

$$\begin{aligned} 1-\frac{C}{\sqrt{L_0}}\le \frac{R_0}{\sqrt{L_0}}\le 1 \text{ as } L_0\rightarrow \infty . \end{aligned}$$

Proof

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

$$\begin{aligned} r^3-L_0r+2M_0L_0=0, \end{aligned}$$

which we write as \(r^3+pr+q=0\) where \(p=-L_0\) and \(q=2M_0L_0\). Set

$$\begin{aligned} D=\left( \frac{p}{3}\right) ^3+\left( \frac{q}{2}\right) ^2=L_0^2\left( -\frac{L_0}{27}+M_0^2\right) . \end{aligned}$$

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

$$\begin{aligned} r_1=u+v,\;\; r_{2,3}=\frac{u+v}{2}\pm \frac{u-v}{2}i\sqrt{3}, \end{aligned}$$

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

$$\begin{aligned} R_0:=r_1=u+v=2\left( \frac{q^2}{4}+|D|\right) ^{1/6}\cos {\frac{\beta }{3}}, \end{aligned}$$

where

$$\begin{aligned} \beta = \frac{\pi }{2}+\arctan {\frac{(q/2)}{\sqrt{|D|}}}. \end{aligned}$$

We have

$$\begin{aligned} 0\le \frac{(q/2)}{\sqrt{|D|}}=\le \frac{M_0}{\sqrt{L_0}\sqrt{\frac{1}{27}-\frac{M_0^2}{L_0}}}\le \frac{C}{\sqrt{L_0}}, \end{aligned}$$

where we used that \(L_0-L_*>0\) since we consider large \(L_0\). Hence there exists a constant \(C>0\) such that

$$\begin{aligned} \frac{\pi }{2}\le \beta \le \frac{\pi }{2}+\frac{C}{\sqrt{L_0}}, \end{aligned}$$

for large \(L_0\). This implies that

$$\begin{aligned} \frac{\sqrt{3}}{2}-\frac{C}{\sqrt{L_0}}\le \cos {\frac{\beta }{3}}\le \frac{\sqrt{3}}{2}, \end{aligned}$$

for some positive constant C. Moreover, we have

$$\begin{aligned} 2\left( \frac{q^2}{4}+|D|\right) ^{1/6}=2\left( M_0^2L_0^2+L_0^2\left( \frac{L_0}{27}-M_0^2\right) \right) ^{1/6}=\frac{2\sqrt{L_0}}{\sqrt{3}}. \end{aligned}$$

Hence we conclude that for large \(L_0\)

$$\begin{aligned} 1-\frac{C}{\sqrt{L_0}}\le \frac{R_0}{\sqrt{L_0}}\le 1. \end{aligned}$$

\(\square \)

We now consider the Einstein–Vlasov system on the domain \(r\ge R_0\) and we prescribe data for \(\mu \) and \(\lambda \) at \(r=R_0\) by letting

$$\begin{aligned} e^{2\mu (R_0)}=e^{-2\lambda (R_0)}=1-\frac{2M_0}{R_0}. \end{aligned}$$

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. [10] 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

$$\begin{aligned} \frac{R_1}{R_0}\rightarrow 1 \text{ as } R_0\rightarrow \infty . \end{aligned}$$
(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

$$\begin{aligned} f(r,w,L)=\Phi (E,L), \end{aligned}$$
(3.13)

where \(\Phi =0\) whenever \(E>E_0\), the matter quantities \(\rho \) and p become functionals of \(\mu \), and we have

$$\begin{aligned} \rho= & {} \frac{2\pi }{r^2}\int _{\frac{\sqrt{L_0}}{r}}^{E_0e^{-\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}}}\mathrm{d}L\mathrm{d}s, \end{aligned}$$
(3.14)
$$\begin{aligned} p= & {} \frac{2\pi }{r^2}\int _{\frac{\sqrt{L_0}}{r}}^{E_0e^{-\mu }}\int _{L_0}^{r^2(s^2-1)}\Phi (e^{\mu }s,L)\sqrt{s^2-1-\frac{L}{r^2}}\mathrm{d}L\mathrm{d}s. \end{aligned}$$
(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

$$\begin{aligned} \gamma =-\mu -\frac{1}{2}\log {\frac{L_0}{r^2}}, \end{aligned}$$

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 \ge 0\)

$$\begin{aligned} \rho= & {} \pi ^{1}_{k,l}r^{2l}\left( \frac{L_0}{r^2}\right) ^{l+2}(e^{\gamma }-1)^{l+k+3/2}P_{l+5/2-k}(e^{\gamma }), \end{aligned}$$
(3.17)
$$\begin{aligned} p= & {} \pi ^{2}_{k,l}r^{2l}\left( \frac{L_0}{r^2}\right) ^{l+2}(e^{\gamma }-1)^{l+k+5/2}P_{l+3/2-k}(e^{\gamma }), \end{aligned}$$
(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\le C\), independently of \(R_0\) as shown below, whereas for other values of the parameters we claim that

$$\begin{aligned} R_1-R_0\sim R_0^{\frac{q-2-2l}{q}}, \end{aligned}$$

where \(q=k+l+5/2\), cf. Sect. 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\))

$$\begin{aligned} \rho (r)=\pi ^2 r\left( \frac{L_0}{r^2}\right) ^{5/2}(t-1)^2\left( \frac{3t^3+6t^2+4t+2}{15}\right) , \end{aligned}$$
(3.19)

and

$$\begin{aligned} p(r)=\pi ^2 r\left( \frac{L_0}{r^2}\right) ^{5/2}(t-1)^3\left( \frac{3t^2+9t+8}{60}\right) . \end{aligned}$$
(3.20)

Let \(R_0\) be large and define

$$\begin{aligned} \delta :=\left( \frac{1}{20\pi ^3}\right) ^{1/3}. \end{aligned}$$
(3.21)

The argument below will be carried out in the interval \(I:=[R_0,R_1]\) where \(R_1\le 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\ge 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\rightarrow \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. Eqs. (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 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 \rightarrow 8/9\) in the extreme limits, i.e. when \(R_0\rightarrow 0\) as in [10] or when \(R_0\rightarrow \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 Sect. 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\) as described above. For \(R_0\) sufficiently large, there exists \(R_1\) such that \(R_0<R_1\le R_0+50\delta \), and 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 \rightarrow \frac{8}{9}\) as \(R_0\rightarrow \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\rightarrow \infty \) as \(L_0 \rightarrow \infty \) in view of Lemma 3.2. Next we note that \(R_1^*/R_0\le 1+C/R_0\). We remark that here and below C denotes a constant which might depend on \(M_0\) and which might change from line to line. Hence, by Lemma 3.2 we have for any \(r\in [R_0,R_1^*]\)

$$\begin{aligned} 1-\alpha (R_0)\le {\frac{L_0}{r^2}}\le 1+\alpha (R_0), \end{aligned}$$
(4.22)

where we have introduced the notation \(\alpha (R_0)\) where \(\alpha (R_0)\ge 0\) and \(\alpha (R_0)\rightarrow 0\) as \(R_0\rightarrow \infty \). By the definition of \(R_0\) we have that \(\gamma (R_0)=0\). Moreover,

$$\begin{aligned} \gamma '(r)=-\mu '(r)+\frac{1}{r}. \end{aligned}$$
(4.23)

Now, since \(\rho =p=0\) and \(m(r)=M_0\) for \(r\le R_0\), we have that that \(\mu '(R_0)<1/R_0\). Here we used that \(M_0/R_0<1/3\) and that

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

The last relation is a consequence of the Einstein equation (2.1). This implies that \(\gamma '(R_0)>0\). Thus \(\gamma (r)>0\) in a right neighbourhood 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 '(r)\ge 0\),

$$\begin{aligned} \gamma (r)=\gamma (r)-\gamma (R_0)\le \log \frac{r}{R_0}=\log \left( 1+\frac{r-R_0}{R_0}\right) , \end{aligned}$$

which implies that \(0\le \gamma \le \alpha (R_0)\) on I if \(R_0\) is sufficiently large. Hence for \(r\in I\)

$$\begin{aligned} \gamma (r)\le e^{\gamma (r)}-1\le \gamma (r)(1+\alpha (R_0)). \end{aligned}$$

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

$$\begin{aligned} 1\le \frac{3t^3+6t^2+4t+2}{15}\le 1+\alpha (R_0), \end{aligned}$$

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

$$\begin{aligned} \pi ^2 r\gamma ^2(r)(1-\alpha (R_0))\le \rho (r)\le \pi ^2 r\gamma ^2(r)(1+\alpha (R_0)), \end{aligned}$$

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

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

and

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

for \(r\in I\). \(\square \)

Lemma 4.1

Let \(\delta \) be as above. Then

$$\begin{aligned} \gamma '(r)\ge \frac{1}{2r}\; \text{ for } r\in [R_0,R_0+\delta ]=:I_1. \end{aligned}$$

Proof

We have by the mean value theorem that for any \(\sigma \le \delta \)

$$\begin{aligned} \gamma (R_0+\sigma )=\gamma (R_0+\sigma )-\gamma (R_0)=\sigma \gamma '(\xi )\le \frac{\delta }{R_0}, \end{aligned}$$

where \(\xi \in [R_0,R_0+\sigma ]\), since \(\gamma (r)\le 1/r\). Hence,

$$\begin{aligned} \rho (r)\le \pi ^2 r\frac{\delta ^2}{R_0^2}\; \text{ on } I_1. \end{aligned}$$

We get for \(\sigma \le \delta \)

$$\begin{aligned} m(R_0+\sigma )\le M_0+\frac{4\pi ^3\delta ^2}{R_0^2}\int _{R_0}^{R_0+\sigma }\eta ^3 \, \mathrm{d}\eta \le M_0+\frac{4\pi ^3\delta ^2(R_0+\sigma )^3}{R_0^2}\sigma . \end{aligned}$$

By taking \(R_0\) large we thus obtain

$$\begin{aligned} \frac{m(R_0+\sigma )}{R_0+\sigma }\le 4\pi ^3\delta ^2\sigma +\alpha (R_0)\le 4\pi ^3\delta ^3+\alpha (R_0)\le \frac{1}{5} +\alpha (R_0). \end{aligned}$$

This implies that for \(r\in I_1\)

$$\begin{aligned} \frac{m(r)}{r}e^{2\lambda (r)}\le \frac{1}{3}+\alpha (R_0). \end{aligned}$$

Moreover, we have for \(r\in I_1\) that

$$\begin{aligned} 4\pi r^2 p(r)e^{2\lambda (r)}\le 4\pi r^2 p(r)\left( \frac{5}{3}+\alpha (R_0)\right) \le \frac{4\pi ^3\delta ^3}{3}\left( \frac{5}{3}+\alpha (R_0)\right) \le \frac{1}{9}+\alpha (R_0). \end{aligned}$$

Hence, for \(r\in I_1\) we get for sufficiently large \(R_0\)

$$\begin{aligned} \mu '(r)\le \left( \frac{1}{3}+\frac{1}{9}+\alpha (R_0)\right) \frac{1}{r}\le \frac{1}{2r}. \end{aligned}$$

This completes the proof of the lemma. \(\square \)

The lemma implies that for \(\sigma \le \delta \)

$$\begin{aligned} \gamma (R_0+\sigma )\ge \gamma (R_0)+\sigma \inf _{0\le s\le \sigma } \gamma '(R_0+s)\ge \frac{\sigma }{2(R_0+\sigma )}. \end{aligned}$$

Let

$$\begin{aligned} \sigma _*:=\delta , \end{aligned}$$

and define

$$\begin{aligned} \gamma _*=\frac{\sigma _*}{2(R_0+\sigma _*)}. \end{aligned}$$

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

$$\begin{aligned} \gamma (R_0+\sigma _*)\ge \gamma _*. \end{aligned}$$

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\le R_0+11\delta \).

Proof

Since \(\gamma (R_0+\sigma _{*})\ge \gamma _{*},\) and since Lemma 4.1 gives that \(\gamma '(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 \ge \gamma _{*}\) on the interval. We will show that there in fact \(\gamma (R_0+\sigma _*+\Delta )=\gamma _*\). We have from (4.24)

$$\begin{aligned} m(R_0+\sigma _*+\Delta )\ge & {} 4\pi ^3\int _{R_0+\sigma _{*}}^{R_0+\sigma _{*}+\Delta }\frac{\sigma _{*}^2}{4(R_0+\sigma _{*})^2}r^3\, \mathrm{d}r\nonumber \\\ge & {} \pi ^3\sigma _*^2(R_0+\sigma _*)\Delta . \end{aligned}$$
(4.26)

Hence,

$$\begin{aligned} \frac{m(R_0+\sigma _{*}+\Delta )}{R_0+\sigma _{*}+\Delta }\ge \frac{\pi ^3\sigma _{*}^2(R_0+\sigma _*)}{R_0+\sigma _{*}+\Delta }\Delta \ge \frac{9}{10} \pi ^3\sigma _{*}^2\Delta , \end{aligned}$$

where we used that \(R_0\ge 100\delta \), so that \((R_0+\sigma _*)/(R_0+\sigma _*+\Delta )\ge 9/10\), since we consider large \(R_0\). Since for \(r\ge R_0\) necessarily

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

we conclude that there is a \(\Delta \) such that

$$\begin{aligned} \Delta \le \frac{40}{81\pi ^3 \sigma _*^2}=\frac{40}{81\pi ^3 \delta ^2}=\frac{40\cdot 20^{2/3}}{81\pi }\le 10\delta , \end{aligned}$$

with the property that \(\gamma (R_0+\sigma _*+\Delta )=\gamma _*\), since otherwise we obtain a contradiction. \(\square \)

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

$$\begin{aligned} \frac{m(r_2)}{r_2}\ge \frac{2}{5}. \end{aligned}$$

Proof

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

$$\begin{aligned} \left( \frac{m(r)}{r^2}+4\pi rp(r)\right) 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)\mathrm{d}\eta .\nonumber \\ \end{aligned}$$
(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–Volkoff equation

$$\begin{aligned} p_r=-\mu _r(p+\rho )-\frac{2}{r}(p-p_T), \end{aligned}$$

cf. [5] for details. 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

$$\begin{aligned} \frac{m(r_2)}{r_2}e^{(\mu +\lambda )(r_2)}\ge & {} \frac{8\pi }{r_2}\int _{R_0}^{r_2} \eta ^2 e^{\mu +\lambda }\rho \,\mathrm{d}\eta -4\pi r_2^2pe^{(\mu +\lambda )(r_2)}. \end{aligned}$$
(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

$$\begin{aligned} \int _{R_0}^{r_2}4\pi \eta \rho e^{\lambda }\,\mathrm{d}\eta= & {} \int _{R_0}^{r_2}\left( -\frac{\mathrm{d}}{\mathrm{d}r}e^{-\lambda }\right) \mathrm{d}\eta +\int _{R_0}^{r_2}\frac{me^{\lambda }}{\eta ^2}\mathrm{d}\eta \nonumber \\\ge & {} \sqrt{1-\frac{2M_0}{R_0}}-\sqrt{1-\frac{2m(r_2)}{r_2}}\nonumber \\= & {} 1-\sqrt{1-\frac{2m(r_2)}{r_2}}-\left( 1-\sqrt{1-\frac{2M_0}{R_0}}\right) \nonumber \\= & {} \frac{2m(r_2)}{r_2\left( 1+\sqrt{1-\frac{2m(r_2)}{r_2}}\right) }-\alpha (R_0). \end{aligned}$$
(4.29)

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

$$\begin{aligned} \frac{m(r_2)}{r_2}e^{(\mu +\lambda )(r_2)}\ge & {} \frac{8\pi }{r_2}\int _{R_0}^{r_2}\eta ^2 e^{\mu +\lambda }\rho \,\mathrm{d}\eta -4\pi r_2^2p(r_2)e^{(\mu +\lambda )(r_2)}\\\ge & {} 8\pi e^{\mu (R_0)}\frac{R_0}{r_2} \int _{R_0}^{r_2}\eta e^{\lambda }\rho \,\mathrm{d}\eta -4\pi r_2^2 p(r_2)e^{(\mu +\lambda )(r_2)}. \end{aligned}$$

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

$$\begin{aligned} \frac{m(r_2)}{r_2}e^{(\mu +\lambda )(r_2)}\ge & {} \frac{e^{\mu (R_0)}R_0}{r_2}\frac{4m(r_2)}{r_2(1+\sqrt{1-\frac{2m(r_2)}{r_2}})}\nonumber \\&-Pe^{(\mu +\lambda )(r_2)}-\alpha (R_0). \end{aligned}$$
(4.30)

Letting

$$\begin{aligned} Y:=\frac{m(r_2)}{r_2}, \end{aligned}$$

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

$$\begin{aligned} Y\ge 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)}. \end{aligned}$$
(4.31)

Next we show a couple of properties of the solutions that we need to proceed with the argument. From Lemma 4.3 we have that \(\gamma \) approaches \(\gamma _{*}\) from above and therefore \(\gamma '(r_2)\le 0,\) which implies that

$$\begin{aligned} \mu '(r_2)\ge \frac{1}{r_2}, \end{aligned}$$
(4.32)

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

$$\begin{aligned} 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}. \end{aligned}$$

Now recall that

$$\begin{aligned} \mu '=\left( \frac{m}{r^2}+4\pi rp\right) e^{2\lambda }. \end{aligned}$$

We will show that \(Y=m(r_2)/r_2\ge 1/4\). Assume the contrary, i.e. assume

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

then

$$\begin{aligned} e^{2\lambda (r_2)}=\frac{1}{1-2Y}<2, \end{aligned}$$

and thus

$$\begin{aligned} 4\pi r_2^2p(r_2)e^{2\lambda (r_2)}\le \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}\le \frac{1}{50}. \end{aligned}$$
(4.34)

Here we used that \(\delta ^3=1/20\pi ^3\) and that \(r_2\le R_0+11\delta \) so that

$$\begin{aligned} \frac{r_2^3}{(R_0+\delta )^3}\le \frac{6}{5} \text{ for } R_0 \text{ large }. \end{aligned}$$

Hence if (4.33) holds, then

$$\begin{aligned} \mu '(r_2)=\left( \frac{m(r_2)}{r_2^2}+4\pi r_2p(r_2)\right) e^{2\lambda (r_2)}\le \left( \frac{1}{2}+\frac{1}{50}\right) \frac{1}{r_2}, \end{aligned}$$

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

$$\begin{aligned} Y\ge \frac{1}{4}. \end{aligned}$$

Let us next consider the difference \(\mu (R_0)-\mu (r_2).\) We have

$$\begin{aligned} \mu (r_2)-\mu (R_0)=\int _{R_0}^{r_2} \left( \frac{m}{\eta ^2}+4\pi \eta p\right) e^{2\lambda }\,\mathrm{d}\eta . \end{aligned}$$

Since \(\gamma '(r)\le 1/r\) we have for \(R_0\le r\le r_2\le R_0+11\delta \)

$$\begin{aligned} \gamma (r)\le \log {\frac{r}{R_0}}\le \frac{11\delta }{R_0}, \end{aligned}$$

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

$$\begin{aligned} 4\pi r^2 p(r)e^{2\lambda (r)}\le 12\pi ^3 r^3 \gamma ^3(r)\le \frac{3}{5}\big (\frac{r}{R_0}\big )^3 11^3\le C. \end{aligned}$$

Hence, since \(m/r\le 4/9\) we get

$$\begin{aligned} \int _{R_0}^{r_2} \left( \frac{m}{\eta ^2}+4\pi \eta p\right) e^{2\lambda }\,\mathrm{d}\eta \le \int _{R_0}^{r_2} \frac{C}{\eta }\,\mathrm{d}\eta \le C\log {\frac{r_2}{R_0}}\le \frac{C}{R_0}. \end{aligned}$$

As a consequence we obtain that

$$\begin{aligned} 1-\alpha (R_0)\le e^{\mu (R_0)-\mu (r_2)}\left( \frac{R_0}{r_2}\right) \le 1+\alpha (R_0). \end{aligned}$$

From this estimate we also have that \(e^{-(\mu +\lambda )(r_2)}\) is bounded independently of \(R_0\) and using that \(Y\ge 1/4\) we can therefore write (4.31) as

$$\begin{aligned} 1\ge \frac{4\sqrt{1-2Y}}{1+\sqrt{1-2Y}}-4P-\alpha (R_0). \end{aligned}$$
(4.35)

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

$$\begin{aligned} 1+\sqrt{1-2Y}, \end{aligned}$$

and then squaring both sides we obtain

$$\begin{aligned} Y\ge \frac{4}{9}-\frac{B\sqrt{1-2Y}(1+\sqrt{1-2Y})}{3}, \end{aligned}$$
(4.36)

where we dropped the term involving \(B^2\) due to sign. Using that \(Y\ge 1/4\) we estimate

$$\begin{aligned} \frac{\sqrt{1-2Y}(1+\sqrt{1-2Y})}{3}\le \frac{1}{3\sqrt{2}}\left( 1+\frac{1}{\sqrt{2}}\right) \le \frac{9}{20}. \end{aligned}$$

Taking \(R_0\) large we derive from (4.36) the inequality

$$\begin{aligned} Y\ge \frac{4}{9}-2P. \end{aligned}$$

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

$$\begin{aligned} Y\ge \frac{4}{9}-\frac{1}{50}\ge \frac{2}{5}. \end{aligned}$$

\(\square \)

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\rightarrow 0\), as in [6], the pressure term \(P\rightarrow 0\) and one can conclude from the argument in the lemma above that as \(R_0\rightarrow 0\) the compactness ratio \(\Gamma \rightarrow 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\rightarrow 1\) as \(R_0\rightarrow \infty \), then we can use the argument above, with \(r_2\) replaced by \(R_1\), to show that the compactness ratio \(\Gamma \rightarrow \frac{8}{9}\) as \(R_0\rightarrow \infty \). This shows the last claim in Theorem 3.10.

Inspired by an idea of T. Makino introduced in [21], we show that \(\gamma \) necessarily must vanish close to the point \(r_2\) if \(R_0\) is sufficiently large. Let

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

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

$$\begin{aligned} rx'=\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)}. \end{aligned}$$

In our case \(r>R_0\) and \(\gamma >0\) and we will show that \(\gamma (r)=0\) for some \(R_1\in [R_0,R_1^*]\). Since \(\gamma >0\) and \(\rho ,p\ge 0\) the first term and the last term can be dropped and we have

$$\begin{aligned} rx'\ge -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 }. \end{aligned}$$
(4.37)

Take \(R_0\) sufficiently large so that \(m(r_2)/r_2\ge 2/5\) by Lemma 4.4. Let \(r\in [r_2,16r_2/15],\) then since m is increasing in r we get

$$\begin{aligned} \frac{m(r)}{r}\ge \frac{m(r_2)}{r}=\frac{r_2}{r}\frac{m(r_2)}{r_2}\ge \frac{15}{16}\cdot \frac{2}{5}=\frac{3}{8}. \end{aligned}$$

Now by the definition of x it follows that

$$\begin{aligned} \frac{x}{1-2\gamma x}=\frac{m}{\gamma r}e^{2\lambda }=\frac{m}{\gamma r(1-2m/r)}\ge \frac{3}{2\gamma }, \text{ when } \frac{m}{r}\ge \frac{3}{8}. \end{aligned}$$

Thus on \([r_2,16r_2/15],\)

$$\begin{aligned} \frac{2x^2}{3(1-2\gamma x)}-\frac{x}{\gamma }\ge 0, \end{aligned}$$

so that on this interval

$$\begin{aligned} rx'\ge \frac{x^2}{3(1-2\gamma x)}-x\ge \frac{4}{3}x^2-x, \end{aligned}$$
(4.38)

where we used that

$$\begin{aligned} \frac{1}{1-2\gamma x}=\frac{1}{1-2m/r}\ge 4 \text{ when } \frac{m}{r}\ge \frac{3}{8}. \end{aligned}$$

The upper bound

$$\begin{aligned} \gamma \le \log \left( 1+\frac{r-R_0}{R_0}\right) \end{aligned}$$

implies that

$$\begin{aligned} x(r_2)=\frac{m(r_2)}{r_2\gamma (r_2)}\ge \frac{2}{5}\frac{R_0}{r_2-R_0}. \end{aligned}$$
(4.39)

In particular

$$\begin{aligned} \frac{x(r_2)}{x(r_2)-3/4}\le \frac{16}{15}, \end{aligned}$$

for \(R_0\) large. Solving (4.38) yields

$$\begin{aligned} x(r)\ge \frac{3}{4}\left( 1-\frac{r(4x(r_2)/3-1)}{4r_2x(r_2)/3}\right) ^{-1}, \text{ on } r\in [r_2,16r_2/15), \end{aligned}$$

and we get that \(x(r)\rightarrow \infty \) as \(r\rightarrow R_1,\) where

$$\begin{aligned} R_1\le r_2\frac{x(r_2)}{x(r_2)-3/4}. \end{aligned}$$
(4.40)

In view of (4.39) we estimate

$$\begin{aligned} R_1\le r_2+\frac{\frac{15}{8}(r_2-R_0)r_2}{R_0-\frac{15}{8}(r_2-R_0)}. \end{aligned}$$

For sufficiently large \(R_0\) we have \(R_0/10\ge \frac{15}{8}(r_2-R_0)\ge 0\), recall that \(r_2-R_0\le 11\delta \), and we obtain

$$\begin{aligned} R_1\le r_2+\frac{\frac{25}{12}(r_2-R_0)r_2}{R_0}\le r_2+3(r_2-R_0)\le R_0+50\delta . \end{aligned}$$

This completes the proof of the theorem. \(\square \)

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 Sect. 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

$$\begin{aligned} R_1-R_0\sim R_0^{\frac{q-2-2l}{q}}, \end{aligned}$$

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 Sect. 5.2.

In Sect. 5.3 we construct solutions with several shells. We do this by using two different strategies. From the set-up in Sect. 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 neighbouring 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 Fig. 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 Fig. 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 Fig. 3 we have used the same parameter values as in Fig. 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 is shown in Fig. 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\rightarrow \infty \), the compactness ratio \(\Gamma \rightarrow 8/9\) (a fact that is confirmed numerically in Figs. 6c, 7c and 8c), 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 Fig. 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\).

Fig. 1
figure 1

Not a proper shell (\(k=0, l=1/2\))

Fig. 2
figure 2

A proper shell (\(k=0, l=1/2\)). \(L_0=1.86, \Gamma =0.74\) and \(r_*=0.75\)

Fig. 3
figure 3

Not a proper shell (\(k=1, l=1/2\))

Fig. 4
figure 4

A proper shell (\(k=1, l=1/2\)). \(L_0=5.5, \Gamma =0.8\) and \(r_*=0.75\)

Fig. 5
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\), depends on the parameters k and l as

$$\begin{aligned} R_1-R_0\sim R_0^{\frac{q-2-2l}{q}}, \end{aligned}$$
(5.41)

where \(q=k+l+5/2\). As shown in Figs. 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

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

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 Fig. 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

$$\begin{aligned} \frac{R_1-R_0}{R_0^{\tau }}, \end{aligned}$$
(5.42)

for the three shells in Fig. 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 Fig. 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 \rightarrow \frac{8}{9}\) as \(R_0\rightarrow \infty \), cf. Figs. 6c, 7c and 8c.

Fig. 6
figure 6

Three shells with constant thickness (\(\tau =0\))

Fig. 7
figure 7

Three shells with growing thickness (\(\tau =1/4\))

Fig. 8
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 Fig. 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.

Remark 5.1

The photon spheres discussed above surround either the black hole or the shells, i.e. they occur in a vacuum region and therefore the radius of the photon sphere satisfy \(r=3m\), where m is the Hawking mass interior of the photon sphere. There are in fact also photon spheres within the shells. The condition for a photon sphere is that \(\mu _r(r)=1/r\), cf. [15]. Since \(\mu _r\ge 1/r\) when \(m/r\ge 1/3\), it is clear that photon spheres also exist within the shells. In this context we also mention that photon spheres are important in astrophysics, cf. e.g. [16, 28].

Fig. 9
figure 9

\(L_0=28.0\) (shell one) and \(L_0=48\) (shell two)

Fig. 10
figure 10

Multi-shells with one ansatz. \(L_0=48\), \(M_0=1.3\), \(k=0\) and \(l=1/2\)

Fig. 11
figure 11

No proper shells. \(L_0=48\), \(M_0=1.3\), \(k=0\) and \(l=1/2\)

Fig. 12
figure 12

Multi-shells with one ansatz. \(L_0=80\), \(M_0=1.3\), \(k=1\) and \(l=1/2\)

Fig. 13
figure 13

Multi-shells with one ansatz. \(L_0=48\), \(M_0=1.3\), \(k=0\) and \(l=3/2\)

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 Fig. 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 neighbouring 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 Fig. 11 we use the same parameters as in Fig. 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 Fig. 12. Finally we consider the case \(l=3/2\) but otherwise the same parameters as in Fig. 10. In this case three proper shells separated by vacuum are obtained as is depicted in Fig. 13.