††thanks: Corresponding author

# Spin foam amplitude of the black-to-white hole transition

Muxin Han<EMAIL_ADDRESS>Department of Physics, Florida Atlantic University,
777 Glades Road, Boca Raton, FL 33431-0991, USA  Department Physik, Institut
für Quantengravitation, Theoretische Physik III, Friedrich-Alexander
Universität Erlangen-Nürnberg, Staudtstr. 7/B2, 91058 Erlangen, Germany
Dongxue Qu<EMAIL_ADDRESS>Perimeter Institute for Theoretical
Physics, 31 Caroline St N, N2L 2Y5 Waterloo, ON, Canada  Cong Zhang
<EMAIL_ADDRESS>Department Physik, Institut für Quantengravitation,
Theoretische Physik III, Friedrich-Alexander Universität Erlangen-Nürnberg,
Staudtstr. 7/B2, 91058 Erlangen, Germany

###### Abstract

It has been conjectured that quantum gravity effects may cause the black-to-
white hole transition due to quantum tunneling. The transition amplitude of
this process is explored within the framework of the spin foam model on a
2-complex containing 56 vertices. We develop a systematic way to construct the
bulk triangulation from the boundary triangulation to obtain the 2-complex. By
using Thiemann’s complexifier coherent state as the boundary state to resemble
the semiclassical geometry, we introduce a procedure to calculate the
parameters labeling the coherent state from the continuous curved geometry.
Considering that triad fields of different orientations, i.e., $e_{i}^{a}$ and
$-e_{i}^{a}$, give the same intrinsic geometry of the boundary, we creatively
adopt the boundary state as a superposition of the coherent states associated
with both orientations. We employ the method of complex critical point to
numerically compute the transition amplitude. Despite the numerical results,
it is interestingly found that the transition amplitude is dominated by the
terms allowing the change in orientation. This suggests that the black-to-
white hole transition should be accompanied by quantum tunneling process of a
change in orientation.

##  I Introduction

The black hole (BH), one of the most mysterious objects in our Universe, has
drawn considerable attention in modern physics Maldacena (1996); Wald (2001);
Page (2005); Dittrich and Loll (2006); Unruh and Wald (2017); Akiyama  _et
al._ (2019); Gambini  _et al._ (2023); Lan  _et al._ (2023). Along the way of
our exploration on BHs, the Oppenheimer-Snyder model emerges as an important
influence in our understanding of the formation of BHs Oppenheimer and Snyder
(1939). While this model plays a crucial role in transitioning black holes
from abstract mathematical concepts to tangible physical realities, it faces
challenges, particularly concerning the presence of two singularities in the
spacetime it describes Joshi and Malafarina (2011). The first singularity
emerges at the final stage of the matter collapse and is referred to as the
cosmological singularity, while the second is the Schwarzschild singularity
existing in the vacuum region outside the collapsing dust. Indeed, the
existence of singularities is a common challenge within classical general
relativity (GR), indicating the invalidity of the classical theory in the
Planckian curvature region. This prompts the exploration of quantum gravity
(QG) to refine our understanding of gravity in these high-curvature regions
Polchinski (1998); Ashtekar and Lewandowski (2004); Han  _et al._ (2007);
Thiemann (2007); Ambjorn  _et al._ (2001); Rovelli and Vidotto (2014); Loll
(2020); Surya (2019); Gambini and Pullin (2013); Ashtekar  _et al._ (2018);
Song  _et al._ (2021); Zhang  _et al._ (2020a); Han and Liu (2022); Zhang  _et
al._ (2022a); Zhang (2021); Husain  _et al._ (2022); Lewandowski  _et al._
(2023); Gan  _et al._ (2022); Giesel  _et al._ (2023a, b).

Among various quantum Oppenheimer-Snyder models Tippett and Husain (2011);
Piechocki and Schmitz (2020); Ben Achour  _et al._ (2020); Lewandowski  _et
al._ (2023); Kiefer and Mohaddes (2023); Bobula and Pawłowski (2023), recent
attention has been directed toward the one that incorporates loop quantum
gravity (LQG) effects Lewandowski  _et al._ (2023). In this particular model,
one considers the collapse of a homogeneous and pressure less dust ball. The
collapse of the ball is proposed to follow the dynamics of loop quantum
cosmology, which predicts a bounce when the density of the collapsing ball
reaches a Planck scale, preventing the formation of the cosmological
singularity. This departure from the classical collapse will influence the
nature of the spacetime outside the dust ball. By assuming the spherical
symmetry and the static conditions, the metric in the outside region results
from the junction condition as follows:

$\displaystyle\differential s^{2}=-f(r)\differential
t^{2}+f(r)^{-1}\differential r^{2}+r^{2}\differential\Omega^{2},$ (1)

where $f(r)$ is given by

$f(r)=1-\frac{2GM}{r}+\frac{\alpha G^{2}M^{2}}{r^{4}},\
\alpha=2\sqrt{3}\beta^{3}\kappa\hbar,$ (2)

with the Barbero-Immirzi parameter $\beta$ and $\kappa=8\pi G$ (see e.g.,
Marto  _et al._ (2015); Parvizi  _et al._ (2022); Husain  _et al._ (2022);
Bobula and Pawłowski (2023); Giesel  _et al._ (2023c); Lin and Zhang (2024)
for similar works). The Penrose diagrams of the spacetime with the metric (1)
for different value of $M$ are discussed in Lewandowski  _et al._ (2023). It
turns out that the diagrams share similar properties to those of the
Reissner–Nordström BH.

While the spacetime given by (1) offers distinct advantages, based on which
observational effects have been studied Yang  _et al._ (2023); Gingrich
(2024); Zhang  _et al._ (2023); Shao  _et al._ (2023); Ye  _et al._ (2023);
Chen  _et al._ (2024), it is not without debates. For instance, the existence
of the Cauchy horizon implies that the spacetime could be unstable under
perturbation Cao  _et al._ (2023, 2024). Moreover, as in the classical
Oppenheimer-Snyder model, this spacetime does not incorporate Hawking
evaporation. These unresolved issues highlight the ongoing complexity of our
understanding and indicate that we have not reached a final conclusion in
exploring the quantum Oppenheimer-Snyder model. An aspect that we may have
overlooked is what happens at the end of the BH evaporation. As discussed in
Haggard and Rovelli (2015); D’Ambrosio  _et al._ (2021); Rignon-Bret and
Rovelli (2022); Soltani  _et al._ (2021); Han  _et al._ (2023a), the
consideration of the Hawking evaporation reveals that the horizon will
eventually shrink to Planckian size. Under such conditions, QG effects may
play a significant role, potentially leading to quantum tunneling from a BH to
a white hole (WH). Our work will delve into the study of this phenomenon
within the framework of the spin foam model in LQG Engle  _et al._ (2007,
2008); Rovelli and Vidotto (2014).

As proposed in Han  _et al._ (2023a), the occurrence of quantum tunneling will
alter the global structure of the spacetime. The revised spacetime exhibits a
global structure similar to that of the Minkowski spacetime. The metric tensor
is locally the same as given in (1), except for a region, called the $B$
region, where quantum tunneling takes place. The embedding of the $B$ region
in the entire spacetime will equip this region with the boundary data, i.e.,
the reduced $3$-metric and the extrinsic curvature on the boundary. Our task
is to study the spin foam amplitude of the $B$ region with this boundary data.

Several challenges prevent the direct application of the spin foam model to
the $B$ region embedded in classical spacetime. The first issue is that the
spin foam amplitude is defined based on discrete graphs, which are usually
dual to a triangulation. This thus needs us to triangulate the $B$ region at
first. We will construct a triangulation adapting to the symmetries of the $B$
region, namely its spherical symmetry and time reversal symmetry. This
triangulation results in a dual 2-complex, denoted by ${\mathcal{K}}$. The
intersection of ${\mathcal{K}}$ with the boundary $\partial B$ of the $B$
region yields a graph $\gamma$, on which one can define the boundary Hilbert
space $\mathcal{H}_{\gamma}$ of spin network states Ashtekar and Lewandowski
(2004); Thiemann (2007). The spin foam amplitude $A$ associated with the
2-complex is a specific functional on $\mathcal{H}_{\gamma}$, namely
$A:\mathcal{H}_{\gamma}\to\mathbb{C}$ Rovelli and Vidotto (2014). To make the
Hilbert space $\mathcal{H}_{\gamma}$ well-defined, it is required that the
boundary $\partial B$ should be spacelike. However, the initial definition of
$\partial B$ in Han  _et al._ (2023a) is null. To resolve this issue, we need
to choose a slightly larger region $B^{\prime}$ in the spacetime, ensuring
that $B^{\prime}$ contains $B$ and that the boundary $\partial B^{\prime}$ is
spacelike. Indeed, all our calculations, including the triangulation
aforementioned, are performed in $B^{\prime}$ rather than $B$. Throughout the
subsequent discussion, we will continue to use $B$ to denote this slightly
larger region $B^{\prime}$ for consistency. Another noteworthy concern is the
mismatch between the amplitude $A$, defined on the Hilbert space
$\mathcal{H}_{\gamma}$, and the classical geometry specifying the boundary
condition on $\partial B$. This requires us to choose a coherent state in
$\mathcal{H}_{\gamma}$ to resemble the classical geometry. Then, the value of
$A$ on this coherent state is interpreted as the probability amplitudes which
describes the transition, occurring in the $B$ region, between the classical
boundary geometries.

The investigation of the dynamics within the $B$ region through the spin foam
amplitude implies the dominance of QG effects in this region. This assigns a
crucial physical interpretation to the size of the $B$ region, indicating when
the influence of QG effects cannot be neglected. Importantly, within a well-
defined dynamics the moment when the significance of QG effects emerges should
be determined by the dynamics itself. From this point of view, the size of the
$B$ region should be viewed as a dynamical parameter. More precisely, we need
to calculate the amplitude $A$ for a family of $B$ regions of various sizes,
and pick out the one at which $A$ takes its extreme value. Moreover, as
discussed in Han  _et al._ (2023a), the size of the $B$ region is also
relevant to other meaningful observables. Therefore, treating the size as a
dynamical parameter can also allow us to obtain the value of these observables
dynamically. This discussion clarifies the predictability of the spin foam
model.

The paper is constructed as follows: In Sec. II, the spin foam model in the
$B$ region is constructed, where we introduce the systematical method to yield
the bulk triangulation from the boundary triangulation, define the boundary
coherent state, and introduce the spin foam amplitude. Sec. III describes how
to calculate the boundary data of the boundary state from the curved geometry
outside the $B$ region. In Sec. IV, the algorithm for computing the complex
critical point is introduced. In Sec. V, we present the numerical results.
Finally, this work is summarized in Sec. VI. To keep the flow of the
discussion, the theoretical derivations forming the foundation of this work
are presented in the appendices.

##  II Spinfoam model in the ${B}$ region

The spacetime structure proposed in Han  _et al._ (2023a) is presented in Fig.
1. The $B$ region is the gray shaded region in the spacetime. The metric
outside the region is the same as the one in (1). In the $B$ region, quantum
tunneling of black-to-white hole transition occurs, and the spin foam model is
used to describe the dynamics. This section will focus on the $B$ region to
construct the spin foam model.

As mentioned in Sec. I, the $B$ region has the time reversal symmetry. The
time reversal symmetry is with respect to a hypersurface $\mathcal{T}$ within
$B$. This hypersurface $\mathcal{T}$, referred to as the transition surface,
divides $B$ into a union of the past region $B_{-}$ and the future region
$B_{+}$. Let $\partial_{-}B$ and $\partial_{+}B$ be the remaining boundary of
$B_{-}$ and $B_{+}$ except $\mathcal{T}$, so that $\partial
B=\partial_{-}B\cup\partial_{+}B$. The time reversal symmetry gives a
diffeomorphism between $B_{+}$ and $B_{-}$ that preserves the transition
surface $\mathcal{T}$ and maps between the boundary conditions on
$\partial_{+}B$ and $\partial_{-}B$ (see Fig. 1).

The first step to construct the model is to triangulate the $B$ region, so
that the spin foam amplitude can be defined on the 2-complex dual to the
triangulation. In our work, we only consider the triangulation invariant under
the time reversal diffeomorphism. To this end, we firstly triangulate the
region $B_{-}$, and then apply the diffeomorphism to transform the
triangulation from $B_{-}$ to $B_{+}$.

Figure 1: The BH spacetime with the $B$ region. The gray shaded region
represents the $B$ region with the boundaries $\partial
B=\partial_{+}B\cup\partial_{-}B$. For our triangulation, we divide the
boundary $\partial_{-}B$ into a union of two parts,
$\partial_{-}B_{1}\cup\partial_{-}B_{2}$.

###  II.1 Three basic triangulations

For clarity, we use the word ‘point’ to refer to the vertex of a triangulation
in this paper.

To understand the triangulation structure in this work, we begin with three
simple examples closely related to the main problem. These examples help
illustrate how to construct the bulk triangulation from the boundary
triangulation. Let us consider a spacetime region denoted as $R$, bounded by
two spatial hypersurfaces $P$ and $P^{\prime}$, as depicted in Fig. 2. The
intersection of $P$ and $P^{\prime}$ forms a 2-surface $S$.

In the first case, we consider a triangulation of the spatial surface $P$
using the tetrahedra $ABCD$, where the triangle $ABC$ lies in the 2-plane $S$,
as depicted in Fig. 2. We call the point $D$ as the dynamical point, which
moves from $P$ to $P^{\prime}$ and eventually yield a new point $D^{\prime}$
in $P^{\prime}$. By connecting the point $D^{\prime}$ with the points $A$,
$B$, and $C$ respectively, we obtain the triangulation $ABCDD^{\prime}$ of the
region $R$. It is worth noting that the tetrahedron $ABCD^{\prime}$ gives the
triangulation of the hypersurface $P^{\prime}$. We call this tetrahedron the
final tetrahedron. Besides the final tetrahedron $ABCD^{\prime}$, the
4-simplex contains four other boundary tetrahedra: the initial tetrahedron
$ABCD$ and the tetrahedra $ACDD^{\prime}$, $ABDD^{\prime}$, and
$BCDD^{\prime}$. The last three tetrahedra are viewed as the world volume of
the triangles $ACD$, $ABD$ and $BCD$ respectively.

Figure 2: Triangulation of the region $R$. The red line $DD^{\prime}$ is the
wordline of the vertex $D$. The triangle $ABC$ lies in the 2-surface $S$. The
initial tetrahedron $ABCD$ triangulates the hypersurface $P$, while the final
tetrahedron $ABCD^{\prime}$ results in a triangulation of the hypersurface
$P^{\prime}$.

In a more complicated case, we triangulate the spatial surface $P$ using the
tetrahedron $ABCD$ where only the segment $AB$ lies in the 2-plane $S$. In
this case, we have two dynamical points $C$ and $D$ as shown in Fig. 3. Both
$C$ and $D$ travel from $P$ to $P^{\prime}$, resulting in the points
$C^{\prime}$ and $D^{\prime}$. To establish the corresponding
triangularization, we firstly examine the evolution of the tetrahedron $ABCD$
with $D$ treated as the dynamical point. This evolution yields a 4-simplex
$ABCDD^{\prime}$, as what we did in the previous case. Within this 4-simplex,
the final tetrahedron $ABCD^{\prime}$ represents the world volume of the
triangle $ABC$. Subsequently, we explore the evolution of the tetrahedron
$ABCD^{\prime}$ with $C$ as the dynamical point, generating a new 4-simplex
$ABCD^{\prime}C^{\prime}$. The cell $ABCDC^{\prime}D^{\prime}$ is then
triangulated by the 4-simplices $ABCDD^{\prime}$ and
$ABCD^{\prime}C^{\prime}$.

Figure 3: The triangulation of the cell $ABCDD^{\prime}C^{\prime}$. The right
side shows the 4-simplices triangulating the cell. The first 4-simplex
$ABCDD^{\prime}$ results from the evolution of the tetrahedron $ABCD$ with $D$
as the dynamical point. The second 4-simplex $ABCC^{\prime}D^{\prime}$ is
generated by the evolution of the tetrahedron $ABCD^{\prime}$ with $C$ treated
as the dynamical point.

In the above two cases, the boundary triangulation for $P$ remains the same
while the number of dynamical points varies. In the subsequent case, we
introduce a more complicated triangulation for $P$ to demonstrate how the
procedure outlined in the previous cases can be applied. Specifically, we
triangulate the spatial surface $P$ using two tetrahedra $ABCD$ and $BCDE$,
where the segment $AB$ lies in the 2-sphere $S$. These tetrahedra share the
triangle $BCD$ (Fig. 4). Treating the points $C$ and $D$ as the dynamical
ones, we obtain the resulting cell $ABCDEC^{\prime}D^{\prime}$, as illustrated
in Fig. 4. Since the cell $ABCDEC^{\prime}D^{\prime}$ contains of two subcells
$ABCDC^{\prime}D^{\prime}$ and $BCDEC^{\prime}D^{\prime}$, the triangulation
of $ABCDEC^{\prime}D^{\prime}$ is formed by combining the triangulations of
these subcells. The triangulation for the subcells can be performed by
following the procedure depicted in Figure 3. As a consequence, the cell
$ABCDEC^{\prime}D^{\prime}$ emerges by gluing together the four 4-simplices
$ABCDD^{\prime}$, $ABCD^{\prime}C^{\prime}$, $BCDED^{\prime}$, and
$BCED^{\prime}C^{\prime}$, as illustrated in Figure 4.

Figure 4: The triangulation of the cell $ABCDED^{\prime}C^{\prime}$ (the
left), where $ABCDED^{\prime}C^{\prime}$ is triangulated into the four
4-simplices $ABCDD^{\prime}$, $ABCD^{\prime}C^{\prime}$, $BCDED^{\prime}$, and
$BCED^{\prime}C^{\prime}$ (the right).

###  II.2 Triangulation of the region $B_{-}$

Let us divide the past boundary $\partial_{-}B$ of $B$ into a union of two
parts $\partial_{-}B_{1}$ and $\partial_{-}B_{2}$. The intersections between
any two among the hypersurfaces $\partial_{-}B_{1}$, $\partial_{-}B_{2}$ and
$\mathcal{T}$ are 2-spheres, denoted by $S_{1}$, $S_{2}$ and $S_{3}$
respectively (refer to Fig. 1). Following the procedures outlined in Sec.
II.1, the triangulation of $B_{-}$ is uniquely determined, if we perform a
specific triangulation for $\partial_{-}B_{1}$ where the points of the
triangulation lie either on $S_{1}$ or $S_{2}$, and the point on $S_{2}$ are
chosen as the dynamical points.

To construct the triangulation of $\partial_{-}B_{1}$, let us embed it into a
$\mathbb{R}^{3}\ni\vec{x}$ as a spherical shell, i.e., a 3D annulus, which is
topologically homomorphism to $\partial_{-}B_{1}$. For convenience, the center
of the spherical shell is put at the origin. The radii of the spheres $S_{1}$
and $S_{2}$ are $r_{1}$ and $r_{2}$ respectively, with $r_{1}<r_{2}$. We adopt
the triangulation of $\partial_{-}B_{1}$ proposed in D’Ambrosio  _et al._
(2021); Soltani  _et al._ (2021), where one introduces a tetrahedron with
points $\mathfrak{p}_{a}$ ($a=1,2,3,4$) to approximate the sphere $S_{1}$.
Even though the shape of the tetrahedron will not influence the structure of
the triangulation, we would choose it to be regular for convenience. Let
$r_{1}\vec{x}(a)\in\mathbb{R}^{3}$ for $a=1,2,3,4$ to be the coordinate of
$\mathfrak{p}_{a}$. Then, select the points $\mathfrak{p}_{a^{\prime}}\in
S_{2}$ taking coordinates $-r_{2}\vec{x}(a)$. The connections between these
points $\mathfrak{p}_{a^{\prime}}$ ($a=1,2,3,4$) give rise to another regular
tetrahedron that approximates the sphere $S_{2}$. Finally, the triangulation
of $\partial_{-}B_{1}$ is obtained by drawing line segments connecting each
point $\mathfrak{p}_{a}$ (for $a=1,2,3,4$) in $S_{1}$ to the corresponding
points $\mathfrak{p}_{b^{\prime}}$, $\mathfrak{p}_{c^{\prime}}$, and
$\mathfrak{p}_{d^{\prime}}$ in $S_{2}$, where $b$, $c$, and $d$ represent
labels that are different from $a$. Importantly, the minus sign in the
coordinate of $\mathfrak{p}_{a^{\prime}}$ ensures that this connectivity
structure yields a well-defined triangulation. The resulting triangulation can
be visualized as shown in Fig. 5. As illustrated in the figure, there are a
total of $14$ tetrahedra, classified into three types:

  * •

A-type: There are 4 tetrahedra having one point on $S_{1}$ and three points on
$S_{2}$. The points of the tetrahedron are
$\mathfrak{p}_{a}\mathfrak{p}_{b^{\prime}}\mathfrak{p}_{c^{\prime}}\mathfrak{p}_{d^{\prime}}$.

  * •

B-type: There are 6 tetrahedra having two points in $S_{1}$ and two points in
$S_{2}$. The points of the tetrahedron are
$\mathfrak{p}_{a}\mathfrak{p}_{b}\mathfrak{p}_{c^{\prime}}\mathfrak{p}_{d^{\prime}}$.

  * •

C-type: There are 4 tetrahedra having three points on $S_{1}$ and one point on
$S_{2}$. The points are
$\mathfrak{p}_{a}\mathfrak{p}_{b}\mathfrak{p}_{c}\mathfrak{p}_{d^{\prime}}$.

Here, $a,b,c,d$ represent distinct elements of the set $\\{1,2,3,4\\}$.

Figure 5: A visualization for the triangulation of $\partial_{-}B_{1}$. To
enhance clarity, both the transparent (left) and untransparent (right)
versions are provided.

To obtain a triangulation for the 4-dimensional region $B_{-}$ from the
triangulation of its boundary $\partial_{-}B_{1}$, we selected the points in
$S_{2}$ to be dynamical. Starting from each of these dynamical points
$\mathfrak{p}_{a^{\prime}}$ for $a=1,2,3,4$, a line
$\mathfrak{p}_{a^{\prime}}\mathfrak{p}_{a^{\prime\prime}}$ is drawn in
$\partial_{-}B_{2}$, intersecting the sphere $S_{3}$ at a chosen point
$\mathfrak{p}_{a^{\prime\prime}}$.
$\mathfrak{p}_{a^{\prime}}\mathfrak{p}_{a^{\prime\prime}}$ is viewed as the
worldline of $\mathfrak{p}_{a^{\prime}}$ in our final triangulation for
$B_{-}$. Following the method outlined in Sec. II.1, one can construct the
triangulation of $B_{-}$. Specially, each A-type tetrahedron contains three
dynamical points, whose evolution results in a total of 12 4-simplices. Each
B-type tetrahedron has two dynamical points, whose evolution leads to a total
of 12 4-simplices. Each C-type tetrahedron has one dynamical point, whose
evolution results in a total of 4 4-simplices. Consequently, there are a total
of $28$ 4-simplices in the triangulation for the region $B_{-}$.

Consider embedding $S_{3}$ into $\mathbb{R}^{3}$ as a sphere centered at the
origin with a radius $r_{3}>r_{2}$. The spherical shell between the two
spheres $S_{2}$ and $S_{3}$ is topologically the same as $\partial_{-}B_{2}$.
Then, $\mathfrak{p}_{a^{\prime\prime}}$ for $a=1,2,3,4$ are chosen so that
their coordinates are $r_{3}\vec{x}(a)\in\mathbb{R}^{3}$. The resulting
triangulation for the transition surface $\mathcal{T}$, as the boundary of the
triangulation of $B_{-}$, is the same as that for $\partial_{-}B_{1}$, i.e.,
the same as shown in Fig. 5.

The 2-complex $\mathcal{K}$ dual to the triangulation of $B_{-}$ is shown in
Fig. 6. This dual complex contains 28 vertices, 90 edges, and 104 faces, which
correspond to the 4-simplices, tetrahedra, and triangles in the original
triangulation, respectively. The structure of the 2-complex is encoded in the
following function, also denoted as $\mathcal{K}$,

$\mathcal{K}:f\mapsto\\{(v_{1},e_{1}),(v_{2},e_{2}),\cdots,(v_{n},e_{n})\\},$
(3)

where the domain of $\mathcal{K}$ consists of all faces $f$, and the sequence
in the right-hand-side contains all half edges $(v,e)$ surrounding $f$.
Importantly, the order of the sequence coincides with the orientation of the
face dual to $f$. Specific details regarding the assignment of the map
$\mathcal{K}$ are provided in Zhang (2024).

Taking advantage of the 2-complex for $B_{-}$ region, we construct the
2-complex for the entire $B$ region easily due to the time reversal symmetry.

Figure 6: 2-complex dual to the triangulation of the region $B_{-}$.
$v_{abcde}$ denotes the vertex dual to the 4-simplex
$(\mathfrak{p}_{a}\mathfrak{p}_{b}\mathfrak{p}_{c}\mathfrak{p}_{d}\mathfrak{p}_{e})$
and $e_{abcd}$ is the edge dual to the tetrahedron
$(\mathfrak{p}_{a}\mathfrak{p}_{b}\mathfrak{p}_{c}\mathfrak{p}_{d})$. In this
figure, we use the black solid lines to represent the edges passing through
the transition surface $\mathcal{T}$, and the blue dashed line for the edges
in $\partial_{-}B$.

###  II.3 Boundary state

The intersection between $\mathcal{K}$ and $\partial B$ gives the boundary
graph $\gamma$. $\gamma$ contains the nodes $n_{e}$ and links $\ell_{f}$,
which result from the intersections of edges $e$ and faces $f$, respectively,
with $\partial B$. All nodes of $\gamma$ are 4-valence. The boundary states
are vectors in the LQG Hilbert space $\mathcal{H}_{\gamma}$ on the graph
$\gamma$. We will call $\mathcal{H}_{\gamma}$ the boundary Hilbert space in
the current work. Specifically, $\mathcal{H}_{\gamma}$ is given by

$\mathcal{H}_{\gamma}=\bigotimes_{\ell\in\gamma}\left(\bigoplus_{j_{\ell}}\mathcal{H}_{j_{\ell}}\otimes\mathcal{H}_{j_{\ell}}^{*}\right)$
(4)

where $\mathcal{H}_{j}$ is the $\mathrm{SU(2)}$ representation space with spin
$j$ and $\mathcal{H}_{j}^{*}$ is its dual space. The spaces
$\mathcal{H}_{j_{\ell}}$ and $\mathcal{H}_{j_{\ell}}^{*}$ are associated with
the starting node $s_{\ell}$ and targeting node $t_{\ell}$ of $\ell$
respectively.

As discussed in Sec. I, the current work uses the coherent states as the
boundary states of spinfoam amplitude. A class of coherent states that has
been well-understood is the Thiemann’s complexifier coherent state Thiemann
(2001); Sahlmann  _et al._ (2001); Thiemann (2006) (see, e.g. Bianchi  _et
al._ (2010a, b) for the earlier applications to the spin foam model). Given a
link $\ell\equiv\ell_{f}$, a coherent state labelled by
$g_{\ell}\in\mathrm{SL}(2,\mathbb{C})$ can be expressed as a function of the
holonomy $h_{\ell}$ as

$\psi_{g_{\ell}}^{t}(h_{\ell})=\mathcal{N}\sum_{j_{\ell}}d_{j_{\ell}}e^{-\frac{t}{2}j_{\ell}(j_{\ell}+1)}\sum_{m}D^{j_{\ell}}_{mm}(g_{\ell}^{-1}h_{\ell}),$
(5)

with $\mathcal{N}$ denoting the normalization factor. It is convenient to
decompose $g_{\ell}$ as

$g_{\ell}=n_{\ell s_{\ell}}e^{(-ip_{\ell}+\phi_{\ell})\tau_{3}}n_{\ell
t_{\ell}}^{-1}$ (6)

with $n_{\ell s_{\ell}},\ n_{\ell t_{\ell}}\in\mathrm{SU(2)}$ belonging to the
starting and target points of $\ell$ and $p_{\ell},\
\phi_{\ell}\in\mathbb{R}$. The geometry interpretation for the parameters
$n_{\ell s_{\ell}}$, $n_{\ell t_{\ell}}$, $p_{\ell}$ and $\phi_{\ell}$ can be
deduced from the expectation values of operators (refer to Appendix C.3 and
Thiemann and Winkler (2001); Dapor and Liegener (2018); Zhang  _et al._
(2022b)). We have

$\displaystyle\langle\psi_{g_{\ell}}^{t}|\hat{p}_{s,\ell}^{k}\tau_{k}|\psi_{g_{\ell}}^{t}\rangle$
$\displaystyle\cong\frac{1}{t}p_{\ell}n_{\ell s_{\ell}}\tau_{3}n_{\ell
s_{\ell}}^{-1},$ (7)
$\displaystyle\langle\psi_{g_{\ell}}^{t}|\hat{p}_{t,\ell}^{k}\tau_{k}|\psi_{g_{\ell}}^{t}\rangle$
$\displaystyle\cong-\frac{1}{t}p_{\ell}n_{\ell t_{\ell}}\tau_{3}n_{\ell
t_{\ell}}^{-1},$
$\displaystyle\langle\psi_{g_{\ell}}^{t}|\hat{h}_{\ell}|\psi_{g_{\ell}}^{t}\rangle$
$\displaystyle\cong n_{\ell s_{\ell}}e^{\phi_{\ell}\tau_{3}}n_{\ell
t_{\ell}}^{-1},$

where $\hat{p}_{s,\ell}^{k}$ and $\hat{p}_{t,\ell}^{k}$ are the flux operators
and $\hat{h}_{\ell}$ denotes the holonomy operator (refer to Appendix C.2).

In Appendix B, the explicit expression of the spin foam amplitude with the
boundary state (5) is derived. Since the current work concerns with the
asymptotics as $t\to 0$, we only need to focus on the critical equation. Then,
as shown in Appendix B, the same critical equation can be obtained if we
replace the Thiemann’s coherent state by

$\displaystyle|g_{\ell}\rangle=$
$\displaystyle\mathcal{N}\sum_{j_{f}}\exp(-\frac{t}{2}\left(j_{f}-\frac{|p_{\ell}|}{t}\right)^{2}+ij_{f}\mathrm{sgn}(p_{\ell})\phi_{\ell})\times$
(8) $\displaystyle|j_{f}\xi_{\ell t_{\ell}}\rangle\otimes\langle
j_{f}\xi_{\ell s_{\ell}}|,$

where $\xi_{x}$ for $x=(\ell s_{\ell})$ or $(\ell t_{\ell})$ is the spinor
given by

$\xi_{x}=\left\\{\begin{array}[]{cc}n_{x}(1,0)^{T},&p_{\ell}>0,\\\
n_{x}(0,1)^{T},&p_{\ell}<0,\end{array}\right.$ (9)

and $|j_{\ell},\xi_{\ell t_{\ell}}\rangle\in\mathcal{H}_{j_{\ell}}$ and
$\langle j_{\ell},\xi_{\ell s_{\ell}}|\in\mathcal{H}_{j_{f}}^{*}$ are the
$\mathrm{SU(2)}$ coherent states Klauder and Skagerstam (1985); Alesci  _et
al._ (2016). The state (8) is the similar coherent state used in the usual
context of the spin foam model Bianchi  _et al._ (2010a), and its
generalization to higher dimensional LQG is also studied Long  _et al._ (2021,
2022). In terms of the state (8), (7) becomes

$\displaystyle\langle g_{\ell}|\hat{p}_{s,\ell}^{k}|g_{\ell}\rangle\cong$
$\displaystyle\frac{1}{t}|p_{\ell}|\xi_{\ell
s_{\ell}}^{\dagger}\sigma^{k}\xi_{\ell s_{\ell}},$ (10) $\displaystyle\langle
g_{\ell}|\hat{p}_{t,\ell}^{k}|g_{\ell}\rangle\cong$
$\displaystyle-\frac{1}{t}|p_{\ell}|\xi_{\ell
t_{\ell}}^{\dagger}\sigma^{k}\xi_{\ell t_{\ell}},$ $\displaystyle\langle
g_{\ell}|\hat{h}_{\ell}|g_{\ell}\rangle\cong$ $\displaystyle n(\xi_{\ell
s_{\ell}})e^{\mathrm{sgn}(p_{\ell})\phi_{\ell}\tau_{3}}n(\xi_{\ell
t_{\ell}})^{-1},$

where $n(\xi)$ for $\xi=(\xi_{1},\xi_{2})^{T}\in\mathbb{C}^{2}$ is given by

$n(\xi)=\frac{1}{\sqrt{\xi^{\dagger}\xi}}\begin{pmatrix}\xi_{1}&-\overline{\xi_{2}}\\\
\xi_{2}&\overline{\xi_{1}}\end{pmatrix}\in\mathrm{SU(2)}.$ (11)

In terms of these parameters, $g_{\ell}$ is expressed as

$g_{\ell}=n(\xi_{\ell
s_{\ell}})e^{(-i|p_{\ell}|+\mathrm{sgn}(p_{\ell})\phi_{\ell})\tau_{3}}n(\xi_{\ell
t_{\ell}})^{-1}.$ (12)

We get the coherent state on the entire graph $\gamma$ as their tensor
product,

$\left|\\{g_{\ell}\\}_{\ell\subset\gamma}\right\rangle=\bigotimes_{\ell\in\gamma}|g_{\ell}\rangle.$
(13)

Even though the boundary state given in (13) lacks gauge invariance, the spin
foam amplitude will inherently projects it into the gauge invariant subspace.
As a consequence, the boundary state utilized in the spin foam amplitude is in
essence the projected gauge invariant the state.

According to our discussion in Appendix C, the parameters $p_{\ell}$ and
$\xi_{\ell}$ in the coherent state carry the information about the densitied
triad $E^{a}_{i}$ on $\partial B$, while $\phi_{\ell}$ is relevant to the
Ashtekar connection $A_{a}^{i}$ on $\partial B$ (see also e.g. Rovelli and
Speziale (2010)). Therefore, the coherent state carries the information about
the 4-dimensional geometry of $\partial B$, relating the boundary geometry to
the boundary coherent state labels.

The intrinsic geometry of the boundary is characterized by the 3-metric
$q_{ab}$, which is related to the triad components through
$q_{ab}=\delta_{ij}e^{i}_{a}e^{j}_{b}$. This relationship suggests an $O(3)$
gauge transformation, wherein the triad undergoes a transformation
$e^{i}_{a}\mapsto\mathcal{O}^{i}{}_{j}e^{j}_{a}$ with $\mathcal{O}\in O(3)$.
However, in LQG, only $\mathrm{SO(3)}$ gauge transformations are considered.
Consequently, the triads $e^{i}_{a}$ and $-e^{i}_{a}$, while yielding the same
intrinsic geometry, are not considered gauge-equivalent in the context of LQG.
Based on this discussion, it seems reasonable to consider a superposition of
states associated with both triads $e^{i}_{a}$ and $-e^{i}_{a}$ to capture a
more comprehensive description of the boundary geometry.

In the current work, the boundary $\partial B$ contains $\partial_{-}B$ and
$\partial_{+}B$. Let $\tilde{e}^{i}_{a}$ represent a triad in $\partial_{+}B$,
and $\bar{e}^{i}_{a}$ represent a triad in $\partial_{-}B$. Then, across the
entire boundary $\partial B$, the triad fields can take the forms
$(\tilde{e}_{a}^{i},\bar{e}^{i}_{a})$, $(\tilde{e}_{a}^{i},-\bar{e}^{i}_{a})$,
$(-\tilde{e}_{a}^{i},\bar{e}^{i}_{a})$, and
$(-\tilde{e}_{a}^{i},-\bar{e}^{i}_{a})$. Here, we employ pairs like
$(e_{+},e_{-})$ to denote the piecewise triad field, where $e(x)=e_{+}(x)$ for
$x\in\partial_{+}B$ and $e(x)=e_{-}(x)$ for $x\in\partial_{-}B$. Let us use
$g_{\ell}^{(\tilde{s}\bar{s})}$, for each link $\ell$, to denote the label of
the coherent state constructed from the triad
$(\tilde{s}\tilde{e}_{a}^{i},\bar{s}\bar{e}^{i}_{a})$ with
$\tilde{s},\bar{s}=\pm$. To keep the discussion flow, we will momentarily skip
detailing the explicit relation between $g_{\ell}^{(\tilde{s}\bar{s})}$ and
$(\tilde{s}\tilde{e}_{a}^{i},\bar{s}\bar{e}^{i}_{a})$. The above discussion
motivates us to consider the boundary state as the superposition

$|\Psi_{\gamma}\rangle=\sum_{\tilde{s},\bar{s}=\pm}\left|\\{g_{\ell}^{(\tilde{s}\bar{s})}\\}_{\ell\subset\gamma}\right\rangle.$
(14)

Then the spin foam transition amplitude is a sum

$\langle A|\Psi_{\gamma}\rangle=\sum_{\tilde{s},\bar{s}=\pm}\left\langle
A\middle|\\{g_{\ell}^{(\tilde{s}\bar{s})}\\}_{\ell\subset\gamma}\right\rangle.$
(15)

###  II.4 Spin foam amplitude and action

The spin foam amplitude is a specific linear functional acting in the boundary
Hilbert space (4) Rovelli and Vidotto (2014). The spinfoam amplitude with the
boundary state (13) has an integral expression. This expression assigns an
$\mathrm{SU(2)}$ representation $j_{f}$ to each face $f$, an
$\mathrm{SL}(2,\mathbb{C})$ element $g_{ve}$ to each half edge $(v,e)$ of the
edge $e$ taking the vertex $v$ as its end point, and a spinor
$z_{vf}\in\mathbb{C}^{2}$ to each corner $(v,f)$ of the face $f$ at the vertex
$v$. The amplitude reads

$\displaystyle A=$
$\displaystyle\sum_{\\{j_{f}\\}}\int\prod_{f}(d_{j_{f}})^{|V_{f}|+1}\prod_{ve}\differential
g_{ve}\prod_{vf}\differential^{2}\Omega_{vf}\times$ (16)
$\displaystyle\exp[S(\\{j_{f}\\},\\{g_{ve}\\},\\{z_{vf}\\})].$

where $|V_{f}|$ is the number of vertices surrounding the face $f$,
$\differential g_{ve}$ is the Haar measure, $\differential^{2}\Omega_{vf}$ is
a measure on $\mathbb{C}\mathbb{P}^{1}$, and $S$ depending on the tuples
$(\\{j_{f}\\},\\{g_{ve}\\},\\{z_{vf}\\})$ is the effective action given by the
spin foam model. The precise definition of the measures $\differential g_{ve}$
and $\differential^{2}\Omega_{vf}$ can be found in Han and Krajewski (2014).
The integration for the spinor $z_{vf}\in\mathbb{C}^{2}$ is over
$\mathbb{CP}^{1}$, relying on the scaling invariance of the action $S$. The
action $S$ can be written as the following summation

$S(\\{j_{f}\\},\\{g_{ve}\\},\\{z_{vf}\\})=\sum_{f}S_{f}(j_{f},g_{ve},z_{vf}),$
(17)

where $S_{f}$ is the action associated with the face $f$. Specifically, for an
internal face $f$, $S_{f}$ reads

$\displaystyle S_{f}=$ $\displaystyle\sum_{e\in\partial
f}j_{f}\Bigg{[}2\ln\frac{\langle
Z_{s_{e}ef},Z_{t_{e}ef}\rangle}{\left\|Z_{s_{e}ef}\right\|\left\|Z_{t_{e}ef}\right\|}+i\beta\ln\frac{\left\|Z_{s_{e}ef}\right\|^{2}}{\left\|Z_{t_{e}ef}\right\|^{2}}\Bigg{]},$
(18)

where $Z_{vef}:=g_{ve}^{\dagger}z_{vf}$, $s_{e}$ and $t_{e}$ represent the
source and target points of edge $e$ respectively, and $\|Z\|$ denotes
$\sqrt{\langle Z,Z\rangle}$ for all $Z\in\mathbb{C}^{2}$. For a boundary face
$f$, as shown in (135) (see Appendixes A and B for more details), we have

$\displaystyle S_{f}=$
$\displaystyle-\frac{t}{2}\left(j_{f}-\frac{|p_{\ell_{f}}|}{t}\right)^{2}+i\,\mathrm{sgn}(p_{\ell_{f}})j_{f}\phi_{\ell_{f}}+$
(19) $\displaystyle\sum_{e\in\partial f}j_{f}\Bigg{[}2\ln\frac{\langle
Z_{s_{e}ef},Z_{t_{e}ef}\rangle}{\left\|Z_{s_{e}ef}\right\|\left\|Z_{t_{e}ef}\right\|}+i\beta\ln\frac{\left\|Z_{s_{e}ef}\right\|^{2}}{\left\|Z_{t_{e}ef}\right\|^{2}}\Bigg{]}$

where the extra term compared to (18) comes from the coherent state. Moreover,
in (19), $Z_{s_{e}ef}=\xi_{\ell_{f}s_{\ell_{f}}}$ when the boundary of the
face $f$ starts from $e$, and $Z_{t_{e}ef}=\xi_{\ell_{f}t_{\ell_{f}}}$ when
the boundary of the face $f$ ends at $e$, with $\xi_{\ell_{f}s_{\ell}}$ and
$\xi_{\ell_{f}t_{\ell}}$ introduced in (8).

Let $F$, $E$, and $V$ denote the number of faces, _half_ edges, and vertices
in the dual 2-complex $\mathcal{K}$. The action $S$ is a function on the
configuration space

$\Gamma=\left(\mathbb{Z}_{\geq
0}/2\right)^{F}\times\mathrm{SL}(2,\mathbb{C})^{E}\times(\mathbb{C}^{2})^{10V},$

with $j_{f}\in\mathbb{Z}_{\geq 0}/2$, $g_{ve}\in\mathrm{SL}(2,\mathbb{C})$ and
$z_{vf}\in\mathbb{C}^{2}$.

Considering arbitrary $\mathrm{SL}(2,\mathbb{C})$-valued field $v\mapsto
h_{v}$ at the vertices, arbitrary complex field $(v,f)\mapsto\lambda_{vf}$
with $\lambda_{vf}\neq 0$, and arbitrary $\mathrm{SU(2)}$-valued field
$e\mapsto u_{e}$ on the edges, we observe that $S$ exhibits invariance under
the transformation:

$\Big{(}j_{f},g_{ve},z_{vf}\Big{)}\mapsto\Big{(}j_{f},h_{v}^{-1}g_{ve}u_{e},\lambda_{vf}h_{v}^{\dagger}z_{vf}\Big{)}.$
(20)

The integration in (16) is carried out over the quotient space $\Gamma/\sim$,
where the equivalence relation $\sim$ is defined by the gauge transformation
(20).

##  III Boundary data

Let us focus on specifying the boundary data for further calculations.

It is convenient to choose $\partial B$ such that it takes the same
$3$-geometry as the hypersurface surface $t=ar$, with $|a|<1$, in the
Minkowski spacetime111Indeed, doing this will introduce a future convenience
in constructing a flat $4$-D geometry whose boundary shares the same intrinsic
geometry as $\partial B$. Using this flat $4$-D geometry, we will establish
real critical points, serving as the initial data for the exploration of
complex critical points (see Sec. IV). . With $\tilde{v}:=t+r$, the Minkowski
metric becomes

$\differential s^{2}_{\rm
M}=-\differential\tilde{v}^{2}+2\differential\tilde{v}\differential
r+r^{2}\differential\Omega^{2}.$ (21)

On the surface $t=ar$, the reduced 3-metric is

$\differential s_{\rm M}^{2}=(1-a^{2})\differential
r^{2}+r^{2}\differential\Omega^{2}.$ (22)

Consider the Eddtington-Finkelsein coordinate where the BH metric reads

$\differential s_{\rm BH}^{2}=-f(r)\differential v^{2}+2\differential
v\differential r+r^{2}\differential\Omega^{2}.$ (23)

Assume that the surface in the BH spacetime taking the same 3-metric as (22)
is given by

$v+F(r)=0.$ (24)

With a straightforward calculation, we get

$\frac{\differential F}{\differential
r}=-\frac{1\pm\sqrt{1-(1-a^{2})f(r)}}{f(r)},$ (25)

where $1-(1-a^{2})f(r)>0$ for all $r\geq r_{b}$ and $a^{2}<1$. Here, $r_{b}$
as the radius of the dust ball at the moment of bounce will be chosen as the
radius of $S_{1}$.

If the positive sign is selected, the right hand side of (25) is singular at
$r=r_{\pm}$ where $f(r_{\pm})=0$. Thus, we should choose

  * (1)

$F_{1}(r)$ for $\partial_{-}B_{1}$ so that

$\frac{\differential F_{1}}{\differential
r}=-\frac{1-\sqrt{1-(1-a_{1}^{2})f(r)}}{f(r)},$ (26)

with some $a_{1}$; and

  * (2)

$F_{2}(r)$ for $\partial_{-}B_{2}$ with

$\frac{\differential F_{2}}{\differential
r}=-\frac{1-\sqrt{1-(1-a_{2}^{2})f(r)}}{f(r)},$ (27)

for some $a_{2}$.

With the past boundaries $\partial_{-}B_{1}$ and $\partial_{-}B_{2}$, we get
the corresponding future boundaries $\partial_{+}B_{1}$ and
$\partial_{+}B_{2}$ by the time reversal symmetry.

It is easy to check the outwards normal $\mathfrak{n}_{i}$ of
$\partial_{-}B_{i}$ for $i=1,2$ is

$\displaystyle\mathfrak{n}_{i}=$
$\displaystyle-\frac{1-\sqrt{1-(1-a_{i}^{2})f(r)}}{f(r)}\,\partial_{v}+$ (28)
$\displaystyle\sqrt{1-(1-a_{i}^{2})f(r)}\,\partial_{r}.$

Their lengths for $\mathfrak{n}_{i}$ is

$\|\mathfrak{n}_{i}\|^{2}=-(1-a_{i}^{2})<0.$ (29)

Thus, the boundary $\partial_{-}B$ is spacelike.

Consider the intersection of $\partial_{-}B_{1}$ and $\partial_{+}B_{1}$. With
their dihedral angle denoted by $\theta_{1}$, we have

$\cosh(\theta_{1})=-1+\frac{2}{(1-a_{1}^{2})f(r_{1})}.$ (30)

Similarly, we get

$\cosh(\theta_{2})=-1+\frac{2}{(1-a_{2}^{2})f(r_{3})},$ (31)

where $\theta_{2}$ is the dihedral angle between $\partial_{-}B_{2}$ and
$\partial_{+}B_{2}$. Eqs. (30) and (31) establish a connection between the
parameters $a_{i}$ and geometrical properties of the region $B$. In our
practical numerical computation, we will select $\theta_{1}$ and $\theta_{2}$
as the parameters to set, rather than $a_{1}$ and $a_{2}$.

###  III.1 Values of $p_{\ell}$ and $\xi_{\ell n}$

In (10), the quantities $p_{\ell}$ and $\xi_{\ell n}$ are related to the
expectation value of $\hat{p}^{k}_{s,\ell}$, which denotes the area vector of
the surface for the flux. In our work, this surface corresponds to a face
within the triangulation of $\partial B$, which serves as a discrete
approximation of the geometry of $\partial B$. Therefore, determining
$p_{\ell}$ and $\xi_{\ell n}$ involves computing the area vectors of surfaces
within the discrete geometry, constructed as outlined below.

Our choice of the boundary $\partial B$ ensures that $\partial_{-}B_{i}$ for
$i=1,2$ have the same intrinsic geometries as that of the hypersurfaces
$t=a_{i}r$ in the Minkowski spacetime. This fact implies the existence of an
embedding $\mathfrak{i}$ of $\partial_{-}B$ into the Minkowski spacetime, such
that the metric of $\mathfrak{i}(\partial_{-}B)$ reduced from the Minkowski
metric is the same as the metric of $\partial_{-}B$ reduced from the BH
metric. On the other hand, the triangulation for $B_{-}$ (refer to Sec. II.2)
will reduce a triangulation for the boundary $\partial_{-}B$, and the boundary
triangulation is completely determined by points $\mathfrak{p}_{a}$,
$\mathfrak{p}_{a^{\prime}}$ and $\mathfrak{p}_{a^{\prime\prime}}$ for
$a=1,2,3,4$ and their connectivity structure. We thus apply the embedding
$\mathfrak{i}$ to these points, mapping them into the Minkowski spacetime.
Then, we connect these points $\mathfrak{i}(\mathfrak{p}_{a})$,
$\mathfrak{i}(\mathfrak{p}_{a^{\prime}})$ and
$\mathfrak{i}(\mathfrak{p}_{a^{\prime\prime}})$ in the Minkowski spacetime
with straight lines, to replicate the connecting structure present in the
triangulation of $\partial_{-}B$. This procedure leads to a discrete geometry
approximating the intrinsic geometry of $\mathfrak{i}(\partial_{-}B)$ and,
consequently, $\partial_{-}B$. Using the discrete geometry, we can calculate
$p_{\ell}$ and $\xi_{\ell f}$ by

$p_{\ell}\xi_{\ell n}^{\dagger}\vec{\sigma}\xi_{\ell
n}=t\frac{\vec{A}_{\ell}}{\beta\kappa\hbar}$ (32)

where $\vec{A}_{\ell}$ is the area vector of the triangle dual to the link
$\ell$. Since $p_{\ell}$ and $\xi_{\ell n}$ for all links $\ell$ connecting to
each node $n$ are determined by a tetrahedron, the closure condition

$\sum_{\ell\text{ at }n}\epsilon_{\ell n}p_{\ell}\xi_{\ell
n}^{\dagger}\vec{\sigma}\xi_{\ell n}=0$ (33)

must be satisfied, where the sign factors $\epsilon_{\ell n}$ depends on the
direction chosen for $\vec{A}_{\ell}$.

Let $\mathfrak{r}:\partial_{-}B\to\partial_{+}B$ denote the time reversal
symmetry. Owing to this symmetry, $\partial_{+}B$ exhibits the same intrinsic
geometry as that of $\partial_{-}B$. Thus, if $\bar{e}_{a}^{i}$ is the cotriad
field in $\partial_{-}B$, the cotriad field in $\partial_{+}B$ could be
$\mathfrak{r}^{*}(\bar{e}_{a}^{i})$, i.e., the push-forward of
$\bar{e}_{a}^{i}$. Here the push-forward of the forms $\bar{e}_{a}^{i}$ is
regarded as their pull-back by $\mathfrak{r}^{-1}$. Despite the time reversal
symmetry, the 3-orientation of $\partial_{+}B$, induced by the outward normal
field, is opposite to $\mathfrak{r}^{*}(o_{-})$, namely, the push forward of
the reduced orientation $o_{-}$ on $\partial_{-}B$. Furthermore, as discussed
in Appendix C, in the spin foam model, the densitized triad should be related
to the triad by $E^{a}_{i}=ee^{a}_{i}$ with $e\equiv\det(e_{b}^{j})$. It is
dual to a 2-form $E^{a}_{i}\varepsilon_{abc}$, whose surface integration,
defining the flux, is 3-orientation dependent (refer to (142)). Consequently,
the flux at each node $\mathfrak{r}(\ell)\mathfrak{r}(n)\subset\partial_{+}B$
becomes the opposite of the flux at the corresponding $(\ell,n)$. In other
words,

$\langle g_{\ell}|\hat{p}_{n,\ell}^{k}|g_{\ell}\rangle=-\langle
g_{\mathfrak{r}(\ell)}|\hat{p}_{\mathfrak{r}(n),\mathfrak{r}(\ell)}^{k}|g_{\mathfrak{r}(\ell)}\rangle$
(34)

with $\hat{p}_{n,\ell}^{k}=\hat{p}_{s,\ell}^{k}$ for $n=s_{\ell}$ and
$\hat{p}_{n,\ell}^{k}=\hat{p}_{t,\ell}^{k}$ for $n=t_{\ell}$. Combining (34)
with (10), we get

$|p_{\mathfrak{r}(\ell)}|\xi_{\mathfrak{r}(\ell)\mathfrak{r}(n)}^{\dagger}\vec{\sigma}\xi_{\mathfrak{r}(\ell)\mathfrak{r}(n)}=-|p_{\ell}|\xi_{\ell
n}^{\dagger}\vec{\sigma}\xi_{\ell n},$ (35)

where $p_{\ell}$ and $\xi_{\ell n}$ are the data associated with
$\ell,n\subset\partial_{-}B$, and $p_{\mathfrak{r}(\ell)}$ and
$\xi_{\mathfrak{r}(\ell)\mathfrak{r}(n)}$ are the corresponding data for
$\mathfrak{r}(\ell)$ and $\mathfrak{r}(n)$. The equation (35) implies

$|p_{\mathfrak{r}(\ell)}|=|p_{\ell}|,\quad\xi_{\mathfrak{r}(\ell)\mathfrak{r}(n)}=e^{ix}J\xi_{\ell
n},$ (36)

where $x$ is an arbitrary phase, and $J\xi$ for
$\xi=(\xi_{1},\xi_{2})^{T}\in\mathbb{C}^{2}$ is defined by

$J\xi=(-\overline{\xi_{2}},\overline{\xi_{1}})^{T}.$ (37)

Given that the phase $x$ in (36) does not affect the outcomes in our paper, we
can set it to zero, namely, we have

$|p_{\mathfrak{r}(\ell)}|=|p_{\ell}|,\quad\xi_{\mathfrak{r}(\ell)\mathfrak{r}(n)}=J\xi_{\ell
n}.$ (38)

In our computation, we choose the surfaces
$\mathfrak{i}(\partial_{-}B)=\mathfrak{i}(\partial_{-}B_{1})\cup\mathfrak{i}(\partial_{-}B_{2})$
as

$\displaystyle\mathfrak{i}(\partial_{-}B_{1}):\ g_{1}(t,r)$
$\displaystyle=t-a_{1}(r-r_{1})=0,$
$\displaystyle\mathfrak{i}(\partial_{-}B_{2}):\ g_{2}(t,r)$
$\displaystyle=t-a_{2}(r-r_{2})-a_{1}(r_{2}-r_{1})=0.$

It is reminded that $\mathfrak{i}$, as defined in the paragraph before (32),
denotes the embedding of $\partial_{-}B$ into the Minkowski spacetime. Then,
we have

$\displaystyle\mathfrak{i}(S_{1})$
$\displaystyle=\\{(t_{1},r_{1},\theta,\phi)\in\Sigma_{1},\theta\in(0,\pi),\phi(0,2\pi)\\}$
(39) $\displaystyle\mathfrak{i}(S_{2})$
$\displaystyle=\\{(t_{2},r_{2},\theta,\phi)\in\Sigma_{1},\theta\in(0,\pi),\phi(0,2\pi)\\}$
$\displaystyle\mathfrak{i}(S_{3})$
$\displaystyle=\\{(t_{3},r_{3},\theta,\phi)\in\Sigma_{2},\theta\in(0,\pi),\phi(0,2\pi)\\}$

with

$\displaystyle t_{1}$ $\displaystyle=0,$ (40) $\displaystyle t_{2}$
$\displaystyle=a_{1}(r_{2}-r_{1}),$ $\displaystyle t_{3}$
$\displaystyle=a_{2}(r_{3}-r_{2})+a_{1}(r_{2}-r_{1}).$

Introducing

$\displaystyle\vec{R}_{1}=$
$\displaystyle\left(\sqrt{\frac{8}{9}},0,-\frac{1}{3}\right),\vec{R}_{2}=\left(-\sqrt{\frac{2}{9}},\sqrt{\frac{2}{3}},-\frac{1}{3}\right),$
(41) $\displaystyle\vec{R}_{3}=$
$\displaystyle\left(-\sqrt{\frac{2}{9}},-\sqrt{\frac{2}{3}},-\frac{1}{3}\right),\vec{R}_{4}=\left(0,0,1\right),$

we set the coordinates of $\mathfrak{i}(\mathfrak{p}_{a})$,
$\mathfrak{i}(\mathfrak{p}_{a^{\prime}})$ and
$\mathfrak{i}(\mathfrak{p}_{a^{\prime\prime}})$ in Minkowski spacetime to be

$\displaystyle\mathfrak{i}(\mathfrak{p}_{a})$
$\displaystyle=(t_{1},r_{1}\vec{R}_{a}),\
\mathfrak{i}(\mathfrak{p}_{a^{\prime}})=(t_{1},-r_{2}\vec{R}_{a}),$ (42)
$\displaystyle\mathfrak{i}(\mathfrak{p}_{a^{\prime\prime}})$
$\displaystyle=(t_{1},-r_{3}\vec{R}_{a}),\ \forall a=1,2,3,4.$

###  III.2 Relate $\phi_{\ell}$ to extrinsic curvature

According to (7), we need to calculate the holonomy within the boundary
geometry to get $\phi_{\ell}$. Appendix C provides a detailed computation,
with key aspects worth noting. At first, as shown in Han and Zhang (2013), the
value of the spin foam action at the real critical points differs from the
Regge action by a factor $\mathrm{sgn}(V_{4}(v))$. This extra factor suggests
that the volume form in classical action corresponding to the spin foam model
is $\det(\mathfrak{e})\differential^{4}x$ where $\mathfrak{e}$ denotes the
coframe field. Starting from the Einstein-Hilbert action with this volume
form, we establish the corresponding Ashtekar variables and finally determine
the holonomy. Additionally, the link $\ell$ in the holonomy is chosen to be a
geodesic for simplicity. Finally, we assume that the extrinsic curvature
$K_{ab}$ along the link $\ell$ takes the form

$K^{ab}(t)=\frac{\xi(t)}{\|\dot{\ell}_{t}\|^{2}}\dot{\ell}^{a}_{t}\dot{\ell}^{b}_{t}$
(43)

with $\xi$ being a field on $\ell$. This assumption makes the computation
tractable. In LQG, quantities are often smeared by integration. Therefore, the
simplified expression (43) can be interpreted as an approximation of the true
behavior of the extrinsic curvature, which reflects the integrated properties
of the actual extrinsic curvature over the link. In the context described
above, we obtain the expression for $\phi_{\ell}$ as follows:

$\phi_{\ell}=\theta_{\ell}-\beta\mathrm{sgn}(p_{\ell})\mathrm{sgn}(e)\Xi_{\ell},$
(44)

where $\theta_{\ell}$ is the twist angle, $\Xi_{\ell}$ is related to the
extrinsic curvature and $\mathrm{sgn}(e)\equiv\mathrm{sgn}(\det(e_{a}^{i}))$
is computed with any right-hand coordinate system.

The twist angle $\theta_{\ell}$ will be fixed by choosing some certain frame
for the two tetrahedra sharing the triangle dual to $\ell$. The way to fix
this angle will be introduced in more details during the discussion of our
numerical computation. For now, we will focus only on $\Xi_{\ell}$ which, due
to the assumption (43), reads

$\Xi_{\ell}=\int_{\ell}K_{ab}(s)\dot{\ell}^{a}_{s}\dot{\ell}^{b}_{s}\differential
s$ (45)

where $s$ is the parameter of length of $\ell$, $\dot{\ell}_{s}^{a}$ is the
tangent vector $\ell$ with respect to the parameter $s$.

Turning our attention to the continuous geometry on $\partial_{-}B$, a naive
approach to relate $\Xi_{\ell}$ to this geometry involves a direct application
of (182), with $K_{ab}$ being the extrinsic curvature of the continuous
geometry. However, this approach cannot fully encode the information of
$K_{ab}$ within $\Xi_{\ell}$. Indeed, in the continuous geometry, the field
$K_{ab}$ undergoes variations at distinct points, in contrast to the Regge
geometry where the extrinsic curvature behaves as a distribution with a
vanishing value everywhere except on the interfacing triangle $t_{\ell}$,
where it remains constant. Thus, integrating $K_{ab}$ along a single edge, as
in (182), is limited to capturing the information of $K_{ab}$ solely on that
specific edge. Recognizing the limitations of this naive approach, we will
abandon it and adopt a more refined strategy. The new approach takes advantage
of the triangulation for $\partial_{-}B$. Consider two tetrahedra
$T_{1}=\mathfrak{q}_{1}\mathfrak{q}_{2}\mathfrak{q}_{3}\mathfrak{q}_{4}$ and
$T_{2}=\mathfrak{q}_{2}\mathfrak{q}_{3}\mathfrak{q}_{4}\mathfrak{p}_{5}$,
sharing the interfacing triangle
$t_{\ell}=\mathfrak{q}_{2}\mathfrak{q}_{3}\mathfrak{q}_{4}$. Here $q_{i}$ for
$i=1,2,3,4,5$ represents a point in the triangulation of $\partial_{-}B$.
Introducing the “centroids” $\mathfrak{c}_{1}\in T_{1}$ and
$\mathfrak{c}_{2}\in T_{2}$, we define a convex region
$\mathfrak{c}_{1}\mathfrak{q}_{2}\mathfrak{q}_{3}\mathfrak{q}_{4}\mathfrak{c}_{2}$.
This region is foliated into a family of two surfaces $\tau\mapsto
A(\tau)\subset\mathbb{R}^{2}$ with $A(0)=t_{\ell}$. The average extrinsic
curvature $\overline{K}(\tau)$ is defined on each surface as

$\overline{K}(\tau):=\|A_{\tau}\|^{-1}\int_{A_{\tau}}(K_{ab}n^{a}n^{b})n^{c}\varepsilon_{cde}\differential
x^{d}\differential x^{e},$ (46)

where $\|A_{\tau}\|$ is the area of $A_{\tau}$ and $n^{a}$ is the unit normal
of $A_{\tau}$. $\overline{K}(\tau)$ in (46) captures the information
$K_{ab}n^{a}n^{b}$ on the surface $A_{\tau}$. Choosing a link $\ell$
connecting $\mathfrak{q}_{1}$ and $\mathfrak{q}_{2}$ and transversely
intersecting each $A_{\tau}$, we relate $\Xi_{\ell}$ to the extrinsic
curvature by

$\Xi_{\ell}:=\int_{\ell}\overline{K}(\tau(s))\differential s.$ (47)

This equation _defines_ a connection between the extrinsic curvature and the
data $\Xi_{\ell}$ inherent in the coherent state. this definition reverts to
(45) in the specific case of Regge geometry (refer to the discussion in
Appendix C.4). Importantly, as (47) defines $\Xi_{\ell}$ through the average
extrinsic curvature, it effectively encapsulates information about the
extrinsic curvature across the entire region
$\mathfrak{c}_{1}\mathfrak{q}_{2}\mathfrak{q}_{3}\mathfrak{q}_{4}\mathfrak{q}_{2}$.

Our actual computation adopts the coordinate give by
$(x,y,z)=(r\sin\theta\cos\phi,r\sin\theta\sin\phi,r\cos\theta)$ for
simplicity. Within this coordinate, we choose
$\mathfrak{c}_{1}:=\frac{1}{4}(\mathfrak{q}_{1}+\mathfrak{q}_{2}+\mathfrak{q}_{3}+\mathfrak{q}_{4})$
and
$\mathfrak{c}_{2}:=\frac{1}{4}(\mathfrak{q}_{2}+\mathfrak{q}_{3}+\mathfrak{q}_{4}+\mathfrak{q}_{5})$.
The region
$\mathfrak{c}_{1}\mathfrak{q}_{2}\mathfrak{q}_{3}\mathfrak{q}_{4}\mathfrak{c}_{2}$
is defined as the polyhedron bounded by plane faces. These plane faces are
expressed as surfaces in the coordinate system $(x,y,z)$ through linear
equations. In the foliation of
$\mathfrak{c}_{1}\mathfrak{q}_{2}\mathfrak{q}_{3}\mathfrak{q}_{4}\mathfrak{c}c_{2}$,
the slices $A_{\tau}$ are chosen as plane surfaces parallel to $t_{\ell}$. The
link $\ell$, connecting $\mathfrak{c}_{1}$ and $\mathfrak{c}_{2}$, is selected
as the straight line $t\mapsto(1-t)\mathfrak{c}_{1}+t\mathfrak{c}_{2}$ in this
coordinate system.

Due to the time reversal symmetry
$\mathfrak{r}:\partial_{-}B\to\partial_{+}B$, the extrinsic curvature
$K_{ab}^{+}$ in $\partial_{+}B$ is just $\mathfrak{r}^{*}K_{ab}^{-}$, i.e.,
the push forward of the extrinsic curvature $K_{ab}^{-}$ in $\partial_{-}B$.
Therefore, for a given $\ell\subset\partial_{-}B$ carrying $\Xi_{\ell}$, the
value of $\Xi_{\mathfrak{r}(\ell)}$ for
$\mathfrak{r}(\ell)\subset\partial_{+}B$ is the same as $\Xi_{\ell}$, i.e.,

$\Xi_{\mathfrak{r}(\ell)}=\Xi_{\ell}.$ (48)

###  III.3 Superpose the coherent states

Let us firstly summarize the results given in Secs. III.1 and III.2. Given a
link $\ell$ in $\partial_{-}B$, we can calculate $p_{\ell}$, $\xi_{\ell
s_{\ell}}$, $\xi_{\ell t_{\ell}}$, $\theta_{\ell}$ and $\Xi_{\ell}$, to get
the coherent state label $\bar{g}_{\ell}$ as

$\bar{g}_{\ell}=n(\xi_{\ell
s_{\ell}})e^{(-ip_{\ell}+\theta_{\ell}-\beta\mathrm{sgn}(\bar{e})\Xi_{\ell})\tau_{3}}n(\xi_{\ell
t_{\ell}})^{-1}$ (49)

where we choose $p_{\ell}>0$ without loss of generality and $\bar{e}_{a}^{i}$
denotes the triad field in $\partial_{-}B$. Then, by choosing

$p_{\mathfrak{r}(\ell)}=-p_{\ell}$

that satisfies (38), we get the coherent state label
$\tilde{g}_{\mathfrak{r}(\ell)}$

$\tilde{g}_{\mathfrak{r}(\ell)}=n(J\xi_{\ell
s_{\ell}})e^{(-ip_{\ell}-\theta_{\ell}+\beta\mathrm{sgn}(\bar{e})\Xi_{\ell})\tau_{3}}n(J\xi_{\ell
t_{\ell}})^{-1}$ (50)

where we combine (38), (44), (48) and (12), and consider the discussion in the
paragraph before (34) showing that
$\mathrm{sgn}(\mathfrak{r}^{*}\bar{e})=-\mathrm{sgn}(\bar{e})$.

To maintain a consistent global orientation for the two complexes dual to the
triangulation of $B$, the boundary face in $\partial_{+}B$ corresponding to a
face $f\subset\partial_{-}B$ is represented by $\mathfrak{r}(f)^{-1}$. In
simpler terms, this means that the corresponding face includes the same points
as $\mathfrak{r}(f)$ but with the opposite orientation. As a result, for a
link $\ell\subset\partial_{-}B$, its counterpart in $\partial_{+}B$ is given
by $\mathfrak{r}(\ell)^{-1}$. By (50), the coherent state associated with
$\mathfrak{r}(\ell)^{-1}$ is labelled by

$\displaystyle\tilde{g}_{\mathfrak{r}(\ell)^{-1}}=\tilde{g}_{\mathfrak{r}(\ell)}^{-1}$
(51) $\displaystyle=$ $\displaystyle n(\xi_{\ell
t_{\ell}})e^{(-ip_{\ell}-\theta_{\ell}+\beta\mathrm{sgn}(\bar{e})\Xi_{\ell})\tau_{3}}n(\xi_{\ell
s_{\ell}})^{-1}.$

As discussed in (14), we are considering superposition of various coherent
states resembling the same geometry. Following the convention introduced in
(14), the values for $g_{\ell}^{(++)}$ are obtained from (49) and (51). Now,
let us shift our attention to $g_{\ell}^{(s_{1}s_{2})}$ for other signs of
$s_{1}$ and $s_{2}$. Consider the transformation
$\bar{e}_{a}^{i}\mapsto-\bar{e}_{a}^{i}$ in $\partial_{-}B$. This
transformation does not alter the densitized triad field $E^{a}_{i}$ in
$\partial_{-}B$, and consequently, it does not change $\xi_{\ell n}$ and
$p_{\ell}$. Furthermore, the spin connection $\Gamma_{a}^{i}$ and the
extrinsic curvature $K_{ab}$ remain unaffected. As a result, $\theta_{\ell}$
and $\Xi_{\ell}$ remains constant under this transformation. However, due to
the relationship
$\mathrm{sgn}(\det(\bar{e}_{a}^{i}))=-\mathrm{sgn}(\det(-\bar{e}^{i}_{a}))$,
$\phi_{\ell}$ will be transformed to

$\phi_{\ell}\mapsto\phi_{\ell}^{(-)}=\theta_{\ell}+\beta\mathrm{sgn}(\bar{e})\Xi_{\ell}.$
(52)

This leads to the following transformation for the coherent state label
$\bar{g}_{\ell}$,

$\bar{g}_{\ell}\mapsto\bar{g}_{\ell}^{(-)}:=n(\xi_{\ell
s_{\ell}})e^{(-ip_{\ell}+\theta_{\ell}+\beta\mathrm{sgn}(\bar{e})\Xi_{\ell})\tau_{3}}n(\xi_{\ell
t_{\ell}})^{-1}.$ (53)

Similarly, under the transformation
$\mathfrak{r}^{*}(\bar{e}_{a}^{i})\mapsto-\mathfrak{r}^{*}(\bar{e}_{a}^{i})$,
$\tilde{g}_{\mathfrak{r}(\ell)^{-1}}$ for all $\mathfrak{r}(\ell)$ transformed
as

$\displaystyle\tilde{g}_{\mathfrak{r}(\ell)^{-1}}\mapsto\tilde{g}_{\mathfrak{r}(\ell)^{-1}}^{(-)}:=n(\xi_{\ell
t_{\ell}})e^{(-ip_{\ell}-\theta_{\ell}-\beta\mathrm{sgn}(\bar{e})\Xi_{\ell})\tau_{3}}n(\xi_{\ell
s_{\ell}})^{-1}.$ (54)

As a summary, the list of coherent state labels
$\\{g_{\ell}^{(\tilde{s}\bar{s})}\\}$ for links lying entirely in either
$\partial_{+}B$ or $\partial_{-}B$ is given by

$\\{\tilde{g}_{\mathfrak{r}(\ell)^{-1}}^{(\tilde{s})}\\}_{\mathfrak{r}(\ell)\subset\partial_{+}B}\cup\\{\bar{g}_{\ell}^{(\bar{s})}\\}_{\ell\subset\partial_{-}B},$
(55)

where we defined
$\tilde{g}_{\mathfrak{r}(\ell)^{-1}}^{(+)}=\tilde{g}_{\mathfrak{r}(\ell)^{-1}}$
and $\bar{g}_{\ell}^{(+)}=\bar{g}_{\ell}$.

###  III.4 Boundary states associated with links at corners

It is crucial to highlight that the links passing through the spheres $S_{1}$
and $S_{3}$ are excluded from (55). Now, let us discuss the coherent state
labels for these links individually. To maintain the time reversal symmetry,
we will designate such a link in a way that it can be divided into two parts,
$\ell=\bar{\ell}\cup\tilde{\ell}$, where
$\tilde{\ell}=\mathfrak{r}(\bar{\ell})^{-1}$. In considering the values of
$g^{(\tilde{s}\bar{s})}_{\ell}$ for different signs of $\tilde{s},\bar{s}$, we
correspond to intrinsic geometries given by
$(\tilde{s}\tilde{e}_{a}^{i},\bar{s}\bar{e}_{a}^{i})=(\tilde{s}\mathfrak{r}^{*}(\bar{e}_{a}^{i}),\bar{s}\bar{e}_{a}^{i})$
(refer to (14) for the convention). For each of the segments $\bar{\ell}$ and
$\tilde{\ell}$, there exist associated coherent states
$\psi^{t}_{\tilde{g}_{\tilde{\ell}}^{(\tilde{s})}}$ and
$\psi^{t}_{\bar{g}_{\bar{\ell}}^{(\bar{s})}}$. Since $\bar{\ell}$ and
$\tilde{\ell}$ entirely lie in $\partial_{-}B$ and $\partial_{+}B$
respectively, the coherent state labels
$\tilde{g}_{\tilde{\ell}}^{(\tilde{s})}$ and
$\bar{g}_{\bar{\ell}}^{(\bar{s})}$ can be obtained applying the algorithm
introduced in Secs. III.1, III.2, and III.3. The results of the labbels can be
parametrized as

$\displaystyle\tilde{g}_{\tilde{\ell}}^{(\tilde{s})}=$ $\displaystyle
n(\xi_{\tilde{\ell}s_{\tilde{\ell}}})e^{(-ip_{\tilde{\ell}}+\phi_{\tilde{\ell}}^{\tilde{s}})}n(\xi_{\tilde{\ell}t_{\tilde{\ell}}}),$
(56) $\displaystyle\bar{g}_{\bar{\ell}}^{(\bar{s})}=$ $\displaystyle
n(\xi_{\bar{\ell}s_{\bar{\ell}}})e^{(-ip_{\bar{\ell}}+\phi_{\bar{\ell}}^{\bar{s}})}n(\xi_{\bar{\ell}t_{\bar{\ell}}}),$

where it is considered that the signs $\tilde{s}$ and $\bar{s}$ only affect
the value of the angle $\phi$, as indicated by (49), (50), (53), and (54). Due
to he time reversal symmetry, we have

$\displaystyle p_{\tilde{\ell}}$ $\displaystyle=p_{\bar{\ell}},\
\xi_{\bar{\ell}s_{\bar{\ell}}}=\xi_{\tilde{\ell}t_{\tilde{\ell}}},\
\xi_{\bar{\ell}t_{\bar{\ell}}}=\xi_{\tilde{\ell}s_{\tilde{\ell}}}.$ (57)

Thanks to the coherent states associated to the segments, it is natural to
consider the boundary state associated with $\ell$ as the combination of these
individual coherent states. Specifically, the boundary state associated with
$\ell$ is given by

$\displaystyle\int_{\mathrm{SU(2)}}\differential\mu_{H}(h)\psi^{t}_{\bar{g}_{\bar{\ell}}^{(\bar{s})}}(h_{\bar{\ell}}h^{-1})\psi^{t}_{\tilde{g}_{\tilde{\ell}}^{(\tilde{s})}}(hh_{\tilde{\ell}})=\psi^{2t}_{\bar{g}_{\bar{\ell}}^{(\bar{s})}\tilde{g}_{\tilde{\ell}}^{(\tilde{s})}}(h_{\ell}).$
(58)

where $h_{\ell}=h_{\bar{\ell}}h_{\tilde{\ell}}$ is the holonomy for the entire
link. Noticing the change of parameter $t\mapsto 2t$, we find that the state
(58) is approximated by $|g_{\ell}^{(\tilde{s}\bar{s})}\rangle$ (refer to (8))
with

$\displaystyle g_{\ell}^{(\tilde{s}\bar{s})}=$ $\displaystyle n(\xi_{\ell
s_{\ell}})e^{(-ip_{\ell}+\phi_{\ell}^{(\tilde{s}\bar{s})})}n(\xi_{\ell
t_{\ell}})^{-1},$ (59) $\displaystyle\xi_{\ell s_{\ell}}=$
$\displaystyle\xi_{\bar{\ell}s_{\bar{\ell}}}=\xi_{\ell t_{\ell}},$
$\displaystyle p_{\ell}=$
$\displaystyle\frac{1}{2}(p_{\tilde{\ell}}+p_{\bar{\ell}})=p_{\tilde{\ell}},$
$\displaystyle\phi_{\ell}^{(\tilde{s}\bar{s})}=$
$\displaystyle\phi_{\bar{\ell}}^{(\bar{s})}+\phi_{\tilde{\ell}}^{(\tilde{s})}.$

Indeed, rather than the state (58), we prefer to take the boundary state
associated with $\ell$ as $\psi^{t}_{g_{\ell}^{(\tilde{s}\bar{s})}}$. Even
though both of the two state can be approximated by
$|g_{\ell}^{(\tilde{s}\bar{s})}\rangle$, the choice of
$\psi^{t}_{g_{\ell}^{(\tilde{s}\bar{s})}}$ aligns with the states chosen for
other edges lying entirely in $\partial_{\pm}B$. Finally, by definition of
$\phi_{\bar{\ell}}^{(\bar{s})}$ and $\phi_{\tilde{\ell}}^{(\tilde{s})}$, we
get

$\displaystyle\phi_{\ell}^{(\tilde{s}\bar{s})}=\left\\{\begin{array}[]{cl}0,&\tilde{s}=\bar{s},\\\
2\tilde{s}\beta\mathrm{sgn}(\bar{e})\Xi_{\bar{\ell}},&\bar{s}=-\tilde{s}.\end{array}\right.$
(60)

##  IV Computing the critical point

In the previous sections, we define the spin foam amplitude and construct the
boundary data to study the dynamics of the $B$ region. In this section, we
will explain how to apply the method of _complex critical point_ to
numerically compute the spin amplitude Han  _et al._ (2020, 2022, 2023b,
2024). This method involves obtaining the complex critical point for each set
$\\{g_{\ell}^{(\tilde{s}\bar{s})}\\}_{\ell\subset\partial B}$ with
$\tilde{s},\,\bar{s}=\pm 1$. Since the procedure for each
$\\{g_{\ell}^{(\tilde{s}\bar{s})}\\}_{\ell}$ is the same, we will use
$g_{\ell}$ to represent either of $g_{\ell}^{(\tilde{s}\bar{s})}$ in the
following.

In the computation of the spin foam action (19), we choose the parameter $t$
as

$t=\frac{\beta\kappa\hbar}{\Delta}.$ (61)

where $\Delta$ is a typical area. Then, the effective action for a boundary
face can be expressed in terms of $\|\vec{A}\|$ as

$\displaystyle S_{f}=$
$\displaystyle-\frac{\beta\kappa\hbar}{2\Delta}\left(j_{f}-\frac{\|\vec{A}_{\ell_{f}}\|}{\beta\kappa\hbar}\right)^{2}$
(62)
$\displaystyle+i\,\mathrm{sgn}(p_{\ell_{f}})j_{f}\phi_{\ell_{f}}+j_{f}\mathcal{S}_{o}^{(f)}$

where (32) is applied, and all of the remaining terms are absorbed into
$\mathcal{S}_{o}^{(f)}$ for convenience. Replacing $j_{f}\mapsto t^{-1}j_{f}$,
we get

$\displaystyle S_{f}=$
$\displaystyle\frac{1}{t}\Bigg{[}-\frac{1}{2}\left(j_{f}-\frac{\|\vec{A}_{\ell_{f}}\|}{\Delta}\right)^{2}$
(63)
$\displaystyle+i\,\mathrm{sgn}(p_{\ell_{f}})j_{f}\phi_{\ell_{f}}+j_{f}\mathcal{S}_{o}^{(f)}\Bigg{]}.$

This allows the action for a boundary face to take the form $t^{-1}S_{f}$. It
is noteworthy that the substitution $j_{f}\to t^{-1}j_{f}$ can similarly
modify the action for an internal face, rendering it in the form
$t^{-1}S_{f}$. Consequently, the amplitude $A$ becomes

$A=\sum_{j}\int\differential g\differential
ze^{\frac{1}{t}S[j,g,z]}\cong\int\differential j\differential g\differential
ze^{\frac{1}{t}S[j,g,z]},$ (64)

where the Poisson resummation formula is employed to approximate the summation
by an integral (see Han  _et al._ (2022) for more details). Thank to (64), we
are allowed to discuss the asymptotic behavior of $A$ as $t^{-1}\to\infty$.
This involves the complex critical point method.

###  IV.1 Complex critical point method

As shown in (64), we are concerned with the integral of the form

$A_{t}=\int_{\mathcal{C}}e^{\frac{1}{t}S(x)}\differential\mu(x)$

over some $N-$dimensional configuration space $\mathcal{C}$ with some
parameter $t$. The behavior of $A_{t}$ as $t\to 0$ is guided by the critical
points of $S$, where a critical point is a point $x_{o}\in\mathcal{C}$
satisfying

$\partial_{x}S(x_{o})=0.$

It is worth noting that, in general, the equations $\partial_{x}S=0$ may lack
solutions when $x$ is restricted to real values. In such cases, the search
extends to complex critical points. A complex critical point lies in the
complexification $\mathcal{C}^{\mathbb{C}}$ of $\mathcal{C}$ and satisfies
$\partial_{z}S=0$, where $S$ needs to be analytically continued to
$\mathcal{C}^{\mathbb{C}}\ni z$.

In our computation, we will introduce some external parameters
$\chi\in[\chi_{o},\chi_{1}]$ as in Han  _et al._ (2022), to rewrite the
integral as

$A_{t}(\chi)=\int_{\mathcal{C}}e^{\frac{1}{t}S(x,\chi)}\differential\mu(x).$
(65)

Here $S(x,\chi)$ is assumed to be an analytical function with respect to
$\chi$. Additionally, $\chi$ should be selected in a way that $S(x,\chi_{1})$
precisely describes our model of interest. Simultaneously, for
$\chi=\chi_{o}$, $S(x,\chi_{o})$ may represent a distinct model, but it
possesses a real critical point. Let $Z(\chi)$ be the solution to
$\partial_{x}S(x,\chi)=0$ such that $Z(\chi_{o})$ is the real critical point
of $S(x,\chi_{o})$. Here, the real critical point is defined as the solution
to both $\partial_{x}S(x,\chi_{o})=0$ and $\mathrm{Re}(S(x,\chi_{o}))=0$. As
$\chi$ deviates away from $\chi_{o}$, $Z(\chi)$ moves away from $\mathcal{C}$
to $\mathcal{C}^{\mathbb{C}}$ analytically. Consequently,
$Z(\chi_{1})\in\mathcal{C}^{\mathbb{C}}$ will be a complex critical point
contribute the integral in our initial model (see Fig. 7 for the
illustration). With the complex critical point, the small-$t$ asymptotic
expansion for the integral in $A_{t}$ can be established as

$\displaystyle A_{t}=\left(2\pi
t\right)^{\frac{N}{2}}\frac{e^{\frac{1}{t}S(Z(\chi),\chi)}}{\sqrt{\det[-H(Z(\chi),\chi)]}}\left[1+O(t)\right],$
(66)

where the right hand of (66) involves the action $S(Z(\chi),\chi)$ and Hessian
$H(Z(\chi),\chi))$ at the complex critical point. Furthermore, the real part
of $S(Z(\chi),\chi)$ satisfies the condition

$\operatorname{Re}\left(S(Z(\chi),\chi)\right)\leq-C|\operatorname{Im}(Z(\chi))|^{2}.$
(67)

where $C$ is a positive constant Hormander (1983). This condition implies that
in (65), the purely imaginary phase can only occur at the real critical point,
where $\operatorname{Im}(Z)=0$ and $\chi=\chi_{o}$. When $\chi$ deviates from
$\chi_{o}$, causing $\operatorname{Im}(Z)$ finite and
$\operatorname{Re}(\mathcal{S})$ to become negative, the result in (65) is
exponentially suppressed as $t$ approaches $0$.

Figure 7: The real and complex critical points $x_{o}$ and $Z(\chi)$.
$S(z,\chi)$ is analytic extended from the real axis to the complex
neighborhood illustrated by the red disk.

Now let us return to the integral (64) to apply this method. As discussed in
Sec. II.4, the spin foam action is a function on the space
$\Gamma\ni(j_{f},g_{ve},z_{vf})$. However, the integral in (64) is merely
performed over the quotient space $\Gamma/\sim$, where the gauge orbit given
by (20) is divided. This fact implies that when considering the critical
equations,, the derivatives along the gauge orbit should not be taken into
account. Namely, the critical point is a solution to

$\displaystyle\partial_{j_{f_{i}}}S$ $\displaystyle=$ $\displaystyle
0,\quad\partial_{g_{ve}}^{\prime}S=0,\quad\partial_{z_{vf}}^{\prime}S=0$ (68)
$\displaystyle\partial_{j_{f_{b}}}S$ $\displaystyle=$ $\displaystyle 0.$ (69)

where $\partial^{\prime}$ denotes the derivative only along the direction
transversal the gauge orbit.

The first issue for applying the complex critical point method is the choice
of the parameter $\chi\in[\chi_{o},\chi_{1}]$. By definition, $\chi$ is chosen
to connect the concerning model with a model possessing a real critical point.
In the spin foam model, a real critical point $x_{o}$ is constructed from flat
geometry and satisfies $\mathrm{Re}(S)=0$ along with the equations in (68).
However, $x_{o}$ may not fulfill (69), which includes boundary data. Thus, we
need to establish a relationship between the parameter $\chi$ and the boundary
data, so that (69) can also be solved in the model of $S(x,\chi_{o})$.

To maintain the flow of the discussion, the subsequent subsections will
disregard the parameter $\chi$ and focus on obtaining $x_{o}$ from geometry,
as well as searching for complex critical point from $x_{o}$. When the
parameter $\chi$ is considered, the procedure remains similar. We still need
to obtain $x_{o}$ from geometry initially, but search for the critical point
$Z(\chi)$ associated with $S(x,\chi)$ from $x_{0}$ for all
$\chi\in[\chi_{o},\chi_{1}]$222In the practical computation, we search for the
critical point $Z(\chi_{o}+\delta\chi)$ associated with
$S(x,\chi_{o}+\delta\chi)$ with $\delta\chi\ll 1$ from $x_{0}$, then search
for the critical point $Z(\chi_{o}+2\delta\chi)$ associated with
$S(x,\chi_{o}+2\delta\chi)$ from $Z(\chi_{o}+\delta\chi)$, and so forth, until
we obtain the critical point associated with $S(x,\chi_{1})$, i.e., the
solution to both (68) and (69).. The exact definition of $\chi$ will be
provided in Sec. V.2. From now on, $x_{o}$ will be referred to as the real
critical point, while it does not fully satisfy the critical equations (68)
and (69).

###  IV.2 Real critical point

Let us begin by explaining the choice of units in the numerical computation.
in our numerical computations, when referring to an area, it is expressed as
the value taking the typical area $\Delta$ as the unit. Under this unit, the
parameter $t$ is equivalent to the Planck area $\beta\kappa\hbar$.

The absence of (69) enable us to disregard the boundary condition arising from
the coherent state data $p_{\ell}$, $\xi_{\ell}$ and $\phi_{\ell}$.
Nevertheless, in the first step to obtain $x_{o}$ which closely approximates
the final solution to the full critical equations, we choose to in incorporate
the data $p_{\ell}$ and $\xi_{\ell}$.

Given a set of boundary data
$\\{(p_{\ell},\xi_{\ell})\\}_{\forall\ell\in\partial B}$, the typical
procedure to find the critical point involves constructing a flat geometry
that is compatible with the boundary data and extracting the spin foam data
from this geometry. However, since the boundary is closed, a concern arises
regarding whether the geometry satisfying the boundary condition given by
$\\{(p_{\ell},\xi_{\ell})\\}_{\forall\ell\in\partial B}$ is flat. In the case
at hand, the answer is negative. As we will see, we actually can construct a
flat geometry in the region $B_{-}$ that is compatible with the boundary data
on $\partial_{-}B$. However, the extension of this flat geometry to the region
$B_{+}$ does not match the boundary condition on $\partial_{+}B$.

In Sec. III.1, to construct the discrete geometry approximating the intrinsic
geometry of $\partial_{-}B$, we introduced the points
$\mathfrak{i}(\mathfrak{p}_{a})$, $\mathfrak{i}(\mathfrak{p}_{a^{\prime}})$
and $\mathfrak{i}(\mathfrak{p}_{a^{\prime\prime}})$ in the Minkowski spacetime
for all $a=1,2,3,4$. Additionally, in Sec. II, we introduce the triangulation
of the region $B_{-}$ that encapsulates the connecting structure of the points
$\mathfrak{p}_{a}$, $\mathfrak{p}_{a^{\prime}}$ and
$\mathfrak{p}_{a^{\prime\prime}}$. We thus replicate the connections between
the points $\mathfrak{p}_{a}$, $\mathfrak{p}_{a^{\prime}}$ and
$\mathfrak{p}_{a^{\prime\prime}}$, to connect the points
$\mathfrak{i}(\mathfrak{p}_{a})$, $\mathfrak{i}(\mathfrak{p}_{a^{\prime}})$
and $\mathfrak{i}(\mathfrak{p}_{a^{\prime\prime}})$. This process results in a
flat Regge geometry in the region $B_{-}$. Utilizing the geometric
interpretation of the real spin-foam critical points Barrett  _et al._ (2009,
2010); Han and Zhang (2013), we derive a set of data $(j_{f},g_{ve},z_{vf})$
from this geometry for the $B_{-}$ region. According to the definition of the
boundary data $p_{\ell}$ and $\xi_{\ell}$, these data are clearly compatible
with the boundary data on $\partial_{-}B$. Nevertheless, these data are
specifically defined on the 2-complex in $B_{-}$ and solely solve the
equations in (68) associated with $j_{f_{i}}$, $g_{ve}$, and $z_{vf}$ within
$B_{-}$. They even do not solve the equations $\partial_{j_{f_{i}}}S=0$ for
the faces $f_{i}$ across the transition surface $\mathcal{T}$.

Given a vertex $v$, an edge $e$ and a face $f$ contained in $B_{-}$, the time
reversal symmetry $\mathfrak{r}:B_{-}\to B_{+}$ defines the corresponding
objects $\mathfrak{r}(v),\mathfrak{r}(e),\mathfrak{r}(f)$ in $B_{+}$. With the
spin foam data on the 2-complex in $B_{-}$, we assign the spin foam data to
the 2-complex in $B_{+}$ as

$(j_{\mathfrak{r}(f)},g_{\mathfrak{r}(v)\mathfrak{r}(e)},z_{\mathfrak{r}(v)\mathfrak{r}(f)})=\left(j_{f},g_{ve},z_{vf}\right).$
(70)

Then, we get a set of spin foam data for all $v,e,f$ in the entire $B$ region,
referred to as the extended spin foam data. As shown in Appendix D, it is
proved that the extended spin foam data are compatible with the boundary data
on the entire boundary $\partial B$ and solve the equations in (68) for all
relevant variables in the entire $B$ region.

The above discussion outlines the process of obtaining the real critical
point. Nest, let us focus on finding the complex critical point.

###  IV.3 Complex critical point

The absence of (69) indicates that $x_{o}$ does not fully satisfy the critical
equations. In other words, the geometry derived from the real critical point
does not align with the boundary data provided by the coherent state. More
specifically, the resulting geometry exhibits a different dihedral angle
compared to $\Xi_{\ell}$ computed in Sec. III.2. Introducing the dihedral
angles $\Xi_{\ell}$ as an additional boundary condition prevents the
construction of a flat geometry, even locally in $B_{-}$. Therefore, following
the flatness property of spin foam model Donà  _et al._ (2020), we cannot
obtain a real solution to (68) and (69), and have to find complex solutions.

The first step in seeking a complex solution involves parameterizing
$(j_{f},g_{ve},z_{vf})$. The intricacies of this parametrization arise from
that (68) contains only the derivative along the direction transversal to the
gauge orbits. This fact needs us to introduce a gauge fixing in the
parametrization, to ensure that the degrees of freedom along the gauge orbits
are not considered. We denote the real critical point as
$(\mathring{j}_{f},\mathring{g}_{ve},\mathring{z}_{vf})$. Then the following
parametrization of $(j_{f},g_{ve},z_{vf})$, with a gauge fixing condition
taken into account, is applied.

At first, for each vertex $v$, we choose an edge $e_{o}(v)$ and fix
$g_{ve_{o}(v)}$ as

$g_{ve_{o}(v)}=\mathring{g}_{ve_{o}(v)}.$ (71)

This fixes the $\mathrm{SL}(2,\mathbb{C})$ gauge transformation in (20). Note
that in the current work we make the choice of $e_{o}(v)$ such that
$e_{o}(v)\neq e_{o}(v^{\prime})$ for all pairs of $v$ and $v^{\prime}$ sharing
the same edge. Second, for each internal edge $e$, we choose a vertex
$v_{o}(e)$ and parametrize $g_{v_{o}(e)e}$ as

$g_{v_{o}(e)e}=\mathring{g}_{v_{o}(e)e}\begin{pmatrix}\frac{x_{v_{o}(e)e}^{1}}{\sqrt{2}}&\frac{x_{v_{o}(e)e}^{2}+y_{v_{o}(e)e}^{2}}{\sqrt{2}}\\\
0&\mu_{v_{o}(e)e}\end{pmatrix}.$ (72)

Here, $\mu_{v_{o}(e)e}$, along with analogous parameters introduced
subsequently, is determined by $\det(g_{v_{o}(e)e})=1$. This fixes the
$\mathrm{SU(2)}$ gauge transformation in (20). Note that (72) needs to be
applied to all internal edges, including those $e_{o}(v)$ that were chosen for
gauge fixing the $\mathrm{SL}(2,\mathbb{C})$ transformation. For each
$e_{o}(v)\equiv\tilde{e}$, its corresponding $v_{o}(\tilde{e})$ is selected as
the other endpoint compared to $v$. This choice is feasible due to the
conditions we imposed on selecting $v_{o}(e)$. Moreover, all of the other
$g_{ve}$ are parameterized as

$g_{ve}=\mathring{g}_{ve}\begin{pmatrix}1+\frac{x_{ve}^{1}+iy_{ve}^{1}}{\sqrt{2}}&\frac{x_{ve}^{2}+y_{ve}^{2}}{\sqrt{2}}\\\
\frac{x_{ve}^{3}+y_{ve}^{3}}{\sqrt{2}}&\mu_{ve}\end{pmatrix}.$ (73)

Furthermore, the spinors $z_{vf}$ are parametrized as

$z_{vf}=\left\\{\begin{array}[]{ll}(1,\frac{\mathring{z}_{vf}^{2}}{\mathring{z}_{vf}^{1}}+x_{vf}),&\text{
if }\mathring{z}_{vf}^{1}\neq 0,\\\ (x_{vf},1),&\text{ if
}\mathring{z}_{vf}^{1}=0.\end{array}\right.$ (74)

which fixes the rescaling transformation of $z_{vf}$ in (20). Finally $j_{f}$
is parametrized as

$j_{f}=\mathring{j}_{f}(1+x_{f}).$ (75)

Substituting the above parametrization into the action (17), we write the
action as a function of the real variables $x^{i}_{ve}$, $y^{i}_{ve}$,
$x_{vf}$ and $x_{f}$. We then extend the domain of these variables from the
real numbers to complex numbers. The critical equations with respect to these
complex variables are given by $\partial_{z}S=0$ and
$\partial_{\overline{z}}S=0$ for $z$ taking the complexification of
$x^{i}_{ve},y^{i}_{ve},x_{vf}$ or $x_{f}$, where $\overline{z}$ denote the
complex conjugate of $z$.

A technique issue appearing here is the incorporation of the value
$\Xi_{\ell}$ provided in Sec. III.2 into our actual computation. This issue is
noteworthy due to the lack of explanation regarding the definition of the
twist angle $\theta_{\ell}$. According to the expression of the action given
in (II.4), $S$ takes the form

$S=S_{o}+\sum_{f\subset\partial
B}i\mathrm{sgn}(p_{\ell_{f}})j_{f}\phi_{\ell_{f}},$ (76)

where $\phi_{\ell}$ contains the twist angle $\theta_{\ell}$ and the dihedral
angle $\Xi_{\ell}$. We introduce $\mathring{\phi}_{\ell}$ as the coherent
label given by

$\mathring{\phi}_{\ell_{f}}:=\left.i\partial_{j_{f}}S_{o}\right|_{x_{o}},\quad\forall
f\subset\partial B.$ (77)

$\mathring{\phi}_{\ell}$ ensures $\partial_{j_{f_{b}}}S$ takes vanishing value
at $x_{o}$ for all boundary face $f_{b}$. As we discussed in (C.6),
$\mathring{\phi}_{\ell_{f}}$ is equal to the twist angle plus a term
proportional to the dihedral angle. Moreover, the twist angles depend only on
the intrinsic geometry which remains consistent at both the real and the
complex critical points. As a consequence, the twist angle in
$\mathring{\phi}_{\ell}$ is identical to that in $\phi_{\ell}$. Combining this
observation with (44), we arrive at the expression

$\phi_{\ell}=\mathring{\phi}_{\ell}+\beta(\mathring{\Xi}_{\ell}-\mathrm{sgn}(p_{\ell_{f}})\mathrm{sgn}(e)\Xi_{\ell}).$
(78)

where $\mathring{\Xi}_{\ell}=-\partial_{\beta}\mathring{\phi}_{\ell}$ denotes
the dihedral angle compatible with the flat Regge geometry associated with
$x_{o}$. Substituting the result into (77), we get

$\displaystyle S=$ $\displaystyle S_{o}+\sum_{f\subset\partial
B}i\mathrm{sgn}(p_{\ell_{f}})j_{f}\mathring{\phi}_{\ell_{f}}$ (79)
$\displaystyle-i\beta\sum_{f\subset\partial
B}(\mathrm{sgn}(e)\Xi_{\ell}-\mathrm{sgn}(p_{\ell_{f}})\mathring{\Xi}_{\ell})$

For a given real critical point
$x_{o}=(\mathring{j}_{f},\mathring{g}_{ve},\mathring{z}_{vf})$, the parity
transformation

$(\mathring{j}_{f},\mathring{g}_{ve},\mathring{z}_{vf})\mapsto(\mathring{j}_{f},\mathring{g}_{ve}^{\dagger}{}^{-1},\frac{\mathring{g}_{ve}\mathring{g}_{ve}^{\dagger}\mathring{z}_{vf}}{\|\mathring{g}_{ve}^{\dagger}\mathring{z}_{vf}\|^{2}})\equiv\tilde{x}_{o}$
(80)

yields another real critical point denoted by $\tilde{x}_{o}$. As shown in Han
and Zhang (2013), under the parity transformation, $\mathring{\Xi}_{\ell}$ is
transformed as

$\mathring{\Xi}_{\ell}\mapsto-\mathring{\Xi}_{\ell}.$ (81)

In the numerical procedure for searching the complex critical point, the
initial point should be chosen from either $x_{o}$ or $\tilde{x}_{o}$ such
that the corresponding $\mathring{\Xi}_{\ell}$ is closer to $\Xi_{\ell}$.

##  V Numerical investigation

In the current work, we consider the case with the parameter,

$GM=2\times 10^{5}\sqrt{\beta\kappa\hbar},\quad\beta=\frac{1}{10}.$

In this case, the parameter $\alpha$ takes the value

$\alpha=\frac{\sqrt{3}}{50}\beta\kappa\hbar.$

With the chosen parameters, we can proceed with the numerical computation step
by step. Some keys issues in the numerical computation are introduced as
follows.

###  V.1 Fix the boundary $\partial B$

For the numerical calculation, it is necessary to fix the boundary $\partial
B$. As outlined in the beginning of Sec. III, this requires us to set the
parameters $r_{a}$ with $a=1,2,3$ and $\theta_{i}$ with $i=1,2$ for
$\partial_{-}B$. Subsequently, the boundary $\partial_{+}B$ can be determined
accordingly, taking into account the time reversal symmetry.

In our computation, we choose

$\displaystyle r_{1}$ $\displaystyle=\left(\alpha GM/2\right)^{1/3},$ (82)
$\displaystyle r_{2}$ $\displaystyle=(2GM+\frac{GM}{100}),\
r_{3}=(2GM+\frac{GM}{10}).$

Here, $r_{1}$ is chosen as the minimal radial of the collapsing dust ball in
the quantum Oppenheimer-Snyder model. $r_{2}<r_{3}$ are choose as two radials
nearby the horizon.

We select the values of $\theta_{1}$ and $\theta_{2}$ to ensure that all the
tetrahedra in the flat geometry corresponding to the real critical point
$x_{o}$ are spacelike. This selection guarantees the applicability of the
spacelike EPRL model. However, it is important to note that the above
criterion alone cannot uniquely determine the values of $\theta_{1}$ and
$\theta_{2}$. In the current work, we choose the following values as an
example

$\displaystyle\theta_{1}$
$\displaystyle=\theta_{1}^{(c)}-10^{-5}=0.6930807205426226,$ (83)
$\displaystyle\theta_{2}$
$\displaystyle=\theta_{2}^{(c)}-10^{-6}=4.5270907503925455.$

with

$\displaystyle\theta_{1}^{(c)}=$
$\displaystyle\mathrm{arccosh}\left(\frac{9(r_{1}-r_{2})^{2}}{2f(r_{1})r_{2}(2r_{2}-3r_{1})}-1\right),$
$\displaystyle\theta_{2}^{(c)}=$
$\displaystyle\mathrm{arccosh}\left(\frac{9}{4f(r_{3})}-1\right).$

When we embed the region $B_{-}$ into Minkowski spacetime using the map
$\mathfrak{i}$, $\theta_{1}^{(c)}$ and $\theta_{2}^{(c)}$ are the critical
values to ensure the following conditions:

  * (1)

All tetrahedra in the transition surface $\mathfrak{i}(\mathcal{T})$ are
spacelike;

  * (2)

The boundaries $\mathfrak{i}(\partial_{-}B_{1})$ and
$\mathfrak{i}(\partial_{-}B_{2})$ are spacelike;

  * (3)

The time coordinates $t_{2}$ and $t_{3}$ of the spheres $\mathfrak{i}(S_{2})$
and $\mathfrak{i}(S_{3})$ respectively satisfy $t_{2}<t_{3}$, so that
$\mathfrak{i}(S_{3})$ occurs in the future of $\mathfrak{i}(S_{2})$.

Once the boundary $\partial B$ is fixed, the unit $\Delta$ is chosen as the
minimal area in the triangulation of $\partial B$. In our case, it is the area
of the triangle $\mathfrak{p}_{1}\mathfrak{p}_{2}\mathfrak{p}_{3}$. We have

$\Delta=\|\vec{A}_{\mathfrak{p}_{1}\mathfrak{p}_{2}\mathfrak{p}_{3}}\|=264.3604304333458856\beta\kappa\hbar,$
(84)

leading to

$t=\frac{1}{264.3604304333458856},$ (85)

as the parameter for our asymptotic expansion in (66).

###  V.2 Action at the complex critical point

Let us use $S_{\rm res}$ to represent the rescaled action, where the spin foam
action is expressed as $S_{\rm res}/t$. As previously mentioned, $S_{\rm res}$
takes the form

$\displaystyle S_{\rm res}=$ $\displaystyle\sum_{f_{b}\subset\partial
B}\Bigg{\\{}-\frac{1}{2}\left(j_{f_{b}}-\frac{\|\vec{A}_{\ell_{f_{b}}}\|}{\Delta}\right)^{2}+j_{f_{b}}\Big{[}i\mathrm{sgn}(p_{\ell_{f_{b}}})\mathring{\phi}_{\ell_{f_{b}}}-$
(86) $\displaystyle
i\chi\,\beta\Big{(}\mathrm{sgn}(e)\Xi_{\ell}-\mathrm{sgn}(p_{\ell_{f}})\mathring{\Xi}_{\ell}\Big{)}+\mathcal{S}_{o}^{(f_{b})}\Big{]}\Bigg{\\}}+$
$\displaystyle\sum_{\text{internal face
}f_{i}}j_{f_{i}}\mathcal{S}_{o}^{(f_{i})}.$

Here, we introduce the parameter $\chi\in[0,1]$. When $\chi=1$, $S_{\rm res}$
represents the action associated with curved-geometry boundary data.
Conversely, for $\chi=0$, $S_{\rm res}$ corresponds to the action associated
with the flat-geometry boundary data. As illustrated in Fig. 7, we will solve
the critical equations for all $\chi\in[0,1]$, to make the critical point flow
from the real one to the complex one.

By (68), the critical point satisfyies $\partial_{j_{i}}S_{\rm res}=0$ and
$\partial_{j_{b}}S_{\rm res}=0$. In addition, (86) indicates that $S_{\rm
res}$ exhibits linearity concerning $j_{f_{i}}$ and quadratic dependence on
$j_{f_{b}}$. As a consequence, we get the value of $S_{\rm res}$ at the
critical point, denoted by $\mathring{S}_{\rm res}$, as

$\mathring{S}_{\rm res}=\sum_{f_{b}\subset\partial
B}\frac{\mathfrak{b}_{f_{b}}}{2}\left(\frac{2\|\vec{A}_{\ell_{f_{b}}}\|}{\Delta}+\mathfrak{b}_{f_{b}}\right)$
(87)

where $\mathfrak{b}_{f_{b}}$ is defined as

$\displaystyle\mathfrak{b}_{f_{b}}=$ $\displaystyle
i\mathrm{sgn}(p_{\ell_{f_{b}}})\mathring{\phi}_{\ell_{f_{b}}}+\mathcal{S}_{o}^{(f_{b})}\big{|}_{c}$
(88)
$\displaystyle-i\chi\,\beta\Bigg{(}\mathrm{sgn}(e)\Xi_{\ell}-\mathrm{sgn}(p_{\ell_{f}})\mathring{\Xi}_{\ell}\Bigg{)},$

with $\mathcal{S}_{o}^{(f_{b})}\big{|}_{c}$ denoting the value of
$\mathcal{S}_{o}^{(f_{b})}$ at the critical point.

It is worth noting that $\mathcal{S}_{o}^{(f_{b})}\big{|}_{c}$ depends on
$\chi$ since the critical points vary with $\chi$. Indeed, each summand in
(87) is just the value of
$-1/2\left(j_{f_{b}}-\|\vec{A}_{\ell_{f_{b}}}\|/\Delta\right)^{2}+j_{f_{b}}\mathfrak{b}_{f_{b}}$
at its turning point. By substituting (88) into (87), $S_{\rm res}$ at the
critical point reads

$\mathring{S}_{\rm res}=a\chi^{2}+b\chi+c$ (89)

with

$\displaystyle a$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\sum_{f_{b}\subset\partial
B}\left(\phi_{\ell_{f_{b}}}-\mathring{\phi}_{\ell_{f_{b}}}\right)^{2},$ (90)
$\displaystyle b$ $\displaystyle=$ $\displaystyle i\sum_{f_{b}\subset\partial
B}\mathrm{sgn}(p_{\ell_{f_{b}}})\left(\phi_{\ell_{f_{b}}}-\mathring{\phi}_{\ell_{f_{b}}}\right)\times$
(91)
$\displaystyle\left(i\mathrm{sgn}(p_{\ell_{f_{b}}})\mathring{\phi}_{\ell_{f_{b}}}+\mathcal{S}_{o}^{(f_{b})}\big{|}_{c}+\frac{\|\vec{A}_{\ell_{f_{b}}}\|}{\Delta}\right),$
$\displaystyle c$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{f_{b}\subset\partial
B}\left(i\mathrm{sgn}(p_{\ell_{f_{b}}})\mathring{\phi}_{\ell_{f_{b}}}+\mathcal{S}_{o}^{(f_{b})}\big{|}_{c}\right)\times$
(92)
$\displaystyle\left(i\mathrm{sgn}(p_{\ell_{f_{b}}})\mathring{\phi}_{\ell_{f_{b}}}+\mathcal{S}_{o}^{(f_{b})}\big{|}_{c}+\frac{2\|\vec{A}_{\ell_{f_{b}}}\|}{\Delta}\right).$

Here, (78) is applied to get the value of $a$ and $b$. It should be remarked
that the dependence of $\mathcal{S}_{o}^{(f_{b})}\big{|}_{c}$ on $\chi$ leads
to the dependence of $b$ and $c$ on $\chi$. Consequently, (89) is not really a
quadratic function of $\chi$.

###  V.3 Property of the real critical point

As mentioned earlier, we need to find the complex critical point for each set
$\\{g_{\ell}^{(\tilde{s}\bar{s})}\\}_{\ell\subset\partial B}$ with
$\tilde{s},,\bar{s}=\pm 1$. However, the following considerations suggest that
we only need to focus on the case where $\tilde{s}=\bar{s}$.

At first, according to the discussion for (66) and (67), the amplitude
decreases rapidly as the complex critical point moves away from the real
critical point $x_{o}$. The distance between the complex critical point and
$x_{o}$ is measured by the disparity in boundary conditions corresponding to
the complex and real critical points. This observation aligns with the outcome
presented in (89), where the coefficient of the quadratic term is proportional
to the sum of $(\phi_{\ell_{f_{b}}}-\mathring{\phi}_{\ell_{f_{b}}})^{2}$ for
all boundary faces $f_{b}$. Consequently, instead of considering all possible
$\\{g_{\ell}^{(\tilde{s}\bar{s})}\\}_{\ell\subset\partial B}$ with
$\tilde{s},\bar{s}=\pm 1$, we can focus specifically on those cases where the
associated boundary conditions are closer to that of the real critical point
$x_{o}$.

The numerical results demonstrate that $\mathring{\phi}_{\ell}$ defined in
(77) satisfies

$\mathring{\phi}_{\ell}=-\mathring{\phi}_{\mathfrak{r}(\ell)^{-1}}$ (93)

and, for each link $\ell_{c}$ at the corners,

$\mathring{\phi}_{\ell_{c}}=0.$ (94)

As a consequence, among the boundary conditions corresponding to
$\\{g_{\ell}^{(\tilde{s}\bar{s})}\\}_{\ell\subset\partial B}$ with
$\tilde{s},,\bar{s}=\pm 1$, those closer to the boundary conditions given by
$x_{o}$ are the ones satisfying

$\phi_{\ell}=-\phi_{\mathfrak{r}(\ell)^{-1}},\quad\phi_{\ell_{c}}=0.$ (95)

As shown by the equations (49), (51), (53) and (54), the condition (95) holds
only for $\\{g_{\ell}^{(\tilde{s}\bar{s})}\\}_{\ell\subset\partial B}$ with
$\tilde{s}=\bar{s}$.

Indeed, given a link $\ell_{c}$ at the corner, for $\tilde{s}=-\bar{s}$, (60)
yields $|\phi_{\ell_{c}}|=2|\Xi_{\ell_{c}}|$, relating $|\phi_{\ell_{c}}|$ to
the extrinsic curvature at the corner. Given that the corner $S_{1}$ is
located in a highly quantum region, and the choice of $\partial B$ ensures its
local intrinsic flatness, the extrinsic curvature nearby $S_{1}$ is expected
to be very large. That implies a large value of $|\phi_{\ell_{c}}|$ for
$\ell_{c}$ passing through $S_{1}$. Namely, $\phi_{\ell_{c}}$ for $\ell_{c}$
passing through will deviate significantly from
$\mathring{\phi}_{\ell_{c}}=0$. Consequently, this leads to a highly negative
value of the real part of the action at the complex critical point for
$\tilde{s}=-\bar{s}$, as indicated in (67).

###  V.4 Final results

As discussed in the last subsection, our attention should be solely directed
towards the complex critical point for the boundary data
$\\{g_{\ell}^{(++)}\\}_{\ell\subset\partial B}$ and
$\\{g_{\ell}^{(--)}\\}_{\ell\subset\partial B}$. The detailed numerical
results are available in Zhang (2024). We denote the values of
$\mathring{S}_{\rm res}$ associated with these boundary data as
$\mathring{S}_{\rm res}^{(+)}$ and $\mathring{S}_{\rm res}^{(-)}$
respectively. Our numerical investigation reveals that both $\mathring{S}_{\rm
res}^{(+)}$ and $\mathring{S}_{\rm res}^{(-)}$ are real values, as illustrated
in Fig. 8. According to the numerical results, for $\chi=1$, we have

$\displaystyle\mathring{S}_{\rm res}^{(+)}\big{|}_{\chi=1}=$
$\displaystyle-0.0458193513442056,$ (96) $\displaystyle\mathring{S}_{\rm
res}^{(-)}\big{|}_{\chi=1}=$ $\displaystyle-0.0458193513442275.$

It is interesting that there exist a tiny difference between
$\mathring{S}_{\rm res}^{(+)}$ and $\mathring{S}_{\rm res}^{(-)}$. According
to the numerical results in Fig. 8, $\ln(-\mathring{S}_{\rm res}^{(\pm)})$
exhibits an approximately linear dependence on $\ln(\chi)$. The numerical
fitting yields

$\displaystyle\ln(-\mathring{S}_{\rm res}^{(+)})=$ $\displaystyle
1.9999999999996853\ln(\chi)$ (97)
$\displaystyle-3.0830487586650438+\epsilon_{+}(\chi),$
$\displaystyle\ln(-\mathring{S}_{\rm res}^{(-)})=$ $\displaystyle
1.9999999999996856\ln(\chi)$
$\displaystyle-3.0830487586645638+\epsilon_{-}(\chi).$

with the residuals satisfying $|\epsilon_{\pm}(\chi)|<2\times 10^{-12}$.

Figure 8: The values of $-\mathring{S}_{\rm res}^{(+)}$ and
$\delta\mathring{S}_{\rm res}=\mathring{S}_{\rm res}^{(+)}-\mathring{S}_{\rm
res}^{(-)}$ as a function of the parameter $\chi\in[0,1]$.

The nature of the geometry corresponding to the complex critical point is
still an open issue. It is expected to deviate from the geometry corresponding
to the real critical point $x_{o}$ with some perturbations. Let us refer to
the geometry associated with $x_{o}$ as the $x_{o}$-geometry for convenience.
Then, the $x_{o}$-geometry can provide insight into geometry associated with
the complex critical point.

An interesting property of the $x_{o}$-geometry is that the orientation of the
4-simplices in $B_{-}$ is dynamical. By computing the deficit angle
$\vartheta_{\mathfrak{p}_{a}\mathfrak{p}_{b^{\prime}}\mathfrak{p}_{c^{\prime\prime}}}$
hinged at the internal triangle
$\mathfrak{p}_{a}\mathfrak{p}_{b^{\prime}}\mathfrak{p}_{c^{\prime\prime}}\subset
B_{-}$, we obtain

$\displaystyle|\vartheta_{\mathfrak{p}_{a}\mathfrak{p}_{b^{\prime}}\mathfrak{p}_{b^{\prime\prime}}}|=1.6756169215214705\neq
0,$ (98) $\displaystyle\forall a,b=1,2,3,4,\ a\neq b.$

Let $\Theta_{eve^{\prime}}$ be the dihedral angle between the 2 tetrahedra $e$
and $e^{\prime}$ in the 4-simplex $v$. Then, we have

$\vartheta_{\mathfrak{p}_{a}\mathfrak{p}_{b^{\prime}}\mathfrak{p}_{b^{\prime\prime}}}=\sum_{v}\Theta_{eve^{\prime}},$
(99)

where the sum is taken over all 4-simplices sharing
$\mathfrak{p}_{a}\mathfrak{p}_{b^{\prime}}\mathfrak{p}_{b^{\prime\prime}}$.
Equation (98) implies

$\sum_{v}\Theta_{eve^{\prime}}\neq 0.$ (100)

Furthermore, the fact that $x_{o}$ is a real critical point implies

$\sum_{v}\mathrm{sgn}(V_{v})\Theta_{eve^{\prime}}=0,$ (101)

where $V_{v}$ represents the volume of the 4-simplex $v$. Consequently,
$\mathrm{sgn}(V_{v})$ varies among the 4-simplexes sharing the triangles
$\mathfrak{p}_{a}\mathfrak{p}_{b^{\prime}}\mathfrak{p}_{b^{\prime\prime}}$.
For example, we consider the case where $a=2$ and $b=1$. The triangle
$\mathfrak{p}_{2}\mathfrak{p}_{1^{\prime}}\mathfrak{p}_{1^{\prime\prime}}$ is
shared by four 4-simplices

$\displaystyle v_{1}$
$\displaystyle=\mathfrak{p}_{2}\mathfrak{p}_{3}\mathfrak{p}_{4}\mathfrak{p}_{1^{\prime}}\mathfrak{p}_{1^{\prime\prime}},\quad
v_{2}=\mathfrak{p}_{2}\mathfrak{p}_{3}\mathfrak{p}_{1^{\prime}}\mathfrak{p}_{4^{\prime}}\mathfrak{p}_{1^{\prime\prime}},$
(102) $\displaystyle v_{3}$
$\displaystyle=\mathfrak{p}_{2}\mathfrak{p}_{1^{\prime}}\mathfrak{p}_{3^{\prime}}\mathfrak{p}_{4^{\prime}}\mathfrak{p}_{1^{\prime\prime}},\quad
v_{4}=\mathfrak{p}_{2}\mathfrak{p}_{4}\mathfrak{p}_{1^{\prime}}\mathfrak{p}_{3^{\prime}}\mathfrak{p}_{1^{\prime\prime}}.$

We find

$\displaystyle-\mathrm{sgn}(V_{v_{1}})=\mathrm{sgn}(V_{v_{2}})=\mathrm{sgn}(V_{v_{3}})=\mathrm{sgn}(V_{v_{4}}).$
(103)

Indeed, for each triangle
$\mathfrak{p}_{a}\mathfrak{p}_{b^{\prime}}\mathfrak{p}_{b^{\prime\prime}}$,
its boundary always contains a unique 4-simplex $v_{1}$ given by
$v_{1}=\mathfrak{p}_{a}\mathfrak{p}_{c}\mathfrak{p}_{d}\mathfrak{p}_{b^{\prime}}\mathfrak{p}_{b^{\prime\prime}}$,
where $c$ and $d$ are neither $a$ nor $b$. The numerical computation confirms
that the sign of $V_{v_{1}}$ is opposite to the signs of all other 4-simplices
sharing the triangle.

The dihedral angles occurring on the transition surface $\mathcal{T}$ are
associated with the triangles in the triangulation of $\mathcal{T}$. In the
current model, the values of these dihedral angles can be classified into two
types. The first type pertains to the triangles where two points lie on the
sphere $S_{1}$, while the third point lies on the sphere $S_{3}$. An example
is the triangle
$\mathfrak{p}_{1}\mathfrak{p}_{2}\mathfrak{p}_{2^{\prime\prime}}$. The second
type relates to the triangles where only one point is on the sphere $S_{1}$,
and the other two points are on the sphere $S_{3}$. An example is the triangle
$\mathfrak{p}_{1}\mathfrak{p}_{3^{\prime\prime}}\mathfrak{p}_{4^{\prime\prime}}$.
In the current numerical computation, the dihedral angle of the first type
takes the value

$|\Theta_{1}|=5.6445575060432456,$ (104)

and the second type takes

$|\Theta_{2}|=0.4095700383700623.$ (105)

##  VI Summary and outlook

In this work, we investigate the quantum dynamics of the $B$-region in the
black-to-white hole transition suggested in Han  _et al._ (2023a). To
construct the spin foam amplitude, we first introduce a systematic procedure
to yield a bulk triangulation from the boundary triangulation in the
$B$-region. Then, we define the spin foam amplitude based on the 2-complex
dual to the triangulation, where the Thiemann’s complexifier coherent state is
chosen as the boundary state to resemble the semiclassical boundary geometry.
Since the triad fields with different orientations, i.e., $e^{a}_{i}$ and
$-e^{a}_{i}$, yield the same boundary geometry, our work considers the
superposition of the coherent states with respect to both orientations as the
boundary state. The superposition switch on quantum tunneling process
involving the change in orientation.

Moreover, we examine how the amplitude behaves as $t$ approaches $0$, where
$t$ is the parameter describing the variance of the coherent state. Within
this context, the amplitude can be computed numerically with the method of
stationary phase approximation. This method involves solving the critical
equation $\delta S=0$, where $S$ represents the effective action in the
amplitude. Following the procedure outlined in the paper to solve the critical
equation, we obtain our numerical results, as shown in Sec. V.

As discussed in Sec. V.3, the amplitude is dominated by the terms allowing the
change of orientation from the past boundary $\partial_{-}B$ to the future
boundary $\partial_{+}B$. This result indicates that the black-to-white hole
transition should be accompanied by quantum tunneling of the change of
orientation from $\partial_{-}B$ to $\partial_{+}B$. It is noted that the
change of orientation was also discovered in the loop quantum reduced model
Giesel  _et al._ (2023c). It will be interesting to explore how the findings
in the two models relate to each other.

In the amplitude, there are two terms that allow for the change of
orientation. They correspond to the boundary data
$\\{g_{\ell}^{(++)}\\}_{\ell\subset\partial B}$ and
$\\{g_{\ell}^{(--)}\\}_{\ell\subset\partial B}$ as described in Sec. III.3.
According to the numerical results presented in Sec. V.4, the values of the
action at the critical point differ slightly for these two sets of boundary
data. The reasons for the difference are still unclear and require future
investigation.

The current study focuses solely on determining the value of the effective
action at the complex critical point. However, to obtain the final value of
the amplitude, additional factors such as the determinant of the Hessian at
the complex critical point, as indicated in equation (66), are required.
Furthermore, our computations need to fix the values of parameters such as the
black hole mass $M$, the radii $r_{1}$, $r_{2}$, and $r_{3}$, as well as the
angles $\theta_{1}$ and $\theta_{2}$. That is to say, the amplitude is
inherently a function of these parameters. To gain further insights, it is
necessary to compute the amplitude for various sets of parameter values and
analyze their relative significance. It is expected that the relative values
of the amplitude with respect to these parameters will be more physically
meaningful than the absolute amplitude value. Addressing these issues will be
part of our future work.

###### Acknowledgements.

This work is supported by the National Natural Science Foundation of China
with Grants No. 12275022. M.H. receives supports from the NSF through grant
PHY-2207763, the College of Science Research Fellowship at Florida Atlantic
University. MH, DQ and CZ are supported by research grants provided by the
Blaumann Foundation. The Government of Canada supports research at Perimeter
Institute through Industry Canada and by the Province of Ontario through the
Ministry of Economic Development and Innovation. This work benefits from the
visitor’s supports from Beijing Normal University, FAU Erlangen-Nürnberg, and
the University of Western Ontario.

##  Appendix A Spin foam amplitude for boundary face with the Thiemann’s
coherent state

Let $f$ be a boundary face surounded by edges $e_{0},e_{1},\cdots,e_{n}$,
where the order of the edges is compatible with the orientation of $f$. Thus
$e_{0}$ and $e_{n}$ are the two edges connected to the boundary. The vertices
in the boundary of $f$ are labeled as $v_{n}=e_{n-1}\cap e_{n}$. The direction
of boundary link $\ell$ is defined to go from $e_{0}$ to $e_{n}$. The
amplitude is a function of the boundary holonomy $h_{\ell}$, reading

$A_{f}(h_{\ell_{f}})=\tr(Y_{\gamma}^{\dagger}g_{e_{0}v_{1}}g_{v_{1}e_{1}}Y_{\gamma}\cdot
Y_{\gamma}^{\dagger}g_{e_{1}v_{2}}g_{v_{2}e_{2}}Y_{\gamma}\cdot\cdots\cdot
Y_{\gamma}^{\dagger}g_{e_{n-1}v_{n}}g_{v_{n}e_{n}}Y_{\gamma}h_{\ell}^{-1})$
(106)

In this work, we are interested in the case where the boundary state is taken
as the coherent state Thiemann (2001); Thiemann and Winkler (2001)

$\psi_{g_{\ell}}(h_{\ell})=\sum_{j_{\ell}}d_{j_{\ell}}e^{-\frac{t}{2}j_{\ell}(j_{\ell}+1)}\sum_{m}D^{j_{\ell}}_{mm}(g_{\ell}^{-1}h_{\ell}),$
(107)

with $g_{\ell}\in\mathrm{SL}(2,\mathbb{C})$. Then, we have the spin foam
amplitude associated with $f$ as

$\displaystyle\langle
A_{f}|\psi_{g_{\ell_{f}}}\rangle=\int\differential\mu_{H}(h_{\ell})A_{f}(h_{\ell})\psi_{g_{\ell_{f}}}.$
(108)

Following the derivations for an internal face Han and Krajewski (2014), we
finally get

$\langle
A_{f}|\psi_{g_{\ell_{f}}}\rangle=\sum_{j_{f}}\int\prod_{f}(d_{j_{f}})^{|V_{f}|+1}\prod_{vf}\differential^{2}\Omega_{vf}e^{S_{f}}.$
(109)

where the action $S_{f}$ is

$\displaystyle S_{f}=$ $\displaystyle j_{f}\Bigg{[}\sum_{\text{ internal edges
}e}i\beta\ln(\frac{\|Z_{s_{e}ef}\|^{2}}{\|Z_{t_{e}ef}\|^{2}})+\sum_{\text{internal
edges }e}2\ln(\frac{\langle
Z_{s_{e}ef},Z_{t_{e}ef}\rangle}{\|Z_{s_{e}ef}\|\|Z_{t_{e}ef}\|})+$ (110)
$\displaystyle
i\beta\ln(\frac{\|Z_{s_{e_{n}}e_{n}f}\|^{2}}{\|Z_{t_{e_{0}}e_{0}f}\|^{2}})+2\ln(\frac{\langle
Z_{s_{e_{n}}e_{n}f},g_{\ell}^{-1}Z_{t_{e_{0}}e_{0}f}\rangle}{\|Z_{s_{e_{n}}e_{n}f}\|\|Z_{t_{e_{0}}e_{0}f}\|})\Bigg{]}-\frac{t}{2}j_{f}(j_{f}+1)$

with $Z_{vef}=g_{ve}^{\dagger}z_{vf}$ and $\|Z\|=\sqrt{\langle Z,Z\rangle}$
for $Z\in\mathbb{C}^{2}$.

##  Appendix B Equation of motion obtained from the spin foam action

In the spin foam model, the amplitude reads

$A=\sum_{j}\int\differential g\differential ze^{S(j,g,z)}.$ (111)

The action $S$, in the EPRL model, takes the form (see, e.g., Han and
Krajewski (2014); Han  _et al._ (2020) and Appendix A)

$S(j,g,z)=\sum_{\text{ internal face }f}j_{f}S_{f}+\sum_{\text{ boundary face
}f}\left(j_{f}S^{\prime}_{f}-\frac{t}{2}j_{f}(j_{f}+1)\right)$ (112)

where $S_{f}$ and $S^{\prime}_{f}$ are functions depending merely on the
$\mathrm{SL}(2,\mathbb{C})$ element $g$ and the spinor $z$.

Applying the Poisson summation formula for the amplitude Han  _et al._ (2022),
one has

$A=\sum_{k=-\infty}^{\infty}\int\differential g\differential z\differential
x^{\prime}e^{S(x^{\prime},g,z)-2\pi ikx^{\prime}}\cong\int\differential
g\differential z\differential x^{\prime}e^{S(x^{\prime},g,z)},$ (113)

where the focus is on the phase in the exponential:

$S=\sum_{\text{ internal face }f}x^{\prime}_{f}S_{f}+\sum_{\text{ boundary
face
}f}\left(x^{\prime}_{f}S^{\prime}_{f}-\frac{t}{2}x^{\prime}_{f}(x^{\prime}_{f}+1)\right).$
(114)

To study the asymptotic behavior of $A$ for $t\to 0$, we rescale $x^{\prime}$
by introducing $x=tx^{\prime}$ so that we get

$S=\frac{1}{t}\left(\sum_{\text{ internal face }f}x_{f}S_{f}+\sum_{\text{
boundary face
}f}\left(x_{f}S^{\prime}_{f}-\frac{1}{2}x_{f}^{2}\right)\right)-\frac{x_{f}}{2}.$
(115)

Then, the asymptotic behavior for $t\to 0$ can be analyzed by the stationary
phase approximation, where we define the effective action as

$S_{\rm eff}=\sum_{\text{ internal face }f}x_{f}S_{f}+\sum_{\text{ boundary
face }f}\left(x_{f}S^{\prime}_{f}-\frac{1}{2}x_{f}^{2}\right)$ (116)

The critical point is the solution to the equation of motion $\delta S_{\rm
eff}=0$. The properties of the critical point have been well studied in cases
where the 2-complex contains no boundary (see, e.g., Han and Zhang (2013); Han
and Krajewski (2014)). We can apply the results obtained from these studies,
with the only exception being the treatment of the boundary faces carrying the
coherent state (107). Therefore, our attention should be directed towards the
effective action associated with a boundary face $f$

$\displaystyle S_{\rm eff}^{(f)}=$ $\displaystyle j_{f}\Bigg{[}\sum_{\text{
internal edges
}e}i\beta\ln(\frac{\|Z_{s_{e}ef}\|^{2}}{\|Z_{t_{e}ef}\|^{2}})+\sum_{\text{internal
edges }e}2\ln(\frac{\langle
Z_{s_{e}ef},Z_{t_{e}ef}\rangle}{\|Z_{s_{e}ef}\|\|Z_{t_{e}ef}\|})+$ (117)
$\displaystyle
i\beta\ln(\frac{\|Z_{s_{e_{n}}e_{n}f}\|^{2}}{\|Z_{t_{e_{0}}e_{0}f}\|^{2}})+2\ln(\frac{\langle
Z_{s_{e_{n}}e_{n}f},g_{\ell}^{-1}Z_{t_{e_{0}}e_{0}f}\rangle}{\|Z_{s_{e_{n}}e_{n}f}\|\|Z_{t_{e_{0}}e_{0}f}\|})-\frac{t}{2}j_{f}^{2}\Bigg{]}.$

Then, (116) is a result of the replacement $j_{f}\mapsto x_{f}/t$.

Treating $j_{f}$ as an integration variable, we need to study the equation
$\partial_{j_{f}}S=0$, indicating that both its real and imaginary part are 0.
Starting with its real part, we get

$\displaystyle\sum_{\text{internal edges }e}\mathrm{Re}\ln(\frac{\langle
Z_{s_{e}ef},Z_{t_{e}ef}\rangle}{\|Z_{s_{e}ef}\|\|Z_{t_{e}ef}\|})+\mathrm{Re}\ln(\frac{\langle
Z_{s_{e_{n}}e_{n}f},g_{\ell}^{-1}Z_{t_{e_{0}}e_{0}f}\rangle}{\|Z_{s_{e_{n}}e_{n}f}\|\|Z_{t_{e_{0}}e_{0}f}\|})-\frac{t}{2}j_{f}=0$
(118)

We impose the following decomposition for $g_{\ell}$

$g_{\ell}=n_{s_{\ell}}e^{(-ip+\phi)\tau_{3}}n_{t_{\ell}}^{-1},$ (119)

where $n_{e_{\ell}},n_{t_{\ell}}\in\mathrm{SU(2)}$, $p,\phi\in\mathbb{R}$.
Then, we could define new variables

$\tilde{Z}_{s_{e_{n}}e_{n}f}=n_{t_{\ell}}^{-1}Z_{s_{e_{n}}e_{n}f},\quad\tilde{Z}_{t_{e_{0}}e_{0}f}=n_{s_{\ell}}^{-1}Z_{t_{e_{0}}e_{0}f}$
(120)

and rewrite (118) as

$\displaystyle\sum_{\text{internal edges }e}\mathrm{Re}\ln(\frac{\langle
Z_{s_{e}ef},Z_{t_{e}ef}\rangle}{\|Z_{s_{e}ef}\|\|Z_{t_{e}ef}\|})+\mathrm{Re}\ln(\frac{\langle\tilde{Z}_{s_{e_{n}}e_{n}f},e^{(ip-\phi)\tau_{3}}\tilde{Z}_{t_{e_{0}}e_{0}f}\rangle}{\|\tilde{Z}_{s_{e_{n}}e_{n}f}\|\|\tilde{Z}_{t_{e_{0}}e_{0}f}\|})-\frac{t}{2}j_{f}=0$
(121)

To solve this equation, we impose an additional requirement that the real part
of the action itself should vanish. The additional requirement is motivated by
the desire to find the critical point where the transition amplitude takes its
maximal absolute value. However, a subtle complication arises due to the non-
normalization of the employed boundary state. The norm of the boundary state
$\psi_{g}^{t}$ takes Zhang  _et al._ (2020b, 2022b)

$\|\psi_{g}^{t}\|=\frac{\sqrt{2\sqrt{\pi}}e^{t/8}\sqrt{|p|}e^{\frac{|p|^{2}}{2t}}}{t^{3/4}\sqrt{\sinh(|p|)}}.$
(122)

Therefore, for effective action corresponding to the transition amplitude
$\langle A_{f}|\psi_{g}^{t}\rangle/\|\psi_{g}^{t}\|$, the condition of
vanishing its real part results in

$\displaystyle\sum_{\text{internal edges }e}\mathrm{Re}\ln(\frac{\langle
Z_{s_{e}ef},Z_{t_{e}ef}\rangle}{\|Z_{s_{e}ef},\|\|Z_{t_{e}ef}\|})+\mathrm{Re}\ln(\frac{\langle\tilde{Z}_{s_{e_{n}}e_{n}f},e^{(ip-\phi)\tau_{3}}\tilde{Z}_{t_{e_{0}}e_{0}f}\rangle}{\|\tilde{Z}_{s_{e_{n}}e_{n}f}\|\|\tilde{Z}_{t_{e_{0}}e_{0}f}\|})$
(123) $\displaystyle-\frac{|p|^{2}}{4j_{f}t}-\frac{t}{4}j_{f}=0$

Comparing this equation with (121), we get

$-\frac{|p|^{2}}{4j_{f}t}-\frac{t}{4}j_{f}=-\frac{t}{2}j_{f},$ (124)

leading to the area matching condition

$j_{f}=\frac{|p|}{t}.$ (125)

Taking the area matching condition into account, one can demonstrate

$\displaystyle
e^{-tj_{f}}\left|\left\langle\frac{\tilde{Z}_{s_{e_{n}}e_{n}f}}{\|\tilde{Z}_{s_{e_{n}}e_{n}f}\|},e^{(ip-\phi)\tau_{3}}\frac{\tilde{Z}_{s_{e_{0}}e_{0}f}}{\|\tilde{Z}_{s_{e_{0}}e_{0}f}\|}\right\rangle^{2}\right|$
(126) $\displaystyle=$
$\displaystyle\left(\sin\left(\frac{\beta_{0}}{2}\right)\sin\left(\frac{\beta_{n}}{2}\right)e^{\frac{p-|p|}{2}}+\cos\left(\frac{\beta_{0}}{2}\right)\cos\left(\frac{\beta_{n}}{2}\right)e^{\frac{-|p|-p}{2}}\right)^{2}-e^{-|p|}\frac{1}{2}\sin(\beta_{0})\sin(\beta_{n})[1-\cos(\gamma_{0}-\gamma_{n}-\phi)]$
$\displaystyle\leq$ $\displaystyle\cos^{2}(\frac{\beta_{n}-\beta_{0}}{2})\leq
1.$

Here, the parametrization is applied as follows

$\displaystyle\frac{\tilde{Z}_{s_{e_{0}}e_{0}f}}{\|\tilde{Z}_{s_{e_{0}}e_{0}f}\|}=$
$\displaystyle\left\\{\sin\left(\frac{\beta_{0}}{2}\right)\exp\left(\frac{i\alpha_{0}}{2}-\frac{i\gamma_{0}}{2}\right),\cos\left(\frac{\beta_{0}}{2}\right)\exp\left(\frac{i\alpha_{0}}{2}+\frac{i\gamma_{0}}{2}\right)\right\\},$
(127)
$\displaystyle\frac{\tilde{Z}_{s_{e_{n}}e_{n}f}}{\|\tilde{Z}_{s_{e_{n}}e_{n}f}\|}=$
$\displaystyle\left\\{\sin\left(\frac{\beta_{n}}{2}\right)\exp\left(\frac{i\alpha_{n}}{2}-\frac{i\gamma_{n}}{2}\right),\cos\left(\frac{\beta_{n}}{2}\right)\exp\left(\frac{i\alpha_{n}}{2}+\frac{i\gamma_{n}}{2}\right)\right\\}.$

The result (126), combined with the inequality $\frac{\langle
Z_{s_{e}ef},Z_{t_{e}ef}\rangle^{2}}{\|Z_{s_{e}ef}\|^{2}\|Z_{t_{e}ef}\|^{2}}\leq
1$, results in the solution of (121),

$\displaystyle\frac{Z_{s_{e}ef}}{\|Z_{s_{e}ef}\|}=$ $\displaystyle
e^{i\alpha_{ef}}\frac{Z_{t_{e}ef}}{\|Z_{t_{e}ef}\|},$ (128)
$\displaystyle\frac{\tilde{Z}_{t_{e_{0}}s_{0}f}}{\|\tilde{Z}_{t_{e_{0}}s_{0}f}\|}=$
$\displaystyle\left\\{\begin{array}[]{ll}e^{i\frac{\alpha_{0}-\gamma_{0}}{2}}(1,0),&p>0,\\\
e^{i\frac{\alpha_{0}+\gamma_{0}}{2}}(0,1),&p<0,\end{array}\right.$
$\displaystyle\frac{\tilde{Z}_{s_{e_{n}}s_{n}f}}{\tilde{Z}_{s_{e_{n}}s_{n}f}}=$
$\displaystyle\left\\{\begin{array}[]{ll}e^{i\frac{\alpha_{n}-\gamma_{n}}{2}}(1,0),&p>0,\\\
e^{i\frac{\alpha_{n}+\gamma_{n}}{2}}(0,1),&p<0.\end{array}\right.$

When (128) is satisfied, we have

$\displaystyle\mathfrak{P}\left(\left\langle\frac{\tilde{Z}_{s_{e_{n}}e_{n}f}}{\|\tilde{Z}_{s_{e_{n}}e_{n}f}\|},e^{(ip-\phi)\tau_{3}}\frac{\tilde{Z}_{s_{e_{0}}e_{0}f}}{\|\tilde{Z}_{s_{e_{0}}e_{0}f}\|}\right\rangle^{2}\right)\Bigg{|}_{\text{critical
point
}}=e^{i(\alpha_{0}-\alpha_{n})-\mathrm{sgn}(p)(\gamma_{0}-\gamma_{n}-\phi)}$
(129)

where $\mathfrak{P}(z):=z/|z|$ for $z\in\mathbb{C}$. In addition, combining
(120) and (128), we have

$\frac{Z_{t_{e_{0}}e_{0}f}}{\|Z_{t_{e_{0}}e_{0}f}\|}=n_{s_{\ell}}e^{i\frac{\alpha_{0}-\mathrm{sgn}(p)\gamma_{0}}{2}}\mathring{Z}(p),\quad\frac{Z_{s_{e_{n}}e_{n}f}}{\|Z_{s_{e_{n}}e_{n}f}\|}=n_{t_{\ell}}e^{i\frac{\alpha_{n}-\mathrm{sgn}(p)\gamma_{n}}{2}}\mathring{Z}(p).$
(130)

where $\mathring{Z}(p)=(1,0)^{T}$ for $p>0$ and $\mathring{Z}(p)=(0,1)^{T}$
for $p<0$.

The above discussion solves the equation $\partial_{j_{f}}S_{\rm eff}=0$ for a
boundary face $f$. Let us summarize the results as follows:

For a boundary face $f$ carrying the coherent state (107) labelled by
$g_{\ell}=n_{s_{\ell}}e^{(-ip+\phi)\tau_{3}}n_{t_{\ell}}^{-1},$ the area
matching condition $j_{f}=|p|/t$ leads to the following solution of
$\partial_{j_{f}}S_{\rm eff}=0$,

$\displaystyle\frac{Z_{s_{e}ef}}{\|Z_{s_{e}ef}\|}=e^{i\alpha_{ef}}\frac{Z_{t_{e}ef}}{\|Z_{t_{e}ef}\|},\quad\frac{Z_{t_{e_{0}}e_{0}f}}{\|Z_{t_{e_{0}}e_{0}f}\|}=e^{i\alpha_{e_{0}f}}\xi_{s_{\ell}},\quad\frac{Z_{s_{e_{n}}e_{n}f}}{\|Z_{s_{e_{n}}e_{n}f}\|}=e^{i\alpha_{e_{n}f}}\xi_{t_{\ell}},\text{
and, }$ (131) $\displaystyle\sum_{\text{ internal }e\in\partial
f}\beta\ln(\frac{\|Z_{s_{e}ef}\|}{\|Z_{t_{e}ef}\|})+\sum_{\text{internal edges
}e}\mathrm{Im}\ln(\frac{\langle
Z_{s_{e}ef},Z_{t_{e}ef}\rangle}{\|Z_{s_{e}ef}\|\|Z_{t_{e}ef}\|})+$
$\displaystyle\beta\ln(\frac{\|Z_{s_{e_{n}}e_{n}f}\|}{\|Z_{t_{e_{0}}e_{0}f}\|})+\mathrm{Im}\ln(\frac{\langle\xi_{s_{\ell}},Z_{t_{e_{0}}e_{0}f}\rangle}{\|\xi_{s_{\ell}}\|\|Z_{t_{e_{0}}e_{0}f}\|})+\mathrm{Im}\ln(\frac{\langle
Z_{t_{e_{n}}e_{n}f},\xi_{t_{\ell}}\rangle}{\|Z_{t_{e_{n}}e_{n}f}\|\|\xi_{t_{\ell}}\|})+\frac{\mathrm{sgn}(p)}{2}\phi=0,$

for $\alpha_{ef}\in\mathbb{R}$ and $k\in\mathbb{Z}$, where we define

$\xi_{s_{\ell}}=n_{s_{\ell}}\mathring{Z}(p),\
\xi_{t_{\ell}}=n_{t_{\ell}}\mathring{Z}(p).$ (132)

The solutions to the equations $\partial_{z}S=0$ and $\partial_{g}S=0$ for
boundary faces exhibit same properties as the solutions of the same equations
associated with internal faces, allowing for the direct application of results
from previous literature Han and Zhang (2013); Han and Krajewski (2014).

It is worth noting that the first three equations in (131) can be recovered
through the usuall model, where the boundary state is selected as
$|j_{f}\xi_{t_{\ell}}\rangle\otimes\langle j_{f}\xi_{s_{\ell}}|$. Let
$S^{(0)}_{f}$ denote the value of the action taken at the critical points in
that conventional model. Then, the last equation in (131) can be expressed as

$S^{(0)}_{f}+ij_{f}\mathrm{sgn}(p)\phi=0.$ (133)

This outcome indicates that an alternative choice for the boundary state is
the Gaussian state given by

$\psi_{\ell}=\sum_{j_{f}}\exp(-\frac{t}{2}\left(j_{f}-\frac{|p|}{t}\right)^{2}+ij_{f}\mathrm{sgn}(p)\phi)|j_{f}\xi_{t_{\ell}}\rangle\otimes\langle
j_{f}\xi_{s_{\ell}}|,$ (134)

allowing one to derive the same critical equation as in (131). The spin foam
action associated with the boundary face $f$ carrying the boundary state (134)
reads

$\displaystyle S_{f}=$ $\displaystyle j_{f}\Bigg{[}\sum_{\text{ internal
}e\in\partial
f}i\beta\ln(\frac{\|Z_{s_{e}ef}\|^{2}}{\|Z_{t_{e}ef}\|^{2}})+\sum_{\text{internal
edges }e}2\,\ln(\frac{\langle
Z_{s_{e}ef},Z_{t_{e}ef}\rangle}{\|Z_{s_{e}ef}\|\|Z_{t_{e}ef}\|})+$ (135)
$\displaystyle
i\beta\ln(\frac{\|\xi_{s_{\ell}}\|^{2}}{\|Z_{t_{e_{0}}e_{0}f}\|^{2}})+i\beta\ln(\frac{\|Z_{s_{e_{n}}e_{n}f}\|^{2}}{\|\xi_{t_{\ell}}\|^{2}})+2\ln(\frac{\langle\xi_{s_{\ell}},Z_{t_{e_{0}}e_{0}f}\rangle}{\|\xi_{s_{\ell}}\|\|Z_{t_{e_{0}}e_{0}f}\|})$
$\displaystyle+2\ln(\frac{\langle
Z_{t_{e_{n}}e_{n}f},\xi_{t_{\ell}}\rangle}{\|Z_{t_{e_{n}}e_{n}f}\|\|\xi_{t_{\ell}}\|})\Bigg{]}-\frac{t}{2}(j_{f}-\frac{|p|}{t})^{2}+i\,\mathrm{sgn}(p)j_{f}\phi$
$\displaystyle=$ $\displaystyle j_{f}\Bigg{[}\sum_{e\subset\partial
f}i\beta\ln(\frac{\|Z_{s_{e}ef}\|^{2}}{\|Z_{t_{e}ef}\|^{2}})+\sum_{e\subset\partial
f}2\ln(\frac{\langle
Z_{s_{e}ef},Z_{t_{e}ef}\rangle}{\|Z_{s_{e}ef}\|\|Z_{t_{e}ef}\|})\Bigg{]}-\frac{t}{2}\left(j_{f}-\frac{|p|}{t}\right)^{2}+i\,\mathrm{sgn}(p)j_{f}\phi$

where in the last line we introduce the following convention: when $e=e_{0}$,
the term $Z_{s_{e}ef}$ becomes $\xi_{s_{\ell}}$, and when $e=e_{n}$, the term
$Z_{t_{e}ef}$ becomes $\xi_{t_{\ell}}$.

##  Appendix C Bridge Continuous Geometry with Coherent States

###  C.1 Classical phase space

As shown in Han and Zhang (2013), the value of the spin foam action at the
critical point takes the form

$S=-i\epsilon\sum_{f}\sum_{v\text{ surrounding
}f}\mathrm{sgn}(V_{4}(v))j_{f}\Theta_{v}+\text{extra term}$ (136)

where $\epsilon=\pm$ is a global sign on the two complex,
$\mathrm{sgn}(V_{4}(v))$ denote the sign of the 4-volume of the simplex dual
to the vertex $v$, and $\Theta_{v}$ is the 4-D dihedral angle between the two
tetrahedra dual to the boundary edges $e$ and $e^{\prime}$ of $f$ with $e\cap
e^{\prime}=v$. The result is different from the Regge action by the factor
$\mathrm{sgn}(V_{4}(v))$. As we know, the Regge action can be regarded as a
discretized version of the Einstein Hilbert action
$\int\sqrt{-g}R\differential^{4}x$. The sign $\mathrm{sgn}(V_{4}(v))$ factor
in (136) makes us believe that $\sqrt{-g}$ should be replaced by the signed
volume $\mathfrak{e}\equiv\det({\mathfrak{e}}_{\mu}^{I})$ in the classical
action corresponding to the spin foam model, where ${\mathfrak{e}}_{\mu}^{I}$
is the coframe field. Due to the extra term in (136), we consider the action
with the Holst term,

$\displaystyle
S_{G}=\frac{1}{2\kappa\beta}\int_{\mathcal{M}}{\mathfrak{e}}\,{\mathfrak{e}}^{\sigma}_{M}{\mathfrak{e}}^{\theta}_{N}P^{MN}{}_{KL}\Omega_{\sigma\theta}^{KL}\differential^{4}x$
(137)

with $\kappa=8\pi G$ and $\beta$ being the Barbero–Immirzi parameter, where
$P^{MN}{}_{KL}=\beta\delta^{[M}_{K}\delta^{N]}_{L}+\frac{1}{2}\epsilon^{MN}{}_{KL}$,
and $\Omega^{KL}_{\sigma\theta}$ is the curvature form of the
$\mathfrak{sl}(2,\mathbb{C})$ connection compatible with
${\mathfrak{e}}_{\mu}^{I}$. Note that the integral in (137) depends on the
orientation of the coordinate system, unlike the case where the action takes
$\sqrt{|g|}$ as the volume form. This needs us to fix the orientation of the
coordinate system to be the right hand when defining the integral in (137).
Then a minus sign should be added if a left hand coordinate system is chosen
to do the integration.

Let $\mathcal{M}$ be diffeomorphic to $\mathbb{R}\times\Sigma$ with $\Sigma$
dnoting the spatial manifold. Performing the canonical analysis on the action
(137), we get the canonical pairs $(A_{a}^{i},E^{b}_{j})$, where $A_{a}^{i}$
is the usual Ashtekar connection but $E^{b}_{j}$ becomes

$E^{b}_{j}=ee^{b}_{j},$ (138)

with $e\equiv\det(e_{a}^{i})$. In this expression and what follows, indices
$a,b,\cdots$ denotes spatial tensors, and $i,j,k,\cdots$ represent tensors
with value in $\mathfrak{su}(2)$. Under the coordinate transformation from $x$
to $y$, the flux is transformed as

$E^{b}_{j}\to\det(\frac{\partial(x^{1},x^{2},x^{3})}{\partial(y^{1},y^{2},y^{3})})E^{b}_{j}.$
(139)

As a consequence, the dual of the flux, i.e., $E^{b}_{j}\varepsilon_{bcd}$, is
a 2-form. This is in contrast to the case where the flux is
$\sqrt{q}e^{b}_{j}=|e|e^{b}_{j}$, resulting in its dual to be a pseudo 2-form.

The Poisson bracket between the canonical pairs is

$\\{A_{a}^{i}(x),E^{b}_{j}(y)\\}=\kappa\beta\delta_{a}^{b}\delta_{j}^{i}\delta^{3}(x,y).$
(140)

Due to the definition (138) of $E$, the $\delta$-function in the equation is
the one such that $\int\delta^{3}(y_{0},y)\differential^{3}y=1$ for a right
hand coordinate system $y^{a}$, and
$\int\delta^{3}(y_{0},y)\differential^{3}y=-1$ for a left hand coordinate
system $y^{a}$. This ensures the independence of the choice of coordinate
systems for the equality $\\{A_{a}^{i}(x),\int
E^{b}_{j}(y)f^{b}_{j}(y)\differential^{3}y\\}=\kappa\beta f^{b}_{j}(x)$, where
$f^{b}_{j}$ is a smeared function.

###  C.2 From the continuous to the discrete phase spaces

The holonomy $h_{\ell}$ along an oriented link $\ell:t\mapsto\ell(t)$ with
$t\in[0,1]$ is the solution to

$\partial_{t}h_{\ell(t)}(A)=h_{\ell(t)}A_{a}(\ell(t))\dot{\ell}^{a}(t),\
h_{\ell(0)}=1,\ h_{\ell}\equiv h_{\ell(1)}.$ (141)

For a link $\ell$, let $S_{\ell}$ be a 2-surface transversally intersecting
$\ell$ at at $p=\ell(1/2)$. The covariant fluxes are defined as

$\displaystyle p^{i}_{s,\ell}(A,E):=$
$\displaystyle\frac{-1}{\beta\kappa\hbar}\tr\left[-2\tau^{i}\int_{S_{\ell}}h(\rho^{s}_{\ell}(\sigma))E^{c}(\sigma)n_{c}(\sigma)h(\rho^{s}_{\ell}(\sigma)^{-1})\differential^{2}\sigma\right]$
(142) $\displaystyle p^{i}_{t,\ell}(A,E):=$
$\displaystyle\frac{1}{\beta\kappa\hbar}\tr\left[-2\tau^{i}\int_{S_{\ell}}h(\rho^{t}_{\ell}(\sigma))E^{c}(\sigma)n_{c}(\sigma)h(\rho^{t}_{\ell}(\sigma)^{-1})\differential^{2}\sigma\right]$

where $\rho_{\ell}^{s}(x)$ is a path starting from the source of $\ell$,
traveling along $\ell$ until reaching $p$, and then running in $S_{\ell}$ to
reach $x\in S_{\ell}$, and $\rho_{\ell}^{t}(x)$ is a similar path but starting
from the target of $\ell$. In order to provide a clearer understanding of the
definition in (142), let us clarify how to do the integration. To begin with,
a right hand coordinate system $\\{x^{0},\sigma^{1},\sigma^{2}\\}$ is chosen
in a region containing the surface $S_{\ell}$. It is required that
$(\sigma^{1},\sigma^{2})$ could be chosen as a coordinate on $S_{\ell}$. With
this coordinate, $n_{c}$ is defined as the conormal with direction compatible
with the orientation of $S_{\ell}$, satisfying $n_{c}(\partial/\partial
x^{0})^{c}=\pm 1$. Without loss of generality, the orientation of $S_{\ell}$
is chosen to be compatible with $\ell$ in what follows. In the general case
where the orientation of $S_{\ell}$ is allowed to be incompatible with $\ell$,
one only needs the replacements
$p^{i}_{s,\ell}\mapsto\mathrm{sgn}(\dot{\ell}^{a}n_{a})p^{i}_{s,\ell}$ and
$p^{i}_{t,\ell}\mapsto\mathrm{sgn}(\dot{\ell}^{a}n_{a})p^{i}_{t,\ell}$ to
ensure the validity of all subsequent conclusions.

According to (142), one has

$p^{i}_{s,\ell}\tau_{i}=-h_{\ell}p^{i}_{t,\ell}\tau_{i}h_{\ell}^{-1}.$ (143)

The non-vanishing Poisson brackets between the covariant fluxes and the
holonomy are

$\displaystyle\\{h_{\ell},p^{i}_{s,\ell}\\}$
$\displaystyle=-\tau^{i}h_{\ell},$ (144)
$\displaystyle\\{h_{\ell},p^{i}_{t,\ell}\\}$ $\displaystyle=h_{\ell}\tau^{i},$
$\displaystyle\\{p^{i}_{s,\ell},p^{j}_{s,\ell}\\}$
$\displaystyle=-\epsilon_{ijl}p^{l}_{s,\ell},$
$\displaystyle\\{p^{i}_{t,\ell},p^{j}_{t,\ell}\\}$
$\displaystyle=-\epsilon_{ijl}p^{l}_{t,\ell},$

where $\tau_{k}:=-i\sigma_{k}/2$ with $\sigma_{k}$ being the Pauli matrix. In
the quantum theory, a wave function $\psi(A)$ is given by a cylindrical
function

$\psi(A)=\psi_{\gamma}(\\{h_{\ell}(A)\\}_{\ell\in E(\gamma)})$ (145)

where $\gamma$ is a graph embedded in $\Sigma$ comprising links $\ell\in
E(\gamma)$ and $\psi_{\gamma}$ is a complex function on
$\mathrm{SU(2)}^{|E(\gamma)|}$ with $|E(\gamma)|$ denoting the number of
elements in $E(\gamma)$. The action of the flux operators on the wave
functions are

$\displaystyle\hat{p}^{i}_{s,\ell}\psi(A)$
$\displaystyle=i\frac{\differential}{\differential
t}\Big{|}_{t=0}\psi_{\gamma}(h_{\ell^{\prime}}(A),\cdots,e^{t\tau_{i}}h_{\ell}(A),\cdots),$
(146) $\displaystyle\hat{p}^{i}_{t,\ell}\psi(A)$
$\displaystyle=i\frac{\differential}{\differential
t}\Big{|}_{t=0}\psi(h_{\ell^{\prime}}(A),\cdots,h_{\ell}(A)e^{-t\tau_{i}},\cdots),$

which implements the Poisson algebra (144) in such a way that
$[\cdot,\cdot]=i\\{\cdot,\cdot\\}$.

###  C.3 Coherent state

In this section, we focus on a graph of a single link $\ell$ for our
calculation. Given $g_{\ell}\in\mathrm{SL}(2,\mathbb{C})$ which is
parametrized by

$g_{\ell}=e^{-i\vec{p}_{\ell}\cdot\vec{\tau}}u_{\ell},$ (147)

the coherent state labeled by $g_{\ell}$, as shown in (107), is given by

$\psi_{g_{\ell}}^{t}(h_{\ell})=\sum_{j}d_{j}e^{-\frac{t}{2}j(j+1)}\sum_{m}D^{j}_{mm}(h_{\ell}g_{\ell}^{-1}).$
(148)

For a unit vector $\vec{v}=(\cos\alpha\sin\eta,\sin\alpha\sin\eta,\cos\eta),$
we introduce $n(\vec{v})\in\mathrm{SU(2)}$ such that

$n(\vec{v})\tau_{3}n(\vec{v})^{-1}=\vec{v}\cdot\vec{\tau}.$ (149)

Using these notions, we can decompose $g_{\ell}$ as

$g_{\ell}=n(\vec{p}_{\ell}/p_{\ell})e^{(-ip_{\ell}+\phi_{\ell})\tau_{3}}n(\vec{p^{\prime}}_{\ell}/p_{\ell})^{-1},$
(150)

where $p_{\ell}\in\mathbb{R}$ is given by $|p_{\ell}|=\|\vec{p}_{\ell}\|$, and
$\vec{p^{\prime}}_{\ell}$ as well as $\phi_{\ell}$ is defined by

$u_{\ell}=n(\vec{p}_{\ell}/p_{\ell})e^{\phi_{\ell}\tau_{3}}n(\vec{p^{\prime}}_{\ell}/p_{\ell})^{-1}.$
(151)

It can be verified that

$\displaystyle
u_{\ell}^{-1}(\vec{p}_{\ell}\cdot\vec{\tau})u_{\ell}=\vec{p^{\prime}}_{\ell}\cdot\vec{\tau}.$
(152)

The geometry interpretation for the parameters $\vec{p}_{\ell}$,
$\vec{p^{\prime}}_{\ell}$ and $\phi_{\ell}$ can be deduced from the
expectation values of operators. We have Zhang  _et al._ (2020a)

$\displaystyle\frac{\langle\psi_{g_{\ell}}^{t}|\hat{p}_{s,\ell}^{k}|\psi_{g_{\ell}}^{t}\rangle}{\|\psi_{g_{\ell}}^{t}\|^{2}}$
$\displaystyle\cong\frac{1}{t}p_{\ell}^{k}$ (153)
$\displaystyle\frac{\langle\psi_{g_{\ell}}^{t}|\hat{p}_{t,\ell}^{k}|\psi_{g_{\ell}}^{t}\rangle}{\|\psi_{g_{\ell}}^{t}\|^{2}}$
$\displaystyle\cong-\frac{1}{t}p^{\prime}_{\ell}{}^{k},$
$\displaystyle\frac{\langle\psi_{g_{\ell}}^{t}|h_{\ell}|\psi_{g_{\ell}}^{t}\rangle}{\|\psi_{g_{\ell}}^{t}\|^{2}}$
$\displaystyle\cong
n(\vec{p}_{\ell}/p_{\ell})e^{\phi_{\ell}\tau_{3}}n(\vec{p^{\prime}}_{\ell}/p_{\ell})^{-1}.$

In terms of the expection values, $g_{\ell}$ is expressed as

$g_{\ell}=\exp(-it\langle\hat{p}_{s,\ell}^{k}\rangle\tau_{k})\langle
h_{\ell}\rangle.$ (154)

###  C.4 Compute $g_{\ell}$ from the intrinsic and extrinsic geometry

The Ashtekar connection can be written in terms of the spin connection
$\Gamma_{a}^{i}$ and the extrinsic curvature for $K_{a}^{i}$ as

$A_{a}=\Gamma_{a}+\beta K_{a}.$ (155)

Given a link $\ell:[0,1]\ni t\mapsto\ell_{t}\in\Sigma$, we consider
$G_{\ell_{t}}$ as the parallel transport of the spin connection, i.e.,
solution to

$\partial_{t}G_{\ell_{t}}=G_{\ell_{t}}\Gamma_{a}(\ell_{t})\dot{\ell}^{a}_{t},\
G_{\ell_{0}}=1.$ (156)

Decomposing $h_{\ell_{t}}$ into

$h_{\ell_{t}}=V_{\ell_{t}}G_{\ell_{t}},$ (157)

we could obtain the equation for $V_{\ell_{t}}$,

$\partial_{t}V_{\ell_{t}}=\beta
V_{\ell_{t}}G_{\ell_{t}}K_{a}(\ell_{t})\dot{\ell}^{a}_{t}G_{\ell_{t}}^{-1},\
V_{\ell_{0}}=1.$ (158)

Let the link $\ell$ be a geodesic with respect to the intrinsic geometry given
by the triad $e^{a}_{i}=\pm\sqrt{|\det E|}^{-1}E^{a}_{i}$. A straightforward
calculation shows that

$\partial_{t}(G_{\ell_{t}}(e_{a}^{i}(\ell_{t})\dot{\gamma}^{a}_{t}\tau_{i})G_{\ell_{t}}^{-1})=0$
(159)

leading to the following property for $G_{\ell}\equiv G_{\ell_{1}}$

$G_{\ell}(e_{a}^{i}(\ell_{1})\dot{\ell}^{a}_{1}\tau_{i})G_{\ell}^{-1}=e_{a}^{i}(\ell_{0})\dot{\ell}^{a}_{0}\tau_{i}.$
(160)

Let $n_{s}$ and $n_{s}$ be given by (149) for $\vec{v}$ being
$-\dot{\ell}^{a}_{0}e_{a}^{i}(\ell_{0})/\|\dot{\ell}_{0}\|$ and
$-\dot{\ell}^{a}_{1}e_{a}^{i}(\ell_{1})/\|\dot{\ell}_{1}\|$ respectively. They
satisfy

$n_{s}\tau_{3}n_{s}^{-1}=-\frac{\dot{\ell}^{a}_{0}e_{a}^{i}(\ell_{0})\tau_{i}}{\|\dot{\ell}_{0}\|},\quad
n_{t}\tau_{3}n_{t}^{-1}=-\frac{\dot{\ell}^{a}_{1}e_{a}^{i}(\ell_{1})\tau_{i}}{\|\dot{\ell}_{1}\|}.$
(161)

Taking advantage of $n_{s}$ and $n_{t}$ and applying (160), we get that
$G_{\ell}$ takes the following form,

$G_{\ell}=n_{s}e^{\theta_{\ell}\tau_{3}}n_{t}^{-1},$ (162)

for some real number $\theta_{\ell}$.

In order to calculate $V_{\ell_{t}}$, we assume that the extrinsic curvature
$K_{ab}$ takes the form,

$K^{ab}(t)=\frac{\xi(t)}{\|\dot{\ell}_{t}\|^{2}}\dot{\ell}^{a}_{t}\dot{\ell}^{b}_{t}$
(163)

along the link $\ell$, for a field $\xi$ on $\ell$. Indeed, (163) is assumed
because the extrinsic curvature in the Regge calculus takes the same form with
a distributional field $\xi$. It is easy to verify that

$\xi(t)=K^{ab}(t)n_{a}(t)n_{b}(t)\equiv K^{nn}(t),$ (164)

where $n_{a}(t)=q_{ab}(t)\dot{\ell}^{b}_{t}/\|\dot{\ell}_{t}\|$ with
$q_{ab}=\delta_{ij}e^{i}_{a}e^{j}_{b}$ being the spatial metric. The
assumption (163) leads to

$\dot{\ell}_{t}^{a}K_{a}^{i}(t)=\dot{\ell}_{t}^{a}K_{a}{}^{b}(t)e_{b}^{i}(t)=\xi(t)\dot{\ell}^{b}_{t}e_{b}^{i}(t).$
(165)

As a consequence, we get

$G_{\ell_{t}}\left(\dot{\ell}^{a}_{t}K_{a}(t)\right)G_{\ell_{t}}^{-1}=\xi(t)\dot{\ell}^{a}_{0}e_{a}^{i}(\ell_{0})\tau_{i}.$
(166)

This equation, together with (158), results in $V_{\ell}\equiv V_{\ell_{1}}$
as follows

$V_{\ell}=\exp(\beta\left(\int_{\ell}\xi(s)\differential
s\right)\frac{\dot{\ell}^{a}_{0}e_{a}^{i}(\ell_{0})}{\|\dot{\ell}_{0}\|}\tau_{i}),$
(167)

where we have change the variable of integration to the curve length $s$.
Combining the equations (157), (162) and (167), we get

$\displaystyle h_{\ell}=\exp(\beta\left(\int_{\ell}\xi(s)\differential
s\right)\frac{\dot{\ell}^{a}_{0}e_{a}^{i}(\ell_{0})}{\|\dot{\ell}_{0}\|}\tau_{i})n_{s}e^{\theta_{\ell}\tau_{3}}n_{t}^{-1}=n_{s}\exp(\left[\theta_{\ell}-\beta\left(\int_{\ell}\xi(s)\differential
s\right)\right]\tau_{3})n_{t}^{-1}.$ (168)

This expression for $h_{\ell}$ ensures

$h_{\ell}\,e_{a}^{i}(\ell_{1})\dot{\ell}^{a}_{1}\tau_{i}\,h_{\ell}^{-1}=e_{a}^{i}(\ell_{0})\dot{\ell}^{a}_{0}\tau_{i},$
(169)

which shares a similar form as (143), but has a different sign on the right
hand side.

Now let us come to the fluxes $p^{i}_{t,\ell}$ and $p^{i}_{s,\ell}$ associated
with a surface $S_{\ell}$ intersecting $\ell$ at $x_{o}=\ell_{t_{o}}$. We
introduce $h_{\ell[0,t_{o}]}$ as the holonomy along the segment
$[\ell_{0},x_{o}]\subset\ell$, and $h_{\ell[t_{o},1]}$ as the one along
$[x_{o},\ell_{1}]\subset\ell$. By definition of the fluxes in (142), we have

$\displaystyle
h_{\ell[0,t_{o}]}^{-1}(p^{i}_{s,\ell}\tau_{i})h_{\ell[0,t_{o}]}=-h_{\ell[t_{o},1]}(p^{i}_{t,\ell}\tau_{i})h_{\ell[t_{o},1]}^{-1}=\frac{-1}{\kappa\beta}\int_{S_{\ell}}h(\rho_{x_{o}}(\sigma))E^{c}(\sigma)n_{c}(\sigma)h(\rho_{x_{o}}(\sigma)^{-1})\differential^{2}\sigma$
(170)

where $\rho_{x_{o}}(\sigma)$ is the path starting at $x_{o}$ and running in
$S_{\ell}$ till $\sigma\in S_{\ell}$.

To align with our discussion in (C.3), we introduce $p_{\ell}$ as we did in
(150). Specifically, $p_{\ell}\in\mathbb{R}$ satisfies

$(p_{\ell})^{2}=\|\vec{p}_{s,\ell}\|^{2}=\|\vec{p}_{t,\ell}\|^{2}$ (171)