# Effective four-dimensional loop quantum black hole with a cosmological
constant

Jianhui Lin  Department of Physics, South China University of Technology,
Guangzhou 510641, China  Xiangdong Zhang  Corresponding author.
<EMAIL_ADDRESS>Department of Physics, South China University of
Technology, Guangzhou 510641, China

###### Abstract

In this paper, we utilize the effective corrections of the $\bar{\mu}$-scheme
in loop quantum black holes to obtain a 4-dimensional spherically symmetric
metric with a cosmological constant. By imposing the areal gauge on the
components of Ashtekar variables in the classical theory and applying the
holonomy corrections, we derive the equations of motion, which can be solved
to obtain the expression for the effective metric in the Painlevé-Gullstrand
coordinates. A comparison with the $\Lambda=0$ case reveals minimal
modifications near the outer horizon, while significant differences are
observed far from the outer horizon. Moreover, the physical properties of
these quantum-corrected solutions are also discussed.

##  I Introduction

The prediction of the singularity at the center of black holes by classical
General Relativity (GR) has led to the widespread belief that classical theory
has limitations and that quantum gravity effects need to be introduced to cure
this type of spacetime singularity. Pursuing a consistent quantum gravity
theory becames one of the greatest challenges since the 20th century.

One approach to investigating how quantum gravity affects the spacetime of
black holes is to start with a specific quantum gravity theory and determine
its model corresponding to spherically symmetric spacetime, and then make
physical predictions based on the model. loop quantum gravity (LQG) is
currently one of the candidates for a theory of quantum gravity. As a
background independent and nonperturbative theory, it has considerable appeal
in this regard (see, e.g., [1, 2, 3, 4]). Since the late 1980s, LQG based on
Ashtekar variables has seen significant development up to the present day.
This includes natural predictions for discrete geometrical spectrum [5, 6, 7,
8, 9, 10] and successful generalizations to metric $f(R)$ theories, higher-
dimensional gravity and so on [11, 12, 24]. It preserves two fundamental
principles of general relativity: diffeomorphism invariance and background
independence. The construction of LQG adheres to mathematical rigor and
physical self-consistency. In situations where current experimental conditions
cannot be satisfied, it is necessary to theoretically consider whether the
classical limit of quantum theory is correct.

The application of LQG to cosmology, known as loop quantum cosmology (LQC),
has been established as an appropriate semiclassical state [14], and the
expectation values of the quantum operators in this state match well with
their corresponding classical values. It has also led to the conclusion that
the "Big Bounce" replaces the big bang, successfully resolving the problem of
the cosmological singularity of the big bang [15, 16, 17]. LQC is the
symmetry-reduced model of LQG [18]. In the classical scenario, before
quantization, we can utilize the symmetries of spatial homogeneity and
isotropy to reduce the phase space of gravitational degrees of freedom from
infinite dimensions to finite dimensions. Then, by using the methods and
techniques of LQG, it is able to proceed with its quantization. For a detailed
review of LQC, see, e.g., [19, 20].

Due to the success of LQC in resolving classical cosmological puzzles, the
attempt to apply this technique to black holes to address the issue of black
hole singularity is a very intriguing idea. As the simplest black hole
solution, the Schwarzschild black hole then serves as an ideal arena to
implement these ideas. Note that the interior of Schwarzschild black hole is
isometric to the Kantowski-Sachs model [17, 21], thus allowing for the
potential adaptation of techniques and ideas from LQC. Therefore, in the past
decade or so, the exploration of loop quantum black hole models has been a
popular direction.

Generally speaking, according to different quantum parameters, these loop
quantum black hole models can be classified into three main types: the
$\mu_{0}$-scheme [17], the $\bar{\mu}$-scheme [21, 22], and the so-called
generalized $\mu_{0}$-scheme [23, 24]. Among these, the $\bar{\mu}$-scheme
encounters an issue of excessive quantum corrections at the event horizon
[26]. In the classical regime, excessive quantum corrections are generally
considered unacceptable. However, recently, as suggested by some authors [27],
the reason of traditional $\bar{\mu}$-scheme cannot be implemented is simply
because we choose the wrong set of coordinates which becomes null at the
horizon. Hence they suggested to implement the $\bar{\mu}$-scheme in terms of
another set of coordinates which will not become null at the horizon. The use
of everywhere spacelike Painlevé-Gullstrand coordinates can be employed for
the $\bar{\mu}$-scheme and avoid the aforementioned problem, leading to an
effective framework for spherically symmetric vacuum solutions with holonomy
corrections from loop quantum gravity [27].

While for the generalized $\mu_{0}$-scheme which usually referred to as the
Ashtekar-Olmedo-Singh (AOS) approach [23]. In AOS scheme, the quantum
regularization parameters are set to be Dirac observables, i.e. being
constants along each dynamics trajectory but allowed to vary from one to
another. AOS model provides a nice effective description for the Kruskal
extension of both the exterior and interior of the Schwarzschild black holes
with large mass.

Up to now, most of the studies on loop quantum black hole models are limited
to the Schwarzschild case. Since the quantization scheme is not unique to
Schwarzschild black holes, generalizing the quantization scheme to other types
of black holes is crucial to test its universality. One of the extensions is
to include a cosmological constant and to study the quantization of
Schwarzschild de-Sitter (anti de-Sitter) case. The inclusion of a cosmological
constant in a black hole solution is important both in the practical and the
theoretical sense. However, generalizing the loop quantization scheme to
Schwarzschild de-Sitter(anti de-Sitter) case is not that easy. For the
generalized $\mu_{0}$-scheme, including a cosmological constant will add a
term $\Lambda\sqrt{q}$ with $\sqrt{q}$ being the determinant of spatial metric
in Hamiltonian constraint and this term does not commute with $p_{b}$ and
$p_{c}$ [23]. Hence it is very difficult to include a cosmological constant in
generalized $\mu_{0}$-scheme. Therefore, in this paper, we will follow the
line of [27] by utilizing $\bar{\mu}$ scheme in Painlevé-Gullstrand
coordinates to study the loop quantization of Schwarzschild de-Sitter(anti de-
Sitter) black holes.

This paper will be organized as follows: We will start from the classical
Hamiltonian framework and, in Section II, present the specific form of the
metric using the components of connection and densitized triad as dynamical
variables. In Section III, the LQG holonomy corrections are employed and
obtain corresponding metric solutions. Additionally, a discussion of the
physical properties of the solutions is provided in Section IV. Finally,
Section V concludes the article with a summary.

Our conventions are the following: space-time indices are denoted by
$\mu,\nu,\rho,\sigma,...$; spatial indices are denoted by $a;b;c;...$; and
internal indices are denoted by $i,j,k,...$ We only set $c=1$ but keep $G$ and
$\hbar$ explicit.

##  II CLASSICAL THEORY

The metric of a 4-dimensional spherically symmetric spacetime can be expressed
in Painlevé-Gullstrand coordinates as

$\displaystyle ds^{2}=-N^{2}dt^{2}+f^{2}(dx+N^{x}dt)^{2}+y^{2}d\Omega^{2},$
(1)

where the lapse $N$, shift vector $N^{x}$ and $f,y$ all are the function of
time $t$ and the radial coordinate $x$, while $d\Omega$, given by
$d\Omega^{2}=d\theta^{2}+\sin^{2}\theta d\phi^{2}$, is the line element on the
unit sphere.

As the foundation of LQG, connection dynamics is a theory based on the
Hamiltonian formulation of GR, described by Ashtekar-Barbero connections and
their conjugate momentum in terms of triads. In loop quantum black hole
models, we also inherit this feature and use them as basic variables. In this
context, the basic variables are consistent with [27], which respectively are
the Ashtekar-Barbero connection components given by $c(x)$, $p(x)$, and the
densitized triad components given by $E^{c}(x)$, $E^{p}(x)$. The following
represent the densitized triads in terms of metric components [27]

$\displaystyle E^{x}_{1}$ $\displaystyle=$ $\displaystyle
y^{2}sin\theta=E^{c}sin\theta,\qquad
E^{\theta}_{2}=fysin\theta=E^{p}sin\theta,$ $\displaystyle E^{\phi}_{3}$
$\displaystyle=$ $\displaystyle fy=E^{p}.$ (2)

Moreover, the Ashtekar-Barbero connection are written as
$A^{i}_{a}=\Gamma^{i}_{a}+\gamma K^{i}_{a}$, where both the spin-connection
$\Gamma^{i}_{a}$ and the extrinsic curvature $K^{i}_{a}$ only have three non-
zero components [27]

$\displaystyle K^{1}_{x}$ $\displaystyle=\frac{c}{\gamma},$
$\displaystyle\quad K^{2}_{\theta}$ $\displaystyle=\frac{p}{\gamma},$
$\displaystyle\quad K^{3}_{\phi}$ $\displaystyle=\frac{psin\theta}{\gamma}$
$\displaystyle\Gamma^{3}_{\theta}$
$\displaystyle=-\frac{\partial_{x}E^{c}}{2E^{p}},$
$\displaystyle\quad\Gamma^{1}_{\phi}$ $\displaystyle=-cos\theta,$
$\displaystyle\quad\Gamma^{2}_{\phi}$
$\displaystyle=\frac{\partial_{x}E^{c}sin\theta}{2E^{p}},$ (3)

where $\gamma$ is the Barbero-Immirzi parameter. The determinant of the
spatial metric is $q=E^{c}(E^{p}sin\theta)^{2}$.

###  II.1 Actions and constraints

The action of GR plus a cosmological constant reads

$\displaystyle S=\frac{1}{16\pi G}\int
d^{4}x\sqrt{-g}\left(R-2\Lambda\right).$ (4)

After performing canonical analysis, the gravitational action becames

$\displaystyle S=\int
dt\int_{\Sigma}\left[\frac{\dot{A}^{i}_{a}E^{a}_{i}}{8\pi
G\gamma}-N\mathcal{H}_{\Lambda}-N^{a}\mathcal{H}_{a}\right],$ (5)

where the dot means the derivative with respect to time $t$; the Hamiltonian
constraint $\mathcal{H}_{\Lambda}$ and the diffeomorphism constraintis reads
respectively as

$\displaystyle\mathcal{H}_{\Lambda}$ $\displaystyle=$
$\displaystyle-\frac{E^{a}_{i}E^{b}_{j}}{16\pi
G\gamma^{2}\sqrt{q}}{\epsilon^{ij}}_{k}({F_{ab}}^{k}-(1+\gamma^{2}){\Omega_{ab}}^{k}))$
(6) $\displaystyle-\frac{\sqrt{q}N}{8\pi G}\Lambda.$
$\displaystyle\mathcal{H}_{a}$ $\displaystyle=$ $\displaystyle\frac{1}{8\pi
G\gamma}E^{b}_{k}{F_{ab}}^{k}.$ (7)

Here the field strength is given by
$F_{ab}{}^{k}=2\partial_{[a}A_{b]}^{k}+\epsilon_{ij}{}^{k}A_{a}^{i}A_{b}^{j}$
and the spatial curvature denotes
$\Omega_{ab}{}^{k}=2\partial_{[a}\Gamma_{b]}^{k}+\epsilon_{ij}{}^{k}\Gamma_{a}^{i}\Gamma_{b}^{j}$.
Due to the spherical symmetry, we can integrate over $d\Omega$ and use the Eq.
(II), then obtain the following expression

$\displaystyle S=\int dt\int
dx{\left[\frac{\dot{c}E^{c}+2\dot{p}E^{p}}{2G\gamma}-N\mathcal{H}_{\Lambda}-N^{x}\mathcal{H}_{x}\right]},$
(8)

where the Hamiltonian constraint reads

$\displaystyle\mathcal{H}_{\Lambda}=$ $\displaystyle-$
$\displaystyle\frac{1}{2G\gamma}\left[\frac{2cp\sqrt{E^{c}}}{\gamma}+\frac{E^{p}}{\gamma\sqrt{E^{c}}}(p^{2}+\gamma^{2})-\frac{\gamma(\partial_{x}E^{c})^{2}}{4E^{p}\sqrt{E^{c}}}\right.$
(9) $\displaystyle-$
$\displaystyle\left.\gamma\sqrt{E^{c}}\partial_{x}\left(\frac{\partial_{x}E^{c}}{E^{p}}\right)-\gamma\sqrt{E^{c}}E^{p}\Lambda\right],$

and the diffeomorphism constraint is

$\displaystyle\mathcal{H}_{x}=\frac{1}{2G\gamma}(2E^{p}\partial_{x}p-c\partial_{x}E^{c}).$
(10)

Here straightforward calculations show only $\mathcal{H}_{x}$ is non-zero and
$\mathcal{H}_{\theta},\mathcal{H}_{\phi}$ both vanish. Moreover, from Eq. (8),
the symplectic structure of the symmetry-reduced black hole models is given by

$\displaystyle\\{c(x_{1}),E^{c}(x_{2})\\}=2G\gamma\delta(x_{1}-x_{2}),$ (11)
$\displaystyle\\{p(x_{1}),E^{p}(x_{2})\\}=G\gamma\delta(x_{1}-x_{2}),$ (12)

with $\delta$ being the Dirac delta function. We now consider the smeared
constraint function $\mathcal{C}[N]=\int dxN\mathcal{H}_{\Lambda}$ and
$\mathcal{D}[N^{x}]=\int dxN^{x}\mathcal{H}_{x}$. Through lengthy but
straightforward calculations, the Poission brackets between constraints
algebra read

$\displaystyle\\{\mathcal{C}[N_{1}],\mathcal{C}[N_{2}]\\}=\mathcal{D}\left[\frac{E^{c}}{(E^{p})^{2}}(N_{1}\partial_{x}N_{2}-N_{2}\partial_{x}N_{1})\right],$
(13)
$\displaystyle\\{\mathcal{D}[N_{1}^{x}],\mathcal{D}[N_{2}^{x}]\\}=\mathcal{D}[(N_{2}^{x}\partial_{x}N_{1}^{x}-N_{1}^{x}\partial_{x}N_{2}^{x})],$
(14)
$\displaystyle\\{\mathcal{C}[N],\mathcal{D}[N^{x}]\\}=-\mathcal{C}[N^{x}\partial_{x}N].$
(15)

As compared to the case without the cosmological constant [27], the symplectic
structure and the constraint algebra have not changed. Using the Hamilton
evolution equations
$\dot{y}=\left\\{y,\mathcal{C}[N]+\mathcal{D}[N^{x}]\right\\}$, the equations
of motion for each component with respect to time are determined by

$\displaystyle\dot{E}^{c}=\frac{2Np}{\gamma}\sqrt{E^{c}}+N^{x}\partial_{x}E^{c},$
(16)
$\displaystyle\dot{E}^{p}=\frac{N}{\gamma\sqrt{E^{c}}}(cE^{c}+pE^{p})+\partial_{x}(N^{x}E^{p}),$
(17) $\displaystyle\dot{c}=$
$\displaystyle\frac{N}{2\gamma\sqrt{E^{c}}}\left[\frac{E^{p}}{E^{c}}(p^{2}+\gamma^{2})-2cp\right]+\gamma\partial_{x}\biggl{(}\frac{\partial_{x}(N\sqrt{E^{c}})}{E^{p}}\biggr{)}$
$\displaystyle+\frac{N\gamma}{2\sqrt{E^{c}}}\left[\partial_{x}\biggl{(}\frac{\partial_{x}E^{c}}{E^{p}}\biggr{)}-\frac{(\partial_{x}E^{c})^{2}}{4E^{c}E^{p}}\right]+\partial_{x}(N^{x}c)$
$\displaystyle-\frac{\gamma}{2}\partial_{x}{\left(N\frac{\partial_{x}E^{c}}{E^{p}\sqrt{E^{c}}}\right)}+\frac{N\gamma
E^{p}\Lambda}{2\sqrt{E^{c}}}$ (18) $\displaystyle\dot{p}=$ $\displaystyle-$
$\displaystyle\frac{N}{2\gamma\sqrt{E^{c}}}(p^{2}+\gamma^{2})-\frac{\gamma}{2}\left[\frac{N}{\sqrt{E^{c}}}\left(\frac{\partial_{x}E^{c}}{2E^{p}}\right)^{2}\right.$
(19) $\displaystyle-$
$\displaystyle\left.\partial_{x}(N\sqrt{E^{c}})\frac{\partial_{x}E^{c}}{(E^{p})^{2}}+N^{x}\partial_{x}c\right]+\frac{N\gamma\sqrt{E^{c}}\Lambda}{2}$

Once the lapse and shift vector are chosen, by solving the above equations of
motion and the diffeomorphism constraint together with the Hamiltonian
constraint, one can obtain the vacuum spherically symmetric metric with the
cosmological constant.

###  II.2 Classical solutions

Following the same line that in [27], we also adapt the areal gauge for
$E^{c}$ as follows:

$\displaystyle E^{c}=x^{2}.$ (20)

Once the areal gauge is imposed, the diffeomorphism constraint (10) then in
turn implies

$\displaystyle\quad c=\frac{E^{p}}{x}\partial_{x}p.$ (21)

The gauge-fixing condition $\zeta=E^{c}-x^{2}$ is second class with the
diffeomorphism constraint $\mathcal{H}_{x}$, which should be preserved by the
equations of motion. That is to say, $\dot{\zeta}=0$. Using the condition
$\dot{E}^{c}=0$, we can obtain the expression

$\displaystyle N^{x}=-\frac{Np}{\gamma},$ (22)

which implies that $N^{x}$ is no longer a Lagrange multiplier. After choosing
the areal gauge, $N^{x}$ is determined by the lapse function $N$. Through this
process, the action can be simplified and given by:

$\displaystyle S_{GF}=\int dt\int
dx\left[\frac{\dot{p}E^{p}}{G\gamma}-N\mathcal{H}_{\Lambda}\right]$ (23)

with

$\displaystyle\mathcal{H}_{\Lambda}=$
$\displaystyle\frac{1}{2G\gamma}\left[\frac{3\gamma x}{E^{p}}-\frac{2\gamma
x^{2}}{(E^{p})^{2}}\partial_{x}E^{p}-\frac{E^{p}}{\gamma
x}\partial_{x}[x(p^{2}+\gamma^{2})]\right]$
$\displaystyle+\frac{xE^{p}}{2G}\Lambda.$ (24)

Since the $E^{c}$ is fixed by the areal gauge, the remaining symplectic
structure of the symmetry-reduced theory reads

$\displaystyle\\{p(x_{1}),E^{p}(x_{2})\\}=G\gamma\delta(x_{1}-x_{2}),$ (25)

and similarly, the remaining constraint algebra is

$\displaystyle\\{\mathcal{C}[N_{1}],\mathcal{C}[N_{2}]\\}$ $\displaystyle=$
$\displaystyle
C{\left[-\frac{1}{\gamma}\left(N_{1}\partial_{x}N_{2}-N_{2}\partial_{x}N_{1}\right)p\right]}$
(26) $\displaystyle=$ $\displaystyle
C\left[N_{1}^{x}\partial_{x}N_{2}-N_{2}^{x}\partial_{x}N_{1}\right],$

where the second step utilizes Eq. (22), which corresponds to the Painlevé-
Gullstrand coordinates. The equations of motion can also be simplified to

$\displaystyle\dot{E}^{p}=$ $\displaystyle\frac{p}{\gamma
x}\left(NE^{p}-x\partial_{x}\left(NE^{p}\right)\right),$ (27)
$\displaystyle\dot{p}=$ $\displaystyle\frac{\gamma
Nx}{2(E^{p})^{2}}+\frac{\gamma
x^{2}}{(E^{p})^{2}}\partial_{x}N-\frac{N}{2x\gamma}\partial_{x}\left(xp^{2}+\gamma^{2}x\right)$
$\displaystyle+\frac{\gamma Nx}{2}\Lambda,$ (28)

which can be obtained by substituting areal gauge 20 and the expression (22)
into the equations of motion (17) and (19). Now we set $N=1$, and take into
account the Hamiltonian constraint $\mathcal{H}_{\Lambda}=0$. Then,
$\dot{E}^{p}=0$ and $\dot{p}=0$ can be solved to obtain

$\displaystyle E^{p}$ $\displaystyle=$ $\displaystyle C_{1}x,$ (29)
$\displaystyle p^{2}$ $\displaystyle=$
$\displaystyle\gamma^{2}\left(\frac{1}{{C_{1}^{2}}}-1\right)+\frac{\gamma^{2}x^{2}\Lambda}{3}+\frac{C_{2}}{x},$
(30)

where $C_{1}$ and $C_{2}$ are two constants of integration. Furthermore, these
solutions (29) and (30) satisfy the Hamiltonian constraint equation
$\mathcal{H}_{\Lambda}=0$ identically. From equation (II), it is worthy noting
that both $f$ and $y$ can be expressed in terms of $E^{c}$ and $E^{p}$ with

$\displaystyle f^{2}$ $\displaystyle=$
$\displaystyle\frac{(E^{p})^{2}}{E^{c}},$ (31) $\displaystyle y^{2}$
$\displaystyle=$ $\displaystyle E^{c}.$ (32)

Using the equations (22), (31), (32) and $N=1$, by inserting the solutions
(29) and (30) into the metric (1), we obtain

$\displaystyle ds^{2}=$ $\displaystyle-$
$\displaystyle\left(1-\frac{2GM}{x}-\frac{x^{2}\Lambda}{3}\right)dt^{2}$ (33)
$\displaystyle+$ $\displaystyle
2\sqrt{\frac{2GM}{x}+\frac{x^{2}\Lambda}{3}}dxdt+dx^{2}+x^{2}d\Omega^{2},$

where we can determine the integration constants $C_{1}=1$ and
$C_{2}=2GM\gamma^{2}$ by compare this solution with the standard Schwarzschild
de-Sitter(anti-de Sitter) solutions and $M$ is being the mass of the black
holes. Moreover, we can also confirm that the square root of $p$ is negative.
Finally, the solutions to the equations of motion are given by

$\displaystyle E^{p}$ $\displaystyle=$ $\displaystyle x,$ (34) $\displaystyle
p$ $\displaystyle=$
$\displaystyle-\sqrt{\frac{2GM\gamma^{2}}{x}+\frac{\gamma^{2}x^{2}\Lambda}{3}}.$
(35)

##  III Effective correction of LQG

We adopt the standard LQC treatment for Schwarzschild black holes which
requires to introduce holonomy corrections. The holonomy of an
$SU(2)$-connection $A_{a}^{i}$ is path-ordering exponential integral along an
edge $e^{a}$

$\displaystyle h(A)=\mathcal{P}exp\int_{e_{i}}dtA_{a}^{j}\tau_{j}e^{a},$ (36)

with $\mathcal{P}$ labeling path-ordering. Although it is not yet fully clear
whether the LQG dynamics of a loop quantum Schwarzschild black hole can also
provide a good approximation for the full LQG dynamics of the Schwarzschild
black hole spacetime, the success of LQC has led people to believe that the
effective dynamics of the loop quantum Schwarzschild black hole model can
serve as a good approximation for the semiclassical quantum dynamics of
observable objects with physical length scales much larger than the Planck
scale. As the basic operators in LQG are the holonomies and triads, it is
necessary to use the holonomies instead of the Ashtekar-Barbero components in
the present paper. In the resulting quantum effective Hamiltonian constraints,
the holonomy correction is simplified by replacing the components of Ashtekar
connection $p$ with [27]

$\displaystyle
p\rightarrow\frac{x}{\sqrt{\Delta}}\sin\left(\frac{\sqrt{\Delta}}{x}p\right)$
(37)

after the areal gauge is imposed. The effective Hamiltonian can be constructed
by replacing all the instances of $p$ in (II.2) with (37) and the result is

$\displaystyle\mathcal{H}^{(\mathrm{LQG})}_{\Lambda}=$ $\displaystyle-$
$\displaystyle\frac{1}{2G\gamma}\left[\frac{E^{p}}{\gamma
x}\partial_{x}\left(\frac{x^{3}}{\Delta}\sin^{2}\frac{\sqrt{\Delta}p}{x}+\gamma^{2}x\right)\right.$
(38) $\displaystyle-$ $\displaystyle\left.\frac{3\gamma
x}{E^{p}}+\frac{2\gamma x^{2}}{(E^{p})^{2}}\partial_{x}E^{p}-x\gamma
E^{p}\Lambda\right].$

And we can obtain the Poisson brackets between $\mathcal{C}^{\mathrm{(LQG)}}$

$\displaystyle\\{\mathcal{C}^{\mathrm{(LQG)}}[N_{1}],\mathcal{C}^{\mathrm{(LQG)}}[N_{2}]\\}$
$\displaystyle=C^{(\mathrm{LQG})}\left[-\frac{x}{\gamma\sqrt{\Delta}}\mathrm{sin}\frac{\sqrt{\Delta}p}{x}\mathrm{cos}\frac{\sqrt{\Delta}p}{x}\bigl{(}N_{1}\partial_{x}N_{2}\right.$
$\displaystyle-\left.N_{2}\partial_{x}N_{1}\bigr{)}\right].$ (39)

The equation (22) also needs to be changed, and at the same time, in order to
give desired constraint algebra (III), The equation (22) has been modified to
read as

$\displaystyle
N^{x}=-\frac{Nx}{\gamma\sqrt{\Delta}}\sin\frac{\sqrt{\Delta}p}{x}{\cos\frac{\sqrt{\Delta}p}{x}}.$
(40)

Finally just like the classical case, once the effective Hamiltonian 38 is in
hand, the modified equations of motion can be obtained

$\displaystyle\dot{E}^{p}=-\frac{x^{2}}{2\gamma\sqrt{\Delta}}\partial_{x}{\left(\frac{NE^{p}}{x}\right)}\sin\frac{\sqrt{\Delta}p}{x}{\cos\frac{\sqrt{\Delta}p}{x}},$
(41) $\displaystyle\dot{p}$ $\displaystyle=$ $\displaystyle\frac{\gamma
Nx}{2(E^{p})^{2}}\left(1+2x\frac{\partial_{x}N}{N}\right)-\frac{\gamma N}{2x}$
(42) $\displaystyle-$ $\displaystyle\frac{N}{2\gamma\Delta
x}\partial_{x}\bigg{(}x^{3}\sin^{2}\frac{\sqrt{\Delta}p}{x}\bigg{)}+\frac{x\gamma
N\Lambda}{2}.$

We now turn to solve these equations of motion in the Painlevé-Gullstrand
coordinates and set $N=1$. Taking into account the Hamiltonian constraint
$\mathcal{H}^{(\mathrm{LQG})}_{\Lambda}$, the complete set of equations is as
follows

$\displaystyle\dot{E}^{p}=0;\quad\dot{p}=0$ $\displaystyle
N=1;\quad\mathcal{H}^{(\mathrm{LQG})}_{\Lambda}=0$ (43)

Then by relating $\dot{E}^{p}=0$ and $N=1$, Eq. (41) implies three
possibilities:

$\displaystyle 1,\quad\partial_{x}{\left(\frac{E^{p}}{x}\right)}=0,$ (44)
$\displaystyle 2,\quad\sin\frac{\sqrt{\Delta}p}{x}=0,$ (45) $\displaystyle
3,\quad\cos\frac{\sqrt{\Delta}p}{x}=0.$ (46)

For the case $1$ we assuming the term
$\partial_{x}{\left(\frac{E^{p}}{x}\right)}$ in equation (41) vanishes, the
solutions is

$\displaystyle E^{p}=C_{3}x,$ (47)

then, by using Eq. (III), we have

$\displaystyle
p=\frac{x}{\sqrt{\Delta}}\arcsin\left[\left(\frac{1}{C^{2}_{3}-1}\right)\gamma^{2}\Delta\frac{1}{x^{2}}+\frac{C_{4}}{x^{3}}\right.$
$\displaystyle+\left.\frac{\gamma^{2}\Delta\Lambda}{3}\right]^{1/2},$ (48)

where $C_{3}$ and $C_{4}$ are constants of integration. It is not hard to
confirm that the solutions (47) and (III) satisfy the Hamiltonian constraint
$\mathcal{H}^{(\mathrm{LQG})}_{\Lambda}=0$. For case $2$ and case $3$, use the
same logic of [27] and after careful calculations, it is determined that
neither of them satisfies the Hamiltonian constraint. Therefore, the solution
has only one form.

Next, we need to determine the values of these integration constants $C_{3}$
and $C_{4}$. Noticing equations (31), (32), (47) and $N=1$, the metric (1) can
now be written in the following form

$\displaystyle ds^{2}=$ $\displaystyle-$
$\displaystyle\left(1-\left(C_{3}N^{x}\right)^{2}\right)dt^{2}+2C_{3}^{2}N^{x}dxdt$
(49) $\displaystyle+$ $\displaystyle C_{3}^{2}dx^{2}+x^{2}d\Omega^{2},$

where we should note that the LQG correction terms are present in $N^{x}$.
Because we want the target metric to recover the classical case (33) when the
LQG corrections are removed, it is easy to know that $C_{3}^{2}=1$ in this
case. Substituting $p$ (III) into $N^{x}$ (40), we obtain the desired result

$\displaystyle
N^{x}=\left(\frac{x^{2}\Lambda}{3}+\frac{C_{4}}{\gamma^{2}\Delta
x}\right)^{1/2}\left(1-\frac{\gamma^{2}\Delta\Lambda}{3}-\frac{C_{4}}{x^{3}}\right)^{1/2},$
(50)

taking into account that

$\cos\frac{\sqrt{\Delta}p}{x}=-\sqrt{1-\sin^{2}\frac{\sqrt{\Delta}p}{x}}$

for consistency with the classical case. When $\Delta\rightarrow 0$, $N^{x}$
should revert to the classical case (which can be given by using equations
(22) and (35)). Hence, we can ascertain that $C_{4}=2GM\gamma^{2}\Delta$.
Then, the resulting LQG-corrected metric with the cosmological constant reads

$\displaystyle ds^{2}=-F(x)dt^{2}+2N^{x}dxdt+dx^{2}+x^{2}d\Omega^{2},$ (51)

with

$\displaystyle F(x)=1-\left(N^{x}\right)^{2},$ (52) $\displaystyle
N^{x}=\left(\frac{x^{2}\Lambda}{3}+\frac{2GM}{x}\right)^{1/2}\left(1-\frac{2GM\gamma^{2}\Delta}{x^{3}}\right.$
$\displaystyle-\left.\frac{\gamma^{2}\Delta\Lambda}{3}\right)^{1/2}.$ (53)

From the above result, it can be easily verified that as $\Delta\rightarrow
0$, the classical Schwarzschild de-Sitter(anti-de Sitter) solutions are
recovered. Another point worth mentioning is that when $\Lambda$ is set to
zero, the metric (51) can be restored to the case without cosmological
constant as that in [27], which further demonstrates the correctness of the
results of the present paper.

##  IV Physical properties of the effective metric

###  IV.1 Geodesic and Effective Range

As stated in [25] [27], the metric corrected by LQG is only valid for
$x>x_{min}$. This is because, in the LQC, there is an upper limit to the
energy density of any matter field. Therefore, in order to generate a
gravitational field with mass $M$, the corresponding matter field cannot be
infinitely small but will have a radius corresponding to the maximum energy
density. Now we will search for the lower bound $x_{min}$.

Investigating geodesics provides a direct method for understanding the
structure of spacetime. To simplify the analysis, we will focus solely on
radial geodesics that satisfy the following equation

$\displaystyle-\kappa=-F(x)\left(\frac{dt}{d\tau}\right)^{2}+2N^{x}\left(\frac{dx}{d\tau}\right)\left(\frac{dt}{d\tau}\right)+\left(\frac{dx}{d\tau}\right)^{2}$
(54)

where $\tau$ represents the geodesic parameter. When $\kappa=1$, it
corresponds to timelike geodesics, and $\tau$ represents proper time. On the
other hand, if $\kappa=0$, it refers to null geodesics, and $\tau$ represents
a chosen affine parameter. It is important to note that in this case of a
stationary spherically symmetric spacetime, $T^{\mu}=(\partial/\partial
t)^{\mu}$ is a Killing vector. Setting the tangent vector of the geodesic as
$\xi^{\mu}=(\partial/\partial\tau)^{\mu}$, then $T^{\mu}\xi_{\mu}$ is constant
along the geodesic, which we refer to as the conserved energy written as

$\displaystyle
E=T^{\mu}\xi_{\mu}=-F(x)\frac{dt}{d\tau}+N^{x}\frac{dx}{d\tau}.$ (55)

For timelike geodesics $(\kappa=1)$ and considering $E=1$, it corresponds to a
stationary particle starting its motion from infinity. The geodesic equation
(54) can be further simplified as follows

$\displaystyle\frac{dx}{d\tau}=\pm\sqrt{1-F(x)},$ (56)

where it can be observed that $x_{min}$ is at $dx/d\tau=0$ and the solution is

$\displaystyle
x_{min}=\left(\frac{6GM\gamma^{2}\Delta}{3-\gamma^{2}\Delta\Lambda}\right)^{1/3}.$
(57)

Compared with the Schwarzschild case [27], it is easy to see that the term
involving the cosmological constant causes a slight variation on $x_{min}$.
However, when this term is neglected, we obtain exactly the same value of
$x_{min}$ as discussed in [27, 28]. Moreover, since the $x_{min}$ should be
greater than zero, then the denominator of Eq. (57) implies an upper bound on
the cosmological constant as $\Lambda<\frac{3}{\gamma^{2}\Delta}$.

Besides, for null geodesic, the situation will be much simpler. Since
$\kappa=0$, equation 54 can be divided by $dt/d\tau$ and written as

$\displaystyle\frac{dx}{dt}=-N^{x}\pm 1.$ (58)

A nontrivial case arises when the outgoing null rays with $dx/dt=-N^{x}+1$
have zeros at $x>0$. However, the ingoing rays with $dx/dt=-N^{x}-1$ are
always less than zero. One point worth mentioning here is that it was
previously stated that $T^{a}$ represents a Killing vector field, and its
corresponding Killing horizon is given by $T^{a}T_{a}=0$, which can be
simplified as $F(x)=0$. But in fact, $F(x)=0$ and $dx/dt=0$ are equivalent to
each other, since the zeros of $dx/dt$ also represent the Killing horizon.
According to [27], $M_{*}=8\gamma\sqrt{\Delta}/\sqrt{27}G$ corresponding to
the critical mass means that when $\Lambda=0$ and, there is only one Killing
horizon. More precisely, as the mass is greater than $M_{*}$, there exist two
Killing horizons, while if the mass is smaller than $M_{*}$, there are no
Killing horizons. The detailed situation can be visualized in the figure 1.

Figure 1:  For outgoing null geodesics, we have plotted the case of the black
hole with mass $M_{*}$ and compared three different spacetimes. AdS spacetime
is taken to have $\Lambda=-0.1$, while dS spacetime has $\Lambda=0.1$. It can
be observed that the three spacetimes exhibit little difference near the
Killing horizon, and only significant differences arise far from the horizon.
The remaining parameters are taken as
$\left\\{\gamma=1,\Delta=0.01,G=1\right\\}$.

Figure 2:  For outgoing null geodesics, we have separately plotted the cases
of the AdS black holes (first figure) and the dS black holes (second figure)
with the mass of $3M_{*}$. In both figures, the red curve represents the LQG-
corrected metric, while the blue curve means the classical metric. For AdS
black holes, we consider $\Lambda=-0.1$, whereas for dS black holes,
$\Lambda=0.1$. All the remaining parameters are taken as
$\left\\{\gamma=1,G=1\right\\}$.

From Fig 2, in the case of the LQG-corrected AdS black hole, the mass exceeds
$M_{*}$, resulting in the presence of both inner and outer horizons. For LQG-
corrected metric, the endpoints of the red curve correspond to the minimum and
maximum values of $x$, where $dx/dt=1$, as can be easily verified by equations
(52), (56), and (58). Conversely, the classical AdS black hole exhibits only
one horizon, and $dx/dt$ diverges rapidly as $x\rightarrow 0$, indicating the
presence of the spacetime singularity. In the case of dS black hole, the LQG
correction also results in the existence of two Killing horizons, avoiding the
singularity issue. Additionally, the cosmological horizon can be observed from
the curves.

###  IV.2 Curvature scalars and surface gravity

Investigating the curvature of this effective metric can provide us with a
deeper understanding of the geometric structure of this spacetime.
Furthermore, it allows for a more straightforward comparison to be made with
the scenario where the cosmological constant is not present. With the
corrected metric 51, various curvature scalars can be straightforwardly
calculated. The results are summarized as follows:

$\displaystyle
R=4\Lambda-\frac{4}{3}\gamma^{2}\Delta\Lambda^{2}-\frac{24(GM\gamma)^{2}\Delta}{x^{6}},$
(59) $\displaystyle
R_{\mu\nu}R^{\mu\nu}=\frac{4}{9}\Lambda^{2}\left(-3+\gamma^{2}\Delta\Lambda\right)^{2}+\frac{1440(GM\gamma)^{4}\Delta^{2}}{x^{12}}$
$\displaystyle+\frac{16(GM\gamma)^{2}\left(-3+\gamma^{2}\right)\Delta\Lambda}{x^{6}},$
(60) $\displaystyle R_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma}$ $\displaystyle=$
$\displaystyle\frac{8}{27}\Lambda^{2}\left(-3+\gamma^{2}\Delta\Lambda\right)^{2}$
(61) $\displaystyle+$
$\displaystyle\frac{320(GM)^{3}\gamma^{2}\left(-3+2\gamma^{2}\Delta\Lambda\right)\Delta}{x^{9}}$
$\displaystyle+$
$\displaystyle\frac{32\left(GM\gamma\right)^{2}\left(-3+\gamma^{2}\Delta\Lambda\right)\Delta\Lambda}{x^{6}}$
$\displaystyle+$ $\displaystyle\frac{7488(GM\gamma)^{4}\Delta^{2}}{{x}^{12}}.$

For physical reasons, the absolute value of curvature scalars should
correspond to a maximum when $x$ approaches $x_{min}$ in vacuum. By
substituting $x_{min}$ into various curvature scalars, we obtain the following
results

$\displaystyle\lim_{x\to
x_{\min}}R=-\frac{6}{\gamma^{2}\Delta}+8\Lambda-2\gamma^{2}\Delta\Lambda^{2},$
(62) $\displaystyle\lim_{x\to x_{\min}}R_{\mu\nu}R^{\mu\nu}$
$\displaystyle=\frac{2\left(-3+\gamma^{2}\Delta\Lambda\right)^{2}\left(5-4\gamma^{2}\Delta\Lambda+\gamma^{4}\Delta^{2}\Lambda^{2}\right)}{\gamma^{4}\Delta^{2}},$
(63) $\displaystyle\lim_{x\to
x_{\min}}R_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma}$
$\displaystyle=\frac{4\left(-3+\gamma^{2}\Delta\Lambda\right)^{2}\left(10-6\gamma^{2}\Delta\Lambda+\gamma^{4}\Delta^{2}\Lambda^{2}\right)}{\gamma^{4}\Delta^{2}}.$
(64)

These are also constants independent of the black hole mass. Unsurprisingly,
when the terms containing the cosmological constant are discarded, these
curvature scalars align with those in [27] [28].

In addition, exploring the trapped surface can further verify the results of
the effective metric we have obtained. Considering the null geodesic
congruence, due to the spherical symmetry, on a 2-dimensional sphere $S$, the
radial coordinate $x$ and the time coordinate $t$ are constant. Therefore,
there are two independent null tangent vectors orthogonal to $S$, referred to
as outgoing null geodesics and ingoing null geodesics. These two geodesics
correspond to the outgoing expansion parameter $\theta_{+}$ and ingoing
expansion parameter $\theta_{-}$ respectively. By definition, the trapped
surface corresponds to the region where both $\theta_{+}$ and $\theta_{-}$ are
negative. Here, we define the tangent vectors of the outgoing null geodesics
and the ingoing null geodesics to be the same as [27], denoted by

$\displaystyle l^{\mu}=\left(\frac{\partial}{\partial
t}\right)^{\mu}+(1-N^{x})\left(\frac{\partial}{\partial x}\right)^{\mu},$ (65)

and

$\displaystyle k^{\mu}=\left(\frac{\partial}{\partial
t}\right)^{\mu}+(-1-N^{x})\left(\frac{\partial}{\partial x}\right)^{\mu}.$
(66)

Normalize these two vector fields such that $l^{\mu}k_{\mu}=-2$ and the
hypersurface metric for $S$ is

$\displaystyle
h_{\mu\nu}=g_{\mu\nu}+\frac{1}{2}\left(l_{\mu}k\nu+k_{\mu}l_{\nu}\right).$
(67)

Based on the definition of the expansions $\theta_{+}$ and $\theta_{-}$, the
simple expression can be obtained through straightforward calculations and the
results are

$\displaystyle\theta_{+}$ $\displaystyle=$ $\displaystyle
h^{\mu\nu}\nabla_{\mu}l_{\nu}=\frac{2}{x}\left(1-N^{x}\right),$ (68)
$\displaystyle\theta_{-}$ $\displaystyle=$ $\displaystyle
h^{\mu\nu}\nabla_{\mu}k_{\nu}=-\frac{2}{x}\left(1+N^{x}\right).$ (69)

Here, $\theta_{-}$ remains negative within the range of $x$, so considering
$\theta_{+}$ alone is sufficient. The main results are displayed in Figure 4
and 4. Besides, the trapped surface tends to increase with the increase of
mass. From the expression of the expansion (68) and the null geodesic (58), it
can be observed that they differ only by a factor of $2/x$. Therefore, it is
conceivable that the zero points of the null geodesics align with the outgoing
expansion and the physical interpretation is consistent with the null
geodesics.

Figure 3: The expansions of the outgoing null geodesics are compared for
different cosmological constants. In this case, the mass is chosen as $M_{*}$,
indicating that when $\Lambda=0$ , there is only one event horizon. More
detailed calculations reveal when $\Lambda=-0.1$, there is no event horizon
and the curve’s endpoint cannot extend indefinitely but reaches a certain
point. This also implies the maximum value that $x$ can take. When
$\Lambda=0.1$, there are two close inner and outer horizons, and a
cosmological horizon can be clearly observed at the end of the curve. These
properties are analogous to those of the outgoing null geodesics. The
remaining parameters are taken as $\left\\{\gamma=1,\Delta=0.01,G=1\right\\}$.

Figure 4:  Comparison of the null expansion $\theta_{+}$ for the LQG-corrected
metric with different values of the cosmological constant. The mass is
$2M_{*}$, compared to the case with mass $M_{*}$, the most notable difference
is that regardless of the value of the cosmological constant, the black hole
exhibits both inner and outer event horizons. The influence of the
cosmological constant is only evident at the endpoint of the curve. The
remaining parameters are taken as $\left\\{\gamma=1,\Delta=0.01,G=1\right\\}$.

When a black hole possesses an outer horizon, it becomes straightforward to
calculate the surface gravity by employing

$\displaystyle\kappa=\left.\frac{1}{2}\frac{\partial F(x)}{\partial
x}\right|_{x=x_{out}}.$ (70)

where $x_{out}$ represents the position of the outer horizon. Due to the
complexity of the expression of $F(x)$, it is not practical to calculate the
analytical solution of the surface gravity directly. Moreover, if we adopt a
first-order approximation for $\Delta$, the distinctions between the three
types of spacetime will not be apparent and the result of the expression will
be the same. Therefore, we adopt a numerical approach to solve the problem.
The detailed results are shown in Figure 5.

Figure 5:  By taking the mass of the black hole as the independent variable,
we compared the surface gravity under different values of the cosmological
constant. The mass is measured in units of $M_{*}$. The minimum value of the
mass is set at $M_{*}$ because when the mass is less than $M_{*}$, there is no
horizon, and any discussion concerning surface gravity becomes meaningless.
The remaining parameters are taken as
$\left\\{\gamma=1,\Delta=0.01,G=1\right\\}$.

The green curve corresponds to the AdS black hole, and the red curve
represents the dS black hole. At smaller masses, due to the LQG correction,
they behave almost identically and both exhibit a peak. However, as the mass
increases, due to the distinct characteristics of the spacetime, the surface
gravity of the dS black hole rapidly tends to zero and eventually reaches an
upper mass limit, whereas the surface gravity of the AdS black hole first
decreases and then slowly increases.

##  V Conclusions and discussions

The main result in this paper is obtained based on the Painlevé-Gullstrand
coordinate system and utilizes the components of the Ashtekar-Barbero
connection and density-triad to describe a 4-dimensional spherically symmetric
spacetime with a cosmological constant. We first perform calculations in the
classical scenario. Starting from the gravitational action, after symmetry
reduction, we consider the constrained system to consist of Hamiltonian
constraint (scalar constraint) and the diffeomorphism constraint (vector
constraint) due to the constrained nature of the system. By exploiting the
spherical symmetry, further simplification of the constraints can be achieved.
The equations of motion can then be obtained.

Then we impose a gauge fixing on the area and select a particular Painlevé-
Gullstrand coordinate system $(N=1)$ [27]. As a result, the number of
dynamical variables is reduced from four to two, and the equations of motion
are also reduced to two. By directly solving the equations, the classical
Schwarzschild de-Sitter(anti-de Sitter) solutions are obtained as expected.

We then modified the variable $p$ of the connection by incorporating holonomy
corrections and adjusted the Hamiltonian constraint, symplectic structure, and
shift vector accordingly. Similarly, by solving the modified effective
equations of motion, we obtained the LQG-corrected spherically symmetric
metric with a cosmological constant.

Due to the LQG correction, the variable $x$ has a restricted range of values.
For AdS and dS spacetimes, we calculated their common minimum value $x_{min}$
using timelike geodesics. The dS spacetime exhibits the expected cosmological
horizon. The effect of the cosmological constant is only noticeable at
distances far from the event horizon which is consistent with the
expectations. We also investigated the curvature scalars
$R,R_{\mu\nu}R^{\mu\nu},R_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma}$ of the LQG-
corrected metric, which are bounded at $x_{min}$ due to quantum gravity
effects and are constants that only depend on $\Lambda$ and $\Delta$.
Moreover, when setting $\Lambda=0$, all the curvature scalars can be recovered
in the form presented in article [27].

It is worth mentioning that the computation of $x_{min}$ can also be directly
obtained from $F(x)$. Since $F(x)=1-\left(N^{x}\right)^{2}$, it can be
observed that $F(x)\leq 1$, and when $F(x)=1$, it represents the boundary of
applicability for $F(x)$. And $x_{min}$ is one of the solutions to $F(x)=1$.

Furthermore, we plotted the figure of the outgoing expansion, which shares
similar properties with the outgoing null geodesics due to their expressions
differing only by a factor of $2/x$. In addition, we also present the plot of
surface gravity, it is observed that in the dS spacetime, as the black hole
mass increases, $\kappa$ rapidly decreases and approaches zero, whereas in the
AdS spacetime, $\kappa$ continues to slowly increase.

###### Acknowledgements.

This work is supported by National Natural Science Foundation of China (NSFC)
with Grants No.12275087.

## References

  * [1] A. Ashtekar and J. Lewandowski, Background Independent Quantum Gravity: A Status Report, Class. Quantum Grav. 21, R53 (2004).
  * [2] C. Rovelli, Quantum Gravity (Cambridge University Press, Cambridge, England, 2004).
  * [3] T. Thiemann, Modern Canonical Quantum General Relativity (Cambridge University Press, Cambridge, United Kingdom, 2007).
  * [4] M. Han, Y. Ma, and W. Huang, Fundamental structure of loop quantum gravity, Int. J. Mod. Phys. D 16, 1397 (2007).
  * [5] C. Rovelli and L. Smolin, Discreteness of Area and Volume in Quantum Gravity, Nuclear Physics B 442, 593 (1995).
  * [6] A. Ashtekar and J. Lewandowski, Quantum Theory of Geometry: I. Area Operators, Class. Quantum Grav. 14, A55 (1997).
  * [7] A. Ashtekar and J. Lewandowski, Quantum Theory of Geometry II: Volume Operators, Adv. Theor. Math. Phys. 1, 388 (1997).
  * [8] J. Yang and Y. Ma, New Volume and Inverse Volume Operators for Loop Quantum Gravity, Phys. Rev. D 94, 044003 (2016).
  * [9] T. Thiemann, A Length Operator for Canonical Quantum Gravity, Journal of Mathematical Physics 39, 3372 (1998).
  * [10] Y. Ma, C. Soo, and J. Yang, New Length Operator for Loop Quantum Gravity, Phys. Rev. D 81, 124026 (2010).
  * [11] X. Zhang and Y. Ma, Extension of Loop Quantum Gravity to f (R) Theories, Phys. Rev. Lett. 106, 171301 (2011).
  * [12] N. Bodendorfer, T. Thiemann, and A. Thurn, New Variables for Classical and Quantum Gravity in All Dimensions I. Hamiltonian Analysis, Class. Quantum Grav. 30, 045001 (2013).
  * [13] X. Zhang, J. Yang, and Y. Ma, Canonical Loop Quantization of the Lowest-Order Projectable Horava Gravity, Phys. Rev. D 102, 124060 (2020).
  * [14] H. Tan and Y. Ma, Semiclassical States in Homogeneous Loop Quantum Cosmology, Class. Quantum Grav. 23, 6793 (2006).
  * [15] M. Bojowald, Absence of a Singularity in Loop Quantum Cosmology, Phys. Rev. Lett. 86, 5227 (2001).
  * [16] L. Modesto, Disappearance of the Black Hole Singularity in Loop Quantum Gravity, Phys. Rev. D 70, 124009 (2004).
  * [17] A. Ashtekar and M. Bojowald, Quantum Geometry and the Schwarzschild Singularity, Class. Quantum Grav. 23, 391 (2005).
  * [18] A. Ashtekar, M. Bojowald, and J. Lewandowski, Mathematical Structure of Loop Quantum Cosmology, Advances in Theoretical and Mathematical Physics 7, 233 (2003).
  * [19] M. Bojowald, Loop Quantum Cosmology, Living Rev. Relativ. 11, 4 (2008).
  * [20] A. Ashtekar and P. Singh, Loop Quantum Cosmology: A Status Report, Class. Quantum Grav. 28, 213001 (2011).
  * [21] C. G. Boehmer and K. Vandersloot, Loop Quantum Dynamics of the Schwarzschild Interior, Phys. Rev. D 76, 104030 (2007).
  * [22] D.-W. Chiou, Phenomenological Loop Quantum Geometry of the Schwarzschild Black Hole, Phys. Rev. D 78, 064040 (2008).
  * [23] A. Ashtekar, J. Olmedo, and P. Singh, Quantum Transfiguration of Kruskal Black Holes, Phys. Rev. Lett. 121, 241301 (2018).
  * [24] C. Zhang and X. Zhang, Quantum geometry and effective dynamics of Janis-Newman-Winicour singularities, Phys. Rev. D 101, 086002 (2020),
  * [25] J. Lewandowski, Y. Ma, J. Yang, and C. Zhang, Quantum Oppenheimer-Snyder and Swiss Cheese Models, Phys. Rev. Lett. 130, 101501 (2023).
  * [26] X. Zhang, Loop Quantum Black Hole, Universe 9, 313 (2023).
  * [27] J. G. Kelly, R. Santacruz, and E. Wilson-Ewing, Effective Loop Quantum Gravity Framework for Vacuum Spherically Symmetric Spacetimes, Phys. Rev. D 102, 106024 (2020).
  * [28] R. Gambini, J. Olmedo, and J. Pullin, Spherically symmetric loop quantum gravity: Analysis of improved dynamics, Classical Quantum Gravity 37, 205012 (2020).


# Providing Safety Assurances for Systems with Unknown Dynamics

Hao Wang1, Javier Borquez1 and Somil Bansal1 1The authors are associated with
the Ming Hsieh Department of Electrical and Computer Engineering, University
of Southern California. {haowwang, javierbo<EMAIL_ADDRESS>work is
supported in part by the NSF CAREER Program under award 2240163, the DARPA
ANSR program, a Ming Hsieh Department of Electrical and Computer Engineering
Fellowship.

###### Abstract

As autonomous systems become more complex and integral in our society, the
need to accurately model and safely control these systems has increased
significantly. In the past decade, there has been tremendous success in using
deep learning techniques to model and control systems that are difficult to
model using first principles. However, providing safety assurances for such
systems remains difficult, partially due to the uncertainty in the learned
model. In this work, we aim to provide safety assurances for systems whose
dynamics are not readily derived from first principles and, hence, are more
advantageous to be learned using deep learning techniques. Given the system of
interest and safety constraints, we learn an ensemble model of the system
dynamics from data. Leveraging ensemble uncertainty as a measure of
uncertainty in the learned dynamics model, we compute a maximal robust control
invariant set, starting from which the system is guaranteed to satisfy the
safety constraints under the condition that realized model uncertainties are
contained in the predefined set of admissible model uncertainty. We
demonstrate the effectiveness of our method using a simulated case study with
an inverted pendulum and a hardware experiment with a TurtleBot. The
experiments show that our method robustifies the control actions of the system
against model uncertainty and generates safe behaviors without being overly
restrictive. The codes and accompanying videos can be found on the project
website 111https://github.com/haowwang/safety_assurances_for_unknown_dynamics.

##  I Introduction

Autonomous systems are playing increasingly important roles in the functioning
of modern society. However, traditional modeling techniques, such as using
first principles, struggle to model these systems, given their increasing
complexities. Recent advances in deep learning have enabled the modeling and
control of systems with highly complex dynamics. While the methods have
demonstrated strong control performance for a number of autonomous systems,
they can lead to unsafe behaviors or even catastrophic failures due to the
predictive uncertainty in the neural network models.

In this work, we are interested in providing safety assurances for systems
whose dynamics are unknown and difficult to model using first principles.
Despite success in safety analysis for models developed from first principles,
it remains difficult to provide safety assurances for systems with unknown or
uncertain dynamics. Many works have utilized safety analysis frameworks, such
as Control Barrier Function (CBF) [2] and Hamilton-Jacobi (HJ) reachability
analysis [4], to provide safety assurances to systems with unknown or
uncertain dynamics. One popular line of works seeks to reduce model
uncertainties to uncertainties in the CBF constraints, which are used to
synthesize safety-critical controls [8, 20, 19]. However, these methods rely
on prior knowledge of the system to construct the CBFs, and more critically
assume the CBFs constructed for the nominal model is valid for the actual
system. This assumption can be easily violated when the model uncertainty is
severe, and it is difficult to determine when the assumption is no longer
valid. Another family of approaches is to generate safety guarantees for
systems under the _worst-case_ model uncertainty using HJ reachability
analysis. This line of work is rooted in differential game theory and its
connection with the Hamilton-Jacobi-Isaacs partial differential equation (HJI
PDE) [6, 10]. A time-dependent HJI PDE [18] is formulated to study the
pursuer-evader game, a type of two-person zero-sum differential game integral
to HJ reachability analysis, and this formulation has been applied to generate
safety assurances for system under dynamics uncertainties. More specifically,
the authors in [11, 14, 1] represent the model uncertainty as a Gaussian
Process (GP) and utilize the predictive uncertainty of the GP to generate
robust safety assurances online. However, the model uncertainties considered
in the works depend only on the state and not the control input. Moreover, it
is not immediately clear how similar techniques can be used for neural
network-based dynamics models that does not predict uncertainties.

In its core, our method attempts to provide robust safety assurances for
systems with unknown dynamics against model uncertainties. We first learn a
nominal dynamics model of the system along with a measure of model
uncertainty. Then, we compute robust safety assurances utilizing results from
HJ reachability analysis. More specifically, given a set of undesirable
states, our method computes the maximal robust control invariant set and a
controller that renders the set control invariant for the model under bounded
model uncertainty. A critical aspect of our work is that model uncertainty is
both _state and control-dependent_. This control dependence of model
uncertainty is often ignored in prior works for simplicity. Specifically, when
the model uncertainty is control-independent, the players in the two-person
game do not interact directly and renders the game much easier to solve.
Instead, we take this interplay between the control and model uncertainty into
account and systematically handle that under the HJ reachability framework.
This leads to a significantly less conservative estimate of the safe set of
the system. Furthermore, we pay special attention to learned neural network
nominal models. Using an ensemble of neural networks with control-affine
architectures as the nominal model, we harness the modeling power of neural
networks and obtain a heuristic measure of model uncertainty from the
ensemble, with which we provide robust safety assurances for the system.

Our method is shown to provide significant benefits as it explicitly
incorporates state and control-dependent model uncertainty in safety analysis
and, as a result, provides more robust, but not overly conservative
controllers for the system. To summarize, the key contributions of this letter
are two-fold: 1) We propose a framework for providing robust safety assurances
to systems with unknown dynamics by incorporating state and control-dependent
model uncertainty in HJ reachability analysis. We also provide a closed-form
solution to the two-player zero-sum differential game between the control and
disturbance player, in order to accommodate the state and control dependent
nature of the model uncertainty, and 2) We provide a concrete instantiation of
our framework by modeling the system with an ensemble of neural networks and
solving a robust optimal control problem to safeguard the learned model
against bounded model uncertainty.

##  II Problem Formulation

In this work, we are interested in providing safety assurances for systems
with unknown but deterministic _control-affine_ dynamics, as many practical
systems are control-affine [16], governed by ordinary differential equation

$\frac{dx}{dt}=\dot{x}=f(x,u)=f_{1}(x)+f_{2}(x)u$ (1)

where $x\in\mathcal{X}\subseteq\mathbb{R}^{n_{x}}$ and
$u\in\mathcal{U}\subseteq\mathbb{R}^{n_{u}}$ are the state and control of the
system, and $f_{1}(x)$ and $f_{2}(x)$ are matrices of appropriate sizes. Given
a set of states $\mathcal{F}\subset\mathcal{X}$, which we refer to as the
_failure set_ , that the system must avoid (e.g., obstacles for a mobile
robot), our goal is ensuring the system $f$ does not enter the failure set
$\mathcal{F}$ for predefined time horizon $[0,T]$. More formally as state in
Prob. 2.

######  Problem 1 (Safety Problem for System $f$).

Given a failure set $\mathcal{F}$ and time horizon $[0,T]$, obtain a safe set
$\mathcal{S}\subseteq\mathcal{X}$, along with a state-feedback controller
$\pi$, such that the system $f$ starting within $\mathcal{S}$ implies that $f$
does not enter $\mathcal{F}$ for $[0,T]$ using controls provided by $\pi$. A
safe set is the maximal safe set, denoted by $\mathcal{S}^{*}$, if it contains
all other safe sets as subsets.

Because we do not have access to the system dynamics $f$, we further limit the
scope of the problem to employing a model-based approach - first model the
system of interest, then provide safety assurances on the system against any
_bounded_ model uncertainty $d\in\mathcal{D}$, where $\mathcal{D}$ is the set
of admissible model uncertainty. We state the proposed safety assurance
problem as follows:

######  Problem 2 (Safety Assurance under Model Uncertainty).

Given a failure set $\mathcal{F}$ and time horizon $[0,T]$, obtain a set
$\mathcal{R}\subseteq\mathcal{X}$, along with a state-feedback controller
$\pi_{\mathcal{R}}$, such that under the condition that all realized model
uncertainties are contained within $\mathcal{D}$, system $f$ starts within
$\hat{\mathcal{S}}$ implies that $f$ does not enter $\mathcal{F}$ for $[0,T]$
using controls from $\hat{\pi}$, or equivalently,
$\mathcal{R}\subseteq\mathcal{S}^{*}$.

######  Definition 1.

(Robust Safe Set) A solution to Problem 2, $\mathcal{R}\subseteq\mathcal{X}$,
is referred to as a robust safe set. A robust safe set is the maximal robust
safe set, denoted as $\mathcal{R}^{*}$, is the robust safe set that contains
any robust safe set as subsets.

##  III Background

In this section, we provide a brief overview of Hamtilon-Jacobi (HJ)
reachability analysis, an approach that can help solve Prob. 2. Let $g$ be a
system described by dynamics $\dot{x}=g(x,u,d)$, where
$x\in\mathcal{X}\subseteq\mathbb{R}^{n_{x}}$,
$u\in\mathcal{U}\subseteq\mathbb{R}^{n_{u}}$, and
$d\in\mathcal{D}\subseteq\mathbb{R}^{n_{d}}$ are the state, control, and
disturbance of the system. Disturbance $d$, in the interest of this letter,
describes the model uncertainty.We use
$\mathbf{u}:[0,T]\rightarrow\mathcal{U}$ and
$\mathbf{d}:[0,T]\rightarrow\mathbb{R}^{n_{d}}$ to denote the control and
disturbance signals. Furthermore, we denote the state trajectory starting from
state $x$ at time $t$ evolved with control and disturbance signals
$\mathbf{u}(\cdot)$ and $\mathbf{d}(\cdot)$ as
$\xi_{x,t}^{\mathbf{u},\mathbf{d}}$. With a slight abuse of the notation, we
use $\xi_{x,t}^{\mathbf{u},\mathbf{d}}(\tau)$ to denote the state at time
$\tau\geq t$ along the trajectory $\xi_{x,t}^{\mathbf{u},\mathbf{d}}$.

Suppose $g(x,u,d)$ is uniformly continuous, bounded, and Lipschitz continuous
in $x$ for fixed $u$ and $d$. We further assume that $\mathcal{U}$ and
$\mathcal{D}$ are compact, and $\mathbf{u}(\cdot)$ as well as
$\mathbf{d}(\cdot)$ are measurable. Let $l(x)$ be the signed distance function
to $\mathcal{F}$. We can obtain $\mathcal{R}^{*}$ by solving the following
robust optimal control problem (Prob. 3) with initial condition $x=x_{0}\
\forall x_{0}\in\mathcal{X}$ and $t=0$. In robust control literature, Prob. 3
is posed as a two-player zero-sum differential game between the control $u$
and disturbance $d$, who uses only nonanticipative strategies [18]. Let
$\Gamma(t)$ be the set of nonanticipative strategies, and $\mathbb{U}(t)$ be
the set of admisable control signals.

######  Problem 3.

$\begin{split}\inf_{\mathbf{d}(\cdot)\in\Gamma(t)}\sup_{\mathbf{u}(\cdot)\in\mathbb{U}(t)}&\quad
J(x,t,\mathbf{u},\mathbf{d})=\min_{\tau\in[t,T]}l(\xi_{x,t}^{\mathbf{u},\mathbf{d}}(\tau))\\\
s.t.&\quad\dot{x}=g(x,u,d)\end{split}$

Let us define the value function $V(x,t)$ to take on the optimal value of
Prob. 3 at state $x$ and time $t$:

$\begin{split}V(x,t)=&\inf_{\mathbf{d}(\cdot)\in\Gamma(t)}\sup_{\mathbf{u}(\cdot)\in\mathbb{U}(t)}J(x,t,\mathbf{u},\mathbf{d})\\\
&=\inf_{\mathbf{d}(\cdot)\in\Gamma(t)}\sup_{\mathbf{u}(\cdot)\in\mathbb{U}(t)}\min_{\tau\in[t,T]}l(\xi_{x,t}^{\mathbf{u},\mathbf{d}}(\tau))\end{split}$
(2)

Then, $\mathcal{R}^{*}$ can be characterized using $V(x,t)$:

$\mathcal{R}^{*}=\\{x\in\mathcal{X}|V(x,0)>0\\}$ (3)

HJ reachability analysis provides a tractable means to compute the value
function $V(x,t)$. It has been shown that $V(x,t)$ is the unique viscosity
solution of the Hamilton-Jacobi-Issacs Variational Inequality (HJI-VI) [17,
18]:

$\displaystyle\min\\{$ $\displaystyle D_{t}V+H(x,t,\nabla V),l(x)-V(x,t)\\}=0$
(4) $\displaystyle V(x,T)=l(x),\quad\text{for}\ t\in\left[0,T\right]$

$H(x,t,\nabla V)=\max_{u}\min_{d}\langle\nabla V,g(x,u,d)\rangle$ is the
Hamiltonian. $D_{t}V$ and $\nabla V$ denote the temporal derivative and the
spatial gradients of $V(x,t)$. It is important to note that HJ reachability
analysis also provides a state-feedback controller $\pi_{\mathcal{R}^{*}}$
that renders $\mathcal{R}^{*}$ control-invariant:

${\color[rgb]{0,0,0}\pi_{\mathcal{R}^{*}}(x,t)=\arg\max_{u\in\mathcal{U}}\min_{d\in\mathcal{D}}\langle\nabla
V(x,t),g(x,u,d)\rangle}$ (5)

##  IV Safety Assurances for Learned Dynamics

At the heart of our framework is solving a robust optimal control problem
(Prob. 3) to provide safety assurances for system $f$, by incorporating the
worst-case model uncertainty in the safety analysis. We obtain the maximal
robust safe set $\mathcal{R}^{*}$ that safeguards the system against model
uncertainty by solving the HJI-VI (4) for $V(x,t)$. In this section, we first
introduce the notion of _uncertain model_ and use it to set up the robust
optimal control problem. Then, we show $\mathcal{R}^{*}$ does in fact confer
safety assurance to $f$. Finally, we present a concrete instantiation of our
framework with $f$ modeled by an ensemble of neural networks, and we discuss a
method to quantify the model uncertainty from the ensemble.

###  IV-A Model Representation and Hamiltonian Formulation

Given a deterministic, continuous-time, control-affine system $f(x,u)$, we
learn a nominal model $\bar{f}(x,u)=\bar{f_{1}}(x)+\bar{f_{2}}(x)u$ of $f$. We
assume there are _bounded, state-dependent_ model uncertainties
$d_{1}(x)\in\mathcal{D}_{1}(x)\subseteq\mathbb{R}^{n_{x}}$ and
$d_{2}(x)\in\mathcal{D}_{2}(x)\subseteq\mathbb{R}^{n_{x}\times n_{u}}$,
arising from errors of model approximation and unmodeled dynamics additive to
$\bar{f_{1}}(x)$ and $\bar{f_{2}}(x)$. Hence, the _uncertain_ model
$\hat{f}(x,u,d_{1},d_{2})$ can be written as

$\hat{f}(x,u,d_{1},d_{2})=\bar{f_{1}}(x)+d_{1}(x)+(\bar{f_{2}}(x)+d_{2}(x))u$
(6)

We refer to $\mathcal{D}_{1}(x)$ and $\mathcal{D}_{2}(x)$ as the _model
uncertainty bounds_ on $\bar{f}_{1}(x)$ and $\bar{f}_{2}(x)$ at state $x$,
respectively. The Hamiltonian $H(x,t,\nabla V)$ formulated using
$\hat{f}(x,u,d_{1},d_{2})$ is given in (7), where $u^{*},d_{1}^{*}$, and
$d_{2}^{*}$ are the solutions to the maximin game in the RHS of (7a).

$\begin{split}H(x,&t,\nabla
V(x,t))=\max_{u\in\mathcal{U}}\min_{d_{1}\in\mathcal{D}_{1}(x)}\min_{d_{2}\in\mathcal{D}_{2}(x)}\Bigl{\langle}\nabla
V(x,t),\\\
&\bar{f}_{1}(x)+d_{1}+\left(\bar{f}_{2}(x)+d_{2}\right)u\Bigl{\rangle}\\\
\end{split}$ (7a) $=\Bigl{\langle}\nabla
V(x,t),\bar{f}_{1}(x)+d_{1}^{*}+\left(\bar{f}_{2}(x)+d_{2}^{*}\right)u^{*}\Bigl{\rangle}$
(7b)

The maximin game in (7a) does not generally have a closed-form solution since
$u$ and $d_{2}$ are multiplied together. We make the following assumptions to
enable tractable computation of (7a): 1) $\mathcal{D}_{1}(x)$ and
$\mathcal{D}_{2}(x)$ are hypercubes containing the origin in their respective
spaces $\mathbb{R}^{n_{x}}$ and $\mathbb{R}^{n_{x}\times n_{u}}$, and 2) the
set of admissible controls $\mathcal{U}$ is a hypercube containing the origin
in $\mathbb{R}^{n_{u}}$. Under the assumptions, we can obtain closed-form
solutions for $u$, $d_{1}$ and $d_{2}$.

######  Note 1.

The total model uncertainty is given by $d(x,u)=d_{1}(x)+d_{2}(x)u$, and hence
the representation of $d$ is state and control-dependent.

Let $d_{1i}^{*}(x)$ and $\nabla V_{i}(x,t)$ denote the $i^{th}$ component of
$d_{1}^{*}(x)\in\mathbb{R}^{n_{x}}$ and $\nabla V(x,t)\in\mathbb{R}^{n_{x}}$,
respectively. Since $\mathcal{D}_{1}(x)$ is a hypercube containing the origin
of $\mathbb{R}^{n_{x}}$, we can write the uncertainty bound on the $i^{th}$
component of $\bar{f}_{1}(x)$ by an interval
$\left[\underline{d_{1i}(x)},\overline{d_{1i}(x)}\right]$. Then,
$d_{1}^{*}(x)$ is given by

$d_{1i}(x)=\begin{cases}\overline{d_{1i}(x)}\ \text{if}\ \nabla
V_{i}(x,t)<0\\\ \underline{d_{1i}(x)}\ \text{if}\ \nabla V_{i}(x,t)\geq
0\end{cases}$ (8)

We now derive $u^{*}(x)$, the optimal safety controller given in (9). Let us
denote the $j^{th}$ component of $u^{*}\in\mathbb{R}^{n_{u}}$ and the $j^{th}$
column of $\bar{f_{2}}(x)$ by $u_{j}^{*}$ and $\bar{f}_{2j}(x)$. Since
$\mathcal{U}$ is a hypercube in $\mathbb{R}^{n_{u}}$, the bound on the
$j^{th}$ component of $u$ is given by an interval
$\left[\underline{u_{j}},\overline{u_{j}}\right]$. Again, since
$\mathcal{D}_{2}(x)$ is a hypercube in $\mathbb{R}^{n_{x}\times n_{u}}$, the
model uncertainty bound on the $ij^{th}$ component of $\bar{f}_{2}(x)$ is an
interval $\left[\underline{d_{2ij}(x)},\overline{d_{2ij}(x)}\right]$. Let
$d_{2}^{+}(x)\in\mathbb{R}^{n_{x}\times n_{u}}$ and
$d_{2}^{-}(x)\in\mathbb{R}^{n_{x}\times n_{u}}$ be the “best effort”
$d_{2}(x)$ that intuitively try to decrease the Hamiltonian for positive or
negative $u$, respectively. We denote the $j^{th}$ column of $d_{2}^{+}(x)$
and $d_{2}^{-}(x)$ by $d_{2j}^{+}(x)$ and $d_{2j}^{-}(x)$. More precisely, the
$i^{th}$ component of $d_{2j}^{+}(x)\in\mathbb{R}^{n_{x}}$ is given by
$\overline{d_{2ij}}$ when $\nabla V_{i}<0$, $\underline{d_{2ij}}$ when $\nabla
V_{i}\geq 0$. The $i^{th}$ component of $d_{2j}^{-}(x)\in\mathbb{R}^{n_{x}}$
is given by $\overline{d_{2ij}}$ when $\nabla V_{i}\geq 0$,
$\underline{d_{2ij}}$ when $\nabla V_{i}<0$. $u^{*}$ is then given in (9).

$u^{*}_{j}(x)=\begin{cases}\vspace{0.1em}\overline{u_{j}},\text{if}\ \nabla
V(x,t)^{\top}\left(\bar{f}_{2j}(x)+d_{2j}^{+}(x)\right)>0\\\
\underline{u_{j}},\text{if}\ \nabla
V(x,t)^{\top}\left(\bar{f}_{2j}(x)+d_{2j}^{-}(x)\right)<0\\\
0,\text{otherwise}\end{cases}$ (9)

Lastly, we provide $d_{2}^{*}(x)$. Let us denote the $ij^{th}$ component of
$d_{2}^{*}(x)\in\mathbb{R}^{n_{x}\times n_{u}}$ by $d_{2ij}^{*}(x)$.
$d_{2}^{*}(x)$ is given below in (10).

$d_{2ij}^{*}(x)=\begin{cases}\overline{d_{2ij}(x)}\ \text{if}\
u_{j}^{*}(x)\nabla V_{i}(x,t)<0\\\ \underline{d_{2ij}(x)}\ \text{if}\
u_{j}^{*}(x)\nabla V_{i}(x,t)\geq 0\end{cases}$ (10)

We solve the HJI-VI (4), with the Hamiltonian (7b) formulated for $\hat{f}$,
to obtain $V(x,t)$ and the maximal robust safe set $\mathcal{R}^{*}$ as well
as its corresponding safety controller $\pi_{\mathcal{R}^{*}}$. Following
directly from the definition of the value function (2), $\mathcal{R}^{*}$ and
$\pi_{\mathcal{R}^{*}}$ confer safety assurances to system $f$, if for any
state $x\in\mathcal{X}$, the realized model uncertainties at state $x$ are
contained within $\mathcal{D}_{1}(x)$ and $\mathcal{D}_{2}(x)$. This result is
formalized in Lemma. 1.

######  Lemma 1.

Given a control-affine system $f(x,u)=f_{1}(x)+f_{2}(x)u$ and its nominal
model $\bar{f}(x,u)=\bar{f}_{1}(x)+\bar{f}_{2}(x)u$, let
$\delta_{1}(x)=f_{1}(x)-\bar{f}_{1}(x)$ and
$\delta_{2}(x)=f_{2}(x)-\bar{f}_{2}(x)$ be the realized model uncertainties.
If $\delta_{1}(x)\in\mathcal{D}_{1}(x)$ and
$\delta_{2}(x)\in\mathcal{D}_{2}(x)$ $\forall x\in\mathcal{X}$, then
$\mathcal{R}^{*}\subseteq\mathcal{S}^{*}$.

###  IV-B Ensemble Dynamics Representation

Although our framework is agnostic to how the nominal model $\bar{f}$ or the
model uncertainties $d_{1}$ and $d_{2}$ are obtained, in this subsection, we
introduce one method to jointly obtain $\bar{f}$, $d_{1}$, $d_{2}$, and their
corresponding bounds $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$.

In this work, we employ an ensemble of neural networks to obtain uncertainty
in learned dynamics, a popular approach in literature to quantify uncertainty
in deep networks [15, 9]. More specifically, given a deterministic,
continuous-time, control-affine dynamical system $f(x,u)$, we learn an
ensemble of $M$ fully connected feed-forward neural networks with control-
affine architectures, and we use the ensemble as the nominal model $\bar{f}$.
More explicitly, the ensemble is given by

$E=\\{NN^{k}(x,u)=NN^{k}_{1}(x)+NN^{k}_{2}(x)u\\}_{k=1}^{N_{m}}$ (11)

where $NN^{k}_{1}$ and $NN^{k}_{2}$ are neural networks. With a slight abuse
of notation, we denote the prediction of the ensemble $E$ by

$\begin{split}E(x,u)&=\frac{1}{N_{m}}\sum_{k=1}^{N_{m}}NN^{k}(x,u)\\\
&=\left(\frac{1}{N_{m}}\sum_{k=1}^{N_{\texttt{m}}}NN^{k}_{1}(x)\right)+\left(\frac{1}{N_{m}}\sum_{k=1}^{N_{m}}NN^{k}_{2}(x)\right)u\\\
&=\bar{f}_{1}(x)+\bar{f}_{2}(x)u=\bar{f}(x,u)\end{split}$ (12)

We sometimes refer to $E(x,u)$ as the _mean dynamics_ , as we are taking the
mean prediction among the neural networks within the ensemble.

Given the setup of the nominal model $\bar{f}(x,u)$ in (12), we would like the
model uncertainties $d_{1}$ and $d_{2}$ to intuitively quantify the
_variations of outputs_ of the sub-nets $NN_{1}^{k}$ and $NN_{2}^{k}$ within
the ensemble. We bound $d_{1}(x)$ and $d_{2}(x)$ using constant multiples of
standard deviation of $\\{NN^{k}_{1}(x)\\}_{k=1}^{N_{m}}$ and
$\\{NN^{k}_{2}(x)\\}_{k=1}^{N_{m}}$. Let us denote the extents of the $i^{th}$
and $ij^{th}$ dimensions of $\mathcal{D}_{1}(x)$ and $\mathcal{D}_{2}(x)$ by
$\mathcal{D}_{1i}(x)$ and $\mathcal{D}_{2ij}(x)$, respectively. By assumption,
$\mathcal{D}_{1}(x)\subseteq\mathbb{R}^{n_{x}}$ and
$\mathcal{D}_{2}(x)\subseteq\mathbb{R}^{n_{x}\times n_{u}}$ are hypercubes in
their respective spaces, and therefore $\mathcal{D}_{1i}(x)$ and
$\mathcal{D}_{2ij}(x)$ are real intervals. More precisely, the intervals are
given by

$\mathcal{D}_{1i}(x)=[-\alpha\sigma_{1i}(x),\alpha\sigma_{1i}(x)]$ (13a)
$\mathcal{D}_{2ij}(x)=[-\gamma\sigma_{2ij}(x),\gamma\sigma_{2ij}(x)]$ (13b)

where
$\sigma_{1i}(x)=\texttt{StdDev}\left(\\{NN^{k}_{1i}(x)\\}_{k=1}^{N_{m}}\right)$,
$\sigma_{2ij}(x)=\texttt{StdDev}\left(\\{NN^{k}_{2ij}(x)\\}_{k=1}^{N_{m}}\right)$,
StdDev is a short hand for “standard deviation”, $NN^{k}_{1i}(x)$ and
$NN^{k}_{2ij}(x)$ denote the $i^{th}$ and $ij^{th}$ outputs of the $k^{th}$
$NN_{1}$ and $NN_{2}$ sub-nets, respectively.

######  Note 2.

$\alpha$ and $\gamma$ are tunable parameters that determine conservativeness
of the resulting safe set $\mathcal{S}$. As $\alpha$ and $\gamma$ increase,
the model uncertainty bounds $\mathcal{D}_{1}(x)$ and $\mathcal{D}_{2}(x)$
increases $\forall x\in\mathcal{X}$. Accordingly, the safe set $\mathcal{S}$
shrinks. For all the experiments in this letter, we use $\alpha=\gamma=3$.

##  V Experiments

###  V-A Inverted Pendulum

In this example, we simulate an inverted pendulum with state
$x=[\theta,\dot{\theta}]^{\top}$ and control $u$, and its dynamics is given by
$\dot{x}=\left[\dot{\theta},\ddot{\theta}\right]^{\top}=\left[\dot{\theta},\frac{-b\dot{\theta}+\frac{1}{2}mgl\sin{\theta}-u}{\frac{ml^{2}}{3}}\right]^{\top}$
, where $l,m,g,$ and $b$ represents the length and mass of the pendulum,
acceleration of gravity, and the friction coefficient, respectively. For the
purposes of this case study, the analytical expression of the dynamics is
assumed to be unknown.

We first train an ensemble dynamics model, consisting of 5 fully connected
feed-forward neural networks with 3 hidden layers and 256 neurons per hidden
layer, using dataset $\\{(x_{i},u_{i}),\dot{x_{i}}\\}_{i=1}^{M}$ parsed from
trajectory rollouts with random controls. In this experiment, we want the
pendulum to avoid deviating from its unstable equilibrium for more than
$0.6\pi$ radians. Equivalently, the failure set is given by
$\mathcal{F}=\\{[\theta,\dot{\theta}]^{\top}|\theta>0.6\pi\\}\cup\\{[\theta,\dot{\theta}]^{\top}|\theta<-0.6\pi\\}$.

Next, we compute the safe set $\mathcal{S}$, or equivalently complement of the
backward reachable tube (BRT) of the failure set $\mathcal{F}$, for a time
horizon of 0.7 seconds, using the learned ensemble dynamics model along with
model uncertainties. We also consider three baselines: 1) computing the safe
set using only the mean dynamics $E(x,u)$ (Baseline 1 Mean Dynamics), 2)
computing the safe set with $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ calculated
using split conformal prediction with 5000 calibration samples and marginal
coverage of $95\%$ [3] (Baseline 2 Conformal Prediction), and 3) computing the
safe set with $d_{2}(x)u$ approximated as $d_{3}(x)\in\mathbb{R}^{n_{x}}$
(Baseline 3 Partial Game). The purpose of Baseline 3 is to remove the
interaction of the control player $u$ and the disturbance player $d_{2}$, and
renders the model uncertainty _control-independent_ , as the action of the
disturbance player $d_{3}$ no longer depends on that of the control player
$u$. Furthermore, the Hamiltonian computation (7a) decouples into three
independent optimizations with respect to $u,d_{1}$, and $d_{3}$. The bound on
the $i^{th}$ component of $d_{3}$ is given by
$\mathcal{D}_{3}(x)_{i}=[-a(x),a(x)]$, where
$a(x)=\sum_{j=1}^{n_{u}}\max\\{|\underline{u_{j}}|,|\overline{u_{j}}|\\}\times\max\\{|\underline{d_{2ij}(x)}|,|\overline{d_{2ij}(x)}|\\}$.

For the number of training samples $M=300$, we visualize the recovered safe
sets for our method, all the baselines, and the ground truth in Fig. 1. The
ground truth safe set, computed using the analytical expression of the
dynamics, is shaded in green. The safe set recovered using our method is
entirely contained within the ground truth safe set, indicating satisfaction
of the safety constraint. On the other hand, Baseline 1 fails to satisfy the
safety constraint, since the recovered safe set is not contained within the
ground truth safe set. The safe set from Baseline 2 is not visualized in Fig.
1 because its recovered safe set is empty due to model uncertaintity bounds
$\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ being too conservative. Specifically,
conformal prediction provides _state-independent_ uncertainty bounds, which
are dictated by the worst-case modeling errors across all states, leading to
overly conservative behaviors. Baseline 3 is less conservative than Baseline
2, as its model uncertainty bounds are state-dependent, but it is more
conservative than our methods since its model uncertainty $d_{3}(x)$ is
control-independent and its bound $\mathcal{D}_{3}(x)$ is an overapproximation
of that of our method.

Figure 1: Recovered safe sets (states starting from which the inverted
pendulum stays within $0.6\pi$ of the upright for 0.7 seconds) for our method,
the Mean Dynamics baseline (Baseline 1), the Partial Game baseline (Baseline
3) in the inverted pendulum experiment, with the ensemble trained with 300
training samples. Note that Baseline 2 (conformal prediction) is not
visualized, because its safe set is empty.

We also perform an ablation study over the number of training samples $M$ to
further highlight the benefit of using state and control-dependent model
uncertainty. The percent safe set recovered, as a function of $M$, is charted
in Fig. 2. Across all experimented $M$, our method consistently outperforms
the baselines, indicating that the state and control-dependent model
uncertainty representation leads to less conservative behaviors. Our model
uncertainty representation can reflect local variations of model
uncertainties, allowing the safe set to expand or shrink according to local
model uncertainty level. Furthermore, the control-dependent nature of our
model uncertainty representation, which taken into consideration of the
interaction between the control and the model uncertainty, also helps reduce
the conservativeness.

Figure 2: Changes in percent safe set recovered with our method, the Conformal
Prediction baseline (Baseline 2), and the Partial Game baseline (Baseline 3),
over the number of training samples.

###  V-B Turtlebot Hardware Experiment

We apply the proposed approach on a real hardware testbed, TurtleBot 4, which
we refer to as the _vehicle_ , in this experiment. We are interested in
providing safety assurances for the vehicle _carrying a payload_ , and we seek
to demonstrate the importance of incorporating model uncertainty in safety
analysis and the benefit of model learning.

Let $\mathcal{A}\subseteq\mathbb{R}^{2}$ be a rectangular experimental space,
centered at $\left[0,0\right]^{\top}$ with side lengths $l_{x}$ and $l_{y}$.
The failure set $\mathcal{F}$ is hence given by
$\mathcal{F}=\\{x\in\mathcal{X}||p_{x}|>\frac{l_{x}}{2},|p_{y}|>\frac{l_{y}}{2}\\}$,
where $p_{x}$ and $p_{y}$ are the $x,y$ positions of the center of mass of the
vehicle. We model the vehicle as a 4-dimensional system with state
$\left[p_{x},p_{y},\theta,\omega\right]^{\top}$, where $\omega$ is the angular
velocity, and control $u$, an angular velocity setpoint. We operate the
vehicle with a constant forward velocity.

We first collect 120 state and control trajectories using a random control at
each time step. Each trajectory is roughly 5 seconds and yields about 100
training samples $(\left[x,u\right],\dot{x})$. Then, we train an ensemble
dynamics model containing 5 fully connected feed-forward neural networks, each
with 3 hidden layers and 512 neurons per layer.

The safe set and the safety controller are computed for our method and the
mean dynamics baseline until convergence (i.e., the time horizon is
$\left[0,\infty\right]$). For visualization purposes, we project the safe set
to the $x,y$ plane at some fixed $\theta$ and $\omega$. In Fig. 3, the safe
sets projected at $\theta=0.7$ and $\omega=0$ are shown. We apply mean
dynamics baseline’s safety controller $\pi^{*}_{m}$, from two states $x_{1}\
\text{and}\ x_{2}$, within the mean dynamics baseline safe set
$\mathcal{S}_{m}$. $\pi^{*}_{m}$ is unable to keep the vehicle inside
$\mathcal{A}$, since $\mathcal{S}_{m}$ and $\pi^{*}_{m}$ do not take into
consideration the model uncertainty and, as a result, are overly optimistic.
On the other hand, we roll out the vehicle from 2 nearby states, $x_{3}\
\text{and}\ x_{4}$, within our method’s safe set $\mathcal{S}$ with our
method’s safety controller $\pi^{*}$, and the vehicle stays within
$\mathcal{A}$, indicating that our method is able to obtain a better estimate
of the actual safe set of the system.

Figure 3: Comparison between projected safe sets $\mathcal{S}$ and
$\mathcal{S}_{m}$ along with TurtleBot rollout trajectories using safety
controllers $\pi^{*}$ and $\pi^{*}_{m}$. The initial state of rollouts using
$\pi^{*}$ are marked with pentagons. All rollouts start with heading
$\theta=0.7$ and angular velocity $\omega=0$. The states at which the
TurtleBot enters failure set $\mathcal{F}$ are marked with crosses.

We now highlight the benefit of model learning with a safety filtering
experiment. When the vehicle is traveling at a constant forward velocity, it
is common to model the vehicle with a three-dimensional Dubins Car (Dubins3D)
with the dynamics
$\left[\dot{p_{x}},\dot{p_{y}},\dot{\theta}\right]^{\top}=\left[V\cos(\theta),V\sin(\theta),u\right]^{\top}$.
However, since the vehicle carries a payload, which introduces factors that
could render the Dubins3D model inaccurate, a learned model can more
accurately represent the system. We compute the safe set $\mathcal{S}_{d}$ and
safety controller $\pi^{*}_{d}$ using the Dubins3D model to convergence. Then,
we use $\pi^{*}_{d}$ to _filter_ a nominal controller $\pi(x)=0$ in a least
restrictive fashion [7]. The filtered controller $\hat{\pi}_{d}$ keeps the
vehicle moving with the current heading unless it is at risk of exiting
$\mathcal{S}_{d}$, in which case $\pi^{*}_{d}$ takes over (i.e.
$\hat{\pi}_{d}(x)=\pi^{*}_{d}(x)$). We similarly filter $\pi$ with $\pi^{*}$,
and filtered controller is denoted as $\hat{\pi}$.

Figure 4: Filtered controller $\hat{\pi}$ and $\hat{\pi}_{d}$ rollout
trajectories. The blue markers indicate the states at which the nominal
controller $\pi$ is on. Whereas the mint markers indicates the states where
safety controllers $\pi^{*}$ and $\pi^{*}_{d}$ intervene. The state at which
the vehicle exits the experiment space under $\hat{\pi}_{d}$ is marked with a
red cross.

Starting from a state $x\in\mathcal{S}\cap\mathcal{S}_{d}$, we roll out the
vehicle with $\hat{\pi}_{d}$ and $\hat{\pi}$. The state trajectories projected
to the $x-y$ plane are shown in Fig. 4. The vehicle eventually exits
$\mathcal{A}$ under $\hat{\pi}_{d}$, indicating that $\pi^{*}_{d}$ does not
actually render $\mathcal{S}_{d}$ control invariant under the actual vehicle
dynamics and is overly optimistic. However, $\hat{\pi}$ does keep the vehicle
inside $\mathcal{A}$ for the entire experiment. For both trajectories, the
mint-colored markers indicate the states where the safety controller $\pi^{*}$
or $\pi^{*}_{d}$ takes over and commands the vehicle to turn maximally to stay
within $\mathcal{A}$. Since $\mathcal{S}_{d}$ is more optimistic than
$\mathcal{S}$, $\pi^{*}_{d}$ intervenes closer to the boundary of
$\mathcal{A}$ than $\pi^{*}$ would, and the resulting trajectory is
uncomfortably close to exiting $\mathcal{A}$. In contrast, there is a healthy
margin for the trajectory filtered by $\pi^{*}$.

##  VI Conclusion

In this letter, we proposed a framework for generating robust safety
assurances for systems with unknown dynamics. Further, we provide a concrete
instantiation of our framework with ensemble neural network models as the
nominal model and safeguard the system against the worst-case model
uncertainty. Though our method is shown to provide robust safety assurances in
the experiments, it faces several challenges, which we look to address in
future works. First of all, our method does not scale well to higher-
dimensional systems. We will investigate the possibility of incorporating
learning-based reachability computation tools [5, 12] into our framework.
Second, the proposed model uncertainty estimation approach might not provide
model uncertainty bounds that accurately reflect the distribution of realized
model uncertainties. We will address this challenge by examining other
uncertainty estimation approaches, such as [15] and [13].

## ACKNOWLEDGMENT

We would like to thank Albert Lin for insightful discussions during various
phases of this project.

## References

  * [1] Anayo K Akametalu, Jaime F Fisac, Jeremy H Gillula, Shahab Kaynama, Melanie N Zeilinger, and Claire J Tomlin.  Reachability-based safe learning with gaussian processes.  In 53rd IEEE Conference on Decision and Control, pages 1424–1431. IEEE, 2014.
  * [2] Aaron D. Ames, Xiangru Xu, Jessy W. Grizzle, and Paulo Tabuada.  Control barrier function based quadratic programs for safety critical systems.  IEEE Transactions on Automatic Control, 62(8):3861–3876, 2017.
  * [3] Anastasios N Angelopoulos and Stephen Bates.  A gentle introduction to conformal prediction and distribution-free uncertainty quantification.  arXiv preprint arXiv:2107.07511, 2021.
  * [4] Somil Bansal, Mo Chen, Sylvia Herbert, and Claire J. Tomlin.  Hamilton-jacobi reachability: A brief overview and recent advances.  In 2017 IEEE 56th Annual Conference on Decision and Control (CDC), pages 2242–2253, 2017.
  * [5] Somil Bansal and Claire J. Tomlin.  Deepreach: A deep learning approach to high-dimensional reachability.  In 2021 IEEE International Conference on Robotics and Automation (ICRA), pages 1817–1824, 2021.
  * [6] Emmanuel Nicholas Barron, Lawrence Craig Evans, and Robert Jensen.  Viscosity solutions of isaacs’ equations and differential games with lipschitz controls.  Journal of Differential Equations, 53(2):213–233, 1984.
  * [7] Javier Borquez, Kaustav Chakraborty, Hao Wang, and Somil Bansal.  On safety and liveness filtering using hamilton-jacobi reachability analysis.  arXiv preprint arXiv:2312.15347, 2023.
  * [8] Fernando Castañeda, Jason J. Choi, Bike Zhang, Claire J. Tomlin, and Koushil Sreenath.  Pointwise feasibility of gaussian process-based safety-critical control under model uncertainty.  In 2021 60th IEEE Conference on Decision and Control (CDC), pages 6762–6769, 2021.
  * [9] Kurtland Chua, Roberto Calandra, Rowan McAllister, and Sergey Levine.  Deep reinforcement learning in a handful of trials using probabilistic dynamics models.  Advances in neural information processing systems, 31, 2018.
  * [10] Lawrence C Evans and Panagiotis E Souganidis.  Differential games and representation formulas for solutions of hamilton-jacobi-isaacs equations.  Indiana University mathematics journal, 33(5):773–797, 1984.
  * [11] Jaime F. Fisac, Anayo K. Akametalu, Melanie N. Zeilinger, Shahab Kaynama, Jeremy Gillula, and Claire J. Tomlin.  A general safety framework for learning-based control in uncertain robotic systems.  IEEE Transactions on Automatic Control, 64(7):2737–2752, 2019.
  * [12] Jaime F. Fisac, Neil F. Lugovoy, Vicenç Rubies-Royo, Shromona Ghosh, and Claire J. Tomlin.  Bridging hamilton-jacobi safety analysis and reinforcement learning.  In 2019 International Conference on Robotics and Automation (ICRA), pages 8550–8556, 2019.
  * [13] Yarin Gal and Zoubin Ghahramani.  Dropout as a bayesian approximation: Representing model uncertainty in deep learning.  In Proceedings of The 33rd International Conference on Machine Learning, volume 48 of Proceedings of Machine Learning Research, pages 1050–1059, New York, New York, USA, 20–22 Jun 2016. PMLR.
  * [14] Sylvia Herbert, Jason J Choi, Suvansh Sanjeev, Marsalis Gibson, Koushil Sreenath, and Claire J Tomlin.  Scalable learning of safety guarantees for autonomous systems using hamilton-jacobi reachability.  In 2021 IEEE International Conference on Robotics and Automation (ICRA), pages 5914–5920. IEEE, 2021.
  * [15] Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell.  Simple and scalable predictive uncertainty estimation using deep ensembles.  Advances in neural information processing systems, 30, 2017.
  * [16] S. M. LaValle.  Planning Algorithms.  Cambridge University Press, Cambridge, U.K., 2006.  Available at http://planning.cs.uiuc.edu/.
  * [17] John Lygeros.  On reachability and minimum cost optimal control.  Automatica, 40(6):917–927, 2004.
  * [18] Ian Mitchell, Alex Bayen, and Claire J. Tomlin.  A time-dependent Hamilton-Jacobi formulation of reachable sets for continuous dynamic games.  IEEE Transactions on Automatic Control (TAC), 50(7):947–957, 2005.
  * [19] Luyao Niu, Hongchao Zhang, and Andrew Clark.  Safety-critical control synthesis for unknown sampled-data systems via control barrier functions.  In 2021 60th IEEE Conference on Decision and Control (CDC), pages 6806–6813. IEEE, 2021.
  * [20] Andrew Taylor, Andrew Singletary, Yisong Yue, and Aaron Ames.  Learning for safety-critical control with control barrier functions.  In Proceedings of the 2nd Conference on Learning for Dynamics and Control, volume 120 of Proceedings of Machine Learning Research, pages 708–717. PMLR, 10–11 Jun 2020.