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

I. INTRODUCTION

The black hole (BH), one of the most mysterious objects in our Universe, has drawn considerable attention in modern physics [1][2][3][4][5][6][7][8]. 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 [9]. The Oppenheimer-Snyder model is a solution to the Einstein equation. It describes the collapse of a spherically symmetric, pressureless and homogeneous dust ball into a BH. During the collapse, the geometry inside the dust ball is the Friedmann-Lemaître-Robertson-Walker cosmological metric, and the geometry outside is the Schwarzschild one. Consequently, the spacetime contains two singularities: the first one is the cosmological singularity emerging at the final stage of the matter collapse and the second one is the Schwarzschild singularity existing in the vacuum region outside the collapsing dust [10]. While this model plays a crucial role in transitioning BHs from abstract mathematical concepts to tangible physical realities, it faces challenges, particularly concerning the presence of the two singularities. 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 .

Among various quantum Oppenheimer-Snyder models [28,[32][33][34][35][36], recent attention has been directed toward the one that incorporates loop quantum gravity (LQG) effects [28].

This model is constructed based on the * Corresponding author; zhang.cong@mail.bnu.edu.cn fact that the metric outside can be derived without applying the Einstein equation. As shown in [37], once the interior metric describing the collapsing dust ball is given, the exterior metric can be obtained by assuming spherical symmetry and stationarity for the outside, and applying Israel junction conditions at the boundary surface between the inside and outside. Due to this fact, the quantum Oppenheimer-Snyder model in [28] considers the interior metric as the one in loop quantum cosmology instead of the classical cosmology, and derives the corresponding spherically symmetric and stationary exterior metric using the junction condition. The metric reads

where f (r) is given by

with the Barbero-Immirzi parameter ´and » = 8ÃG (see e.g., [27,36,[38][39][40][41] for similar works).

As the dust ball follows the dynamics of loop quantum cosmology, which predicts a bounce when the density of the collapsing ball reaches the Planck scale, the cosmological singularity is prevented. This will influence the nature of the Schwarzschild singularity outside the dust ball. The Penrose diagram of the Oppenheimer-Snyder model for a large BH mass M is presented in Fig. 1. As shown there, the classical singularity is replace by a transition region T connecting a white hole (WH) with the BH. (see [28] for discussion on Penrose diagrams for different value of M .)

While the spacetime illustrated in Fig. (1) offers distinct advantages, based on which observational effects have been studied [42][43][44][45][46][47], it is not without debates. For instance, the existence of the inner horizon implies that the spacetime could be unstable under perturbation [48,49]. In addition, the spacetime is not singularity free. As shown in Fig. 1, there is a timelike singularity at r = 0. 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 corner C (see Fig. 1) where the QG effects may play a significant role. As proposed in [50], the QG effects at this corner could potentially lead to quantum transition from the BH horizon to the WH horizon (see, e.g., [51][52][53][54] for previous investigation on this topic).

As proposed by [50], we cut and glue the spacetime of the quantum Oppenheimer-Snyder model as follows. First, draw a blue line as shown in Fig. (1). This blue line is symmetric under time reversal, and has null portion joining ³ and ´L and a spacelike portion joining ´L to spacelike infinity along a constant Schwarzschild time. Next, delete the entire spacetime region enclosed within the blue line, and glue its two space-like portions. This is can be done since they are both constant Schwarzschild time surfaces. The resulting spacetime, as shown in Fig. 2, exhibits a global structure with a single asymptotic region.

The metric outside the collapsing dust in revised spacetime is locally the same as given in (1.1), except for a region, called the B region, where quantum transition takes place. This fact means that the dynamics outside the B region is governed by the classical equations of motion with quantum corrections, as proposed by the LQG model for symmetric spacetimes. A question then arises regarding the dynamics within the B region: is it governed by quantum gravity? Notably, it is not viable to describe the dynamics within the B region as corrections to the classical equations of motion, since quantum tunnelling, like the one occurring in the B region, is not analytical in ℏ (see [55] for more discussion). A fundamental way to describe the dynamics of a purely quantum process is using the transition amplitude, which can be done in the framework of the covariant LQG, i.e., spin foam model [16,56,57]. Our task is to study the spin foam amplitude for for the quantum transition from the geometry in the past boundary of the B region to the geometry of its future boundary.

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 K. The intersection of K with the boundary ∂B of the B region yields a graph µ, on which one can define the boundary Hilbert space H µ of spin network states [13,15]. The spin foam amplitude A associated with the 2-complex is a specific functional on H µ , namely A : H µ → C [16]. To make the Hilbert space H µ well-defined, it is required that the boundary ∂B should be spacelike. However, the initial definition of ∂B in [50] is null. To resolve this issue, we need to choose a slightly larger region B ′ in the spacetime, ensuring that B ′ contains B and that the boundary ∂B ′ is spacelike. Indeed, all our calculations, including the triangulation aforementioned, are performed in B ′ rather than B. Throughout the subsequent discussion, we will continue to use B to denote this slightly larger region B ′ for consistency.

Another noteworthy concern is the mismatch between the amplitude A, defined on the Hilbert space H µ , and the classical geometry specifying the boundary condition on ∂B. This mismatch needs us to encode the curved continuous geometry on the discrete boundary graph µ, capturing both its intrinsic and extrinsic curvature to properly describe the physics. To this end, we choose a coherent state in H µ 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.

In our current computation, the position of the ∂B is set manually. The position of ∂B is interpreted as the moment when the QG effects cannot be neglected. Importantly, within a well-defined dynamics, this moment should be determined by the theory rather than being fixed arbitrarily. Thus, we would expect the spin foam model to determine the position of the ∂B by itself. This can be implemented as follows. We calculate the amplitude A for a family of B regions with different boundaries, and pick up the one which maximizes the amplitude. Since the value of A is related to the probability of the transition, the B region we pick up correspond to the moment when the transition is most likely to occur. The above discussion is valid for any observables which should be determined by theory itself, and explains how the spin foam model makes prediction.

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 [50] is presented in Fig. 2. The B region is the gray shaded region in the spacetime. The metric outside the region is the same as the one in (1.1). In the B region, quantum tunneling of black-to-white hole transition is supposed to occur, 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 T within B. This hypersurface T , referred to as the transition surface, divides B into a union of the past region B -and the future region B + . Let ∂ -B and ∂ + B be the remaining boundary of B - and B + except T , so that ∂B = ∂ -B ∪ ∂ + B. The time reversal symmetry gives a diffeomorphism between B + and B -that preserves the transition surface T and maps between the boundary data on ∂ + B and ∂ -B (see Fig. 2).

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. It has been noted that several previous works (e.g., [52,54]) have proposed 2-complexes for this model. However, these 2-complexes are found to be degenerate, meaning they are not dual 2-complexes of a 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 of the time reversal symmetry to transform the triangulation from B -to B + .

In the current work, we will apply a systematic way to construct the bulk triangulation for B -from triangulation for the boundary ∂ -B 1 , where we divide ∂ -B into the union ∂ -B 1 ∪ ∂ -B 2 as shown in Fig. 2. To illustrate how to construct the bulk triangulation from the boundary triangulation, in the next section, we will first introduce how the systematic way works, by using three simple examples closely related to the main problem of triangulating B -. Then, in the subsequent section II B, we will directly apply the systematic way to construct the triangulation for B -from the boundary triangulation for ∂ -B 1 .

Third example

In the above two examples, the boundary triangulation for P remains the same where P is triangulated with a single tetrahedron, while the number of dynamical points varies. In the subsequent example, we introduce a more complicated triangulation for P to demonstrate A B C D E C2 D2 A B C D D2 A B C D2 C2 E B C D D2 E B C2 D2 C FIG. 5. The triangulation of the cell ABCDED ′ C ′ (the left), where ABCDED ′ C ′ is triangulated into the four 4-simplices ABCDD ′ , ABCD ′ C ′ , BCDED ′ , and BCED ′ C ′ (the right).

how the procedure outlined in the previous examples 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. 5). Treating the points C and D as the dynamical ones, we obtain the resulting cell ABCDEC ′ D ′ , as illustrated in Fig. 5. Since the cell ABCDEC ′ D ′ contains of two subcells ABCDC ′ D ′ and BCDEC ′ D ′ , the triangulation of ABCDEC ′ D ′ is formed by combining the triangulations of these subcells. The triangulation for the subcells can be performed by following the procedure depicted in Figure 4. As a consequence, the cell ABCDEC ′ D ′ emerges by gluing together the four 4-simplices ABCDD ′ , ABCD ′ C ′ , BCDED ′ , and BCED ′ C ′ , as illustrated in Figure 5.

According the three examples, our systematic way to triangulate a 4-dimensional region can be summarized as follows: first, we perform a 3-dimensional triangulation for the past boundary which is P in the examples; second, select some points as the dynamical points and map these dynamical points to the future boundary which is P ′ in the examples, with the worlds lines connecting each dynamical point with its image point; finally, connect the image points, the non-dynamical points and the dynamical points accordingly, as illustrated by the three examples, to construct the triangulation of a 4-dimensional region.

triangulation for ∂-B1

To construct the triangulation of ∂ -B 1 , let us embed it into R 3 ∋ ⃗ x as a spherical shell, i.e., a 3D annulus, which is topologically homomorphism to ∂ -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 ∂ -B 1 proposed in [52,54], where one introduces a tetrahedron with points p a (a = 1, 2, 3, 4) to approximate the sphere S 1 . Let r 1 ⃗ x(a) ∈ R 3 for a = 1, 2, 3, 4 to be the coordinate of p a . Then, select the points p a ′ ∈ S 2 taking coordinates -r 2 ⃗ x(a). The connections between these points p a ′ (a = 1, 2, 3, 4) give rise to another tetrahedron that approximates the sphere S 2 . Finally, the triangulation of ∂ -B 1 is obtained by drawing line segments connecting each point p a (for a = 1, 2, 3, 4) in S 1 to the corresponding points p b ′ , p c ′ , and p d ′ in S 2 , where b, c, and d represent labels that are different from a. Importantly, the minus sign in the coordinate of p a ′ ensures that this connectivity structure yields a well-defined triangulation. The resulting triangulation can be visualized as shown in Fig. 6. 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 p a p b ′ p c ′ p d ′ .

• B-type: There are 6 tetrahedra having two points in S 1 and two points in S 2 . The points of the tetrahedron are p a p b p c ′ p d ′ . • C-type: There are 4 tetrahedra having three points on S 1 and one point on S 2 . The points are p a p b p c p d ′ . Here, a, b, c, d represent distinct elements of the set {1, 2, 3, 4}. Once the triangulation is done in R 3 , one can pull back it to ∂ -B 1 to get the triangulation for ∂ -B 1 .

In what follows, the pulling-back of the points p a and p a ′ will still be called p a and p a ′ respectively.

triangulation for B-

To obtain a triangulation for the 4-dimensional region B -from the triangulation of its boundary ∂ -B 1 , we selected the points in S 2 to be dynamical. Starting from each of these dynamical points p a ′ for a = 1, 2, 3, 4, a line p a ′ p a ′′ is drawn in ∂ -B 2 , intersecting the sphere S 3 at a chosen point p a ′′ . p a ′ p a ′′ is viewed as the worldline of p a ′ in our final triangulation for B -. Following the method outlined in Sec. II A, 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 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 ∂ -B 2 . Then, p a ′′ for a = 1, 2, 3, 4 are chosen so that their coordinates are r 3 ⃗ x(a) ∈ R 3 . The resulting triangulation for the transition surface T , as the boundary of the triangulation of B -, is the same as that for ∂ -B 1 , i.e., the same as shown in Fig. 6.

The 2-complex K dual to the triangulation of B -is shown in Fig. 7. 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 K,

where the domain of 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 K are provided in [58].

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.

C. Boundary state

The intersection between K and ∂B gives the boundary graph µ. µ contains the nodes n e and links ℓ f , which result from the intersections of edges e and faces f , respectively, with ∂B. All nodes of µ are 4-valence. The boundary states are vectors in the LQG Hilbert space H µ on the graph µ. We will call H µ the boundary Hilbert space in the current work. Specifically, H µ is given by

where H j is the SU(2) representation space with spin j and H * j is its dual space. The spaces H j ℓ and H * j ℓ are associated with the starting node s ℓ and targeting node t ℓ of ℓ respectively.

As discussed in Sec. I, the current work uses the coherent states as the boundary states of spinfoam amplitude. In LQG, there are several types of coherent state such as the coherent Livine-Speziale intertwiners [59] and the Thiemann's complexifier coherent state [60][61][62]. They all have good properties and well studied in literate. However, the coherent Livine-Speziale intertwiners do not contain information about the embedding of the boundary. We thus adopt the Thiemann's complexifier coherent state (see, e.g. [63,64] for the earlier applications to the spin foam model). Given a link ℓ ≡ ℓ f , a Thiemann's complexifier coherent state labelled by g ℓ ∈ SL(2, C) can be expressed as a function of the holonomy h ℓ ∈ SU(2) along ℓ as

with N denoting the normalization factor and t being the classicality parameter which is dimensionless and equal to ´»ℏ divided by some area. Here the Wigner-D matrix D j nm (g ℓ ) for g ℓ ∈ SL(2, C) is defined as

where we decompose g ℓ ∈ SL(2, C) as

with n ℓs ℓ , n ℓt ℓ ∈ SU(2) belonging to the starting and target points of ℓ and p ℓ , ϕ ℓ ∈ R. It can be checked that the state (2.3) is peaked on the elementary operators (refer to Appendix C 3 and [65][66][67]):

where pk s,ℓ and pk t,ℓ are the flux operators and ĥℓ denotes the holonomy operator (refer to Appendix C 2). In addition, the spread of the state goes like t, because the uncertainties of the elementary operators are

In Appendix B, the explicit expression of the spin foam amplitude with the boundary state (2.3) is derived. Since the current work concerns with the asymptotics as t → 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

) where À x for x = (ℓs ℓ ) or (ℓt ℓ ) is the spinor given by

with n x for x = (ℓs ℓ ) or (ℓt ℓ ) given in (2.4), and |j ℓ , À ℓt ℓ ð ∈ H j ℓ and ïj ℓ , À ℓs ℓ | ∈ H * j f are the SU(2) Perelomov coherent states [68][69][70]. Geometrically, the SU(2) coherent state |jÀð ∈ H j corresponds to the vector jÀ ⃗ ÃÀ ∈ R 3 . The state (2.6) is the similar coherent state used in the usual context of the spin foam model [63], and its generalization to higher dimensional LQG is also studied [71,72]. In terms of the state (2.6), (2.5) becomes

where n(À) for À = (À 1 , À 2 ) T ∈ C 2 is given by

(2.9)

In terms of these parameters, g ℓ is expressed as

(2.10)

We get the coherent state on the entire graph µ as their tensor product,

(2.11)

Even though the boundary state given in (2.11) 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 ℓ and À ℓ in the coherent state carry the information about the densitied triad E a i on ∂B, while ϕ ℓ is relevant to the Ashtekar connection A i a on ∂B (see also e.g. [73]). Therefore, the coherent state carries the information about the 4-dimensional geometry of ∂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 = ¶ ij e i a e j b . This relationship suggests an O(3) gauge transformation, wherein the triad undergoes a transformation e i a → O i j e j a with O ∈ O(3). However, in LQG, only 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 , for each link ℓ, to denote the label of the coherent state constructed from the triad (sẽ i a , sē i a ) with s, s = ±. To keep the discussion flow, we will momentarily skip detailing the explicit relation between g (ss) ℓ and (sẽ i a , sē i a ). The above discussion motivates us to consider the boundary state as the superposition

(2.12)

Then the spin foam transition amplitude is a sum

(2.13)

D. Spin foam amplitude and action

The spin foam amplitude is a specific linear functional acting in the boundary Hilbert space (2.2) [16]. The spinfoam amplitude with the boundary state (2.11) has an integral expression. This expression assigns an SU(2) representation j f to each face f , an SL(2, 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 ∈ C 2 to each corner (v, f ) of the face f at the vertex v. The amplitude reads

where |V f | is the number of vertices surrounding the face f , dg ve is the Haar measure, d 2 Ω vf is a measure on CP 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 dg ve and d 2 Ω vf can be found in [74]. The integration for the spinor z vf ∈ C 2 is over CP 1 , relying on the scaling invariance of the action S. The action S can be written as the following summation

where S f is the action associated with the face f . Specifically, for an internal face f , S f reads

(2.16) where Z vef := g ve z vf , s e and t e represent the source and target points of edge e respectively, ∥Z∥ denotes ïZ, Zð for all Z ∈ C 2 , and ´is the Barbero-Immirizi parameter as introduced in (1.2). For a boundary face f , as shown in (B25) (see Appendixes A and B for more details), we have

(2.17) where the extra term compared to (2.16) comes from the coherent state. Moreover, in (2.17), Z seef = À ℓ f s ℓ f when the boundary of the face f starts from e, and Z teef = À ℓ f t ℓ f when the boundary of the face f ends at e, with À ℓ f s ℓ and À ℓ f t ℓ introduced in (2.6).

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

Considering arbitrary SL(2, C)-valued field v → h v at the vertices, arbitrary complex field (v, f ) → ¼ vf with ¼ vf ̸ = 0, and arbitrary SU(2)-valued field e → u e on the edges, we observe that S exhibits invariance under the transformation:

The integration in (2.14) is carried out over the quotient space Γ/ ∼, where the equivalence relation ∼ is defined by the gauge transformation (2.18).

III. BOUNDARY DATA

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

Before going to the specific calculation, let us outline the logic of this section. The final goal of the current work is to numerically solve the equation of motion ¶S = 0 with the action given in (II D). To this end, a starting point for numerical root finding is needed. This starting point can be constructed based on a flat geometry which has the intrinsic and extrinsic boundary geometries approximating those of ∂B. Since the boundary ∂B, as a hypersurface in the spacetime depicted in Fig. 1, is not fixed so far, we can thus choose it such that ∂B takes the same intrinsic geometry as the hypersurface surface t = ar, with an arbitrary real number a satisfying |a| < 1, in the Minkowski spacetime R 4 . This choice will help us achieve the flat geometry needed for constructing the starting point for root finding. Here we use the spherical coordinate (t, r, ¹, ϕ) in the Minkowski spacetime such that the Minkowski metric takes ds 2 = -dt 2 + dr 2 + r 2 d¹ 2 + r 2 sin 2 ¹dϕ 2 .

With ṽ := t + r, the Minkowski metric becomes

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

Consider the Eddtington-Finkelsein coordinate where the BH metric reads

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

for some function F . With a straightforward calculation, we get

where 1 -(1a 2 )f (r) > 0 for all r g 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 . It is reminded that Region B is outside the dustball.

If the positive sign is selected, the right hand side of (3.5) is singular at r = r ± where f (r ± ) = 0. Thus, we should choose

with some a 1 ; and

(2)

for some a 2 .

With the past boundaries ∂ -B 1 and ∂ -B 2 , we get the corresponding future boundaries ∂ + B 1 and ∂ + B 2 by the time reversal symmetry.

It is easy to check the outwards normal

Their lengths for n i is

Thus, the boundary

(3.10)

Similarly, we get

where ¹ 2 is the dihedral angle between ∂ -B 2 and ∂ + B 2 .

Eqs. (3.10) and (3.11) establish a connection between the parameters a i and geometrical properties of the region B.

In our practical numerical computation, we will select ¹ 1 and ¹ 2 as the parameters to set, rather than a 1 and a 2 .

A. Values of p ℓ and ξ ℓn

In (2.8), the quantities p ℓ and À ℓn are related to the expectation value of pk s,ℓ , which denotes the area vector of the surface for the flux. In our work, this surface corresponds to a face within the triangulation of ∂B, which serves as a discrete approximation of the geometry of ∂B. Therefore, determining p ℓ and À ℓn involves computing the area vectors of surfaces within the discrete geometry, constructed as outlined below.

Our choice of the boundary ∂B ensures that ∂ -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 i of ∂ -B into the Minkowski spacetime, such that the 3-metric of i(∂ -B) reduced from the Minkowski metric is the same as the 3-metric of ∂ -B reduced from the BH metric. Since ∂ -B and i(∂ -B) have the same intrinsic geometry, we can calculate the p ℓ and À ℓn , encoding the information of the intrinsic geometry of ∂ -B, in the Minkowski spacetime.

To calculate p ℓ and À ℓn , we first embed the triangulation for B -constructed in Sec. II B into i(B -), a region in the Minkowski spacetime. As discussed at the beginning of Sec. II B, we initially focused on the topological structure of the triangulation. Now, the geometric details, such as the shape of the tetrahedra in the embedded triangulation, must be focused on, to calculate p ℓ and À ℓb . It is observed that the triangulation constructed in Sec. II B is completely determined by the points p a , p a ′ and p a ′′ for all a = 1, 2, 3, 4. Thus, we need to specify the position of points i(p a ), i(p a ′ ) and i(p a ′′ ) in the embedded triangulation corresponding to p a , p a ′ and p a ′′ . Specifically, we require that the tetrahedra i(p

), are equilateral. In addition, Let r 1 , r 2 and r 3 be the radii of the spheres i(S 1 ), i(S 2 ) and i(S 3 ), and t 1 , t 2 and t 3 be the time coordinate of the spheres. Under the Cartesian coordinate system, we re-

The above requirements maximally respect the spherically symmetry. The embedded boundary triangulation provided a discrete geometry approximating the intrinsic geometry of i(∂ -B). Using the discrete geometry, we can calculate p ℓ and À ℓf by

where ⃗ A ℓ is the area vector of the triangle dual to the link ℓ. Since p ℓ and À ℓn for all links ℓ connecting to each node n are determined by a tetrahedron, the closure condition

must be satisfied, where the sign factors ϵ ℓn depends on the direction chosen for ⃗ A ℓ . Let r : ∂ -B → ∂ + B denote the time reversal symmetry. Owing to this symmetry, ∂ + B exhibits the same intrinsic geometry as that of ∂ -B. Thus, if ēi a is the cotriad field in ∂ -B, the cotriad field in ∂ + B could be r * (ē i a ), i.e., the push-forward of ēi a . Here the push-forward of the forms ēi a is regarded as their pull-back by r -1 . Despite the time reversal symmetry, the 3-orientation of ∂ + B, induced by the outward normal field, is opposite to r * (o -), namely, the push forward of the reduced orientation o - on ∂ -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 ≡ det e j b . It is dual to a 2-form E a i ε abc , whose surface integration, defining the flux, is 3-orientation dependent (refer to (C7)). Consequently, the flux at each node r(ℓ)r(n) ¢ ∂ + B becomes the opposite of the flux at the corresponding (ℓ, n). In other words,

with pk n,ℓ = pk s,ℓ for n = s ℓ and pk n,ℓ = pk t,ℓ for n = t ℓ . Combining (3.14) with (2.8), we get

where p ℓ and À ℓn are the data associated with ℓ, n ¢ ∂ -B, and p r(ℓ) and À r(ℓ)r(n) are the corresponding data for r(ℓ)

and r(n). The equation (3.15) implies

where x is an arbitrary phase, and

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

In our computation, we choose the surfaces i(

It is reminded that i, as defined in the paragraph before (3.12), denotes the embedding of

we set the Cartesian coordinates of i(p a ), i(p a ′ ) and i(p a ′′ ) in Minkowski spacetime to be

B. Relate φ ℓ to extrinsic curvature

According to (2.5), we need to calculate the holonomy within the boundary geometry to get ϕ ℓ . Appendix C provides a detailed computation, with key aspects worth noting. At first, as shown in [75], the value of the spin foam action at the real critical points differs from the Regge action by a factor sgn(V 4 (v)). This extra factor suggests that the volume form in classical action corresponding to the spin foam model is det(e)d 4 x where 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 ℓ in the holonomy is chosen to be a geodesic for simplicity. Finally, we assume that the extrinsic curvature K ab along the link ℓ takes the form

with À being a field on ℓ and Ä being the parameter of ℓ. This assumption makes the computation tractable. In LQG, quantities are often smeared by integration. Therefore, the simplified expression (3.23) 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 ϕ ℓ as follows:

where ¹ ℓ is the twist angle, Ξ ℓ is related to the extrinsic curvature and sgn(e) ≡ sgn(det e i a ) is computed with any right-hand coordinate system.

The twist angle ¹ ℓ will be fixed by choosing some certain frame for the two tetrahedra sharing the triangle dual to ℓ. 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 Ξ ℓ which, due to the assumption (3.23), reads

where s is the parameter of length of ℓ, la s is the tangent vector ℓ with respect to the parameter s.

Turning our attention to the continuous geometry on ∂ -B, a naive approach to relate Ξ ℓ to this geometry involves a direct application of (C47) where one would integrate the extrinsic curvature K ab of the continuous geometry along a single edge. However, this approach cannot fully encode the information of K ab within Ξ ℓ . 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 ℓ , where it remains constant. Thus, integrating K ab along a single edge, as in (C47), 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 ∂ -B. Consider two tetrahedra T 1 = q 1 q 2 q 3 q 4 and T 2 = q 2 q 3 q 4 p 5 , sharing the interfacing triangle t ℓ = q 2 q 3 q 4 . Here q i for i = 1, 2, 3, 4, 5 represents a point in the triangulation of ∂ -B. Introducing the "centroids" c 1 ∈ T 1 and c 2 ∈ T 2 , we define a convex region c 1 q 2 q 3 q 4 c 2 . This region is foliated into a family of two surfaces Ä → A(Ä ) ¢ R 2 with A(0) = t ℓ . The average extrinsic curvature K(Ä ) is defined on each surface as

where ∥A Ä ∥ is the area of A Ä and n a is the unit normal of A Ä . K(Ä ) in (3.26) captures the information K ab n a n b on the surface A Ä . Choosing a link ℓ connecting q 1 and q 2 and transversely intersecting each A Ä , we relate Ξ ℓ to the extrinsic curvature by

This equation defines a connection between the extrinsic curvature and the data Ξ ℓ inherent in the coherent state. this definition reverts to (3.25) in the specific case of Regge geometry (refer to the discussion in Appendix C 4). Importantly, as (3.27) defines Ξ ℓ through the average extrinsic curvature, it effectively encapsulates information about the extrinsic curvature across the entire region c 1 q 2 q 3 q 4 q 2 . Our actual computation adopts the coordinate give by (x, y, z) = (r sin ¹ cos ϕ, r sin ¹ sin ϕ, r cos ¹) for simplicity. Within this coordinate, we choose c 1 := 1 4 (q 1 + q 2 + q 3 + q 4 ) and c 2 := 1 4 (q 2 +q 3 +q 4 +q 5 ). The region c 1 q 2 q 3 q 4 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 c 1 q 2 q 3 q 4 cc 2 , the slices A Ä are chosen as plane surfaces parallel to t ℓ . The link ℓ, connecting c 1 and c 2 , is selected as the straight line t → (1t)c 1 + tc 2 in this coordinate system.

Due to the time reversal symmetry r :

(3.28)

C. Superpose the coherent states

Let us firstly summarize the results given in Secs. III A and III B. Given a link ℓ in ∂ -B, we can calculate p ℓ , À ℓs ℓ , À ℓt ℓ , ¹ ℓ and Ξ ℓ , to get the coherent state label ḡℓ as

where we choose p ℓ > 0 without loss of generality and ēi a denotes the triad field in ∂ -B. Then, by choosing p r(ℓ) = -p ℓ that satisfies (3.18), we get the coherent state label gr(ℓ)

where we combine (3.18), (3.24), (3.28) and (2.10), and consider the discussion in the paragraph before (3.14) showing that sgn(r * ē) = -sgn(ē).

To maintain a consistent global orientation for the two complexes dual to the triangulation of B, the boundary face in ∂ + B corresponding to a face f ¢ ∂ -B is represented by r(f ) -1 . In simpler terms, this means that the corresponding face includes the same points as r(f ) but with the opposite orientation. As a result, for a link ℓ ¢ ∂ -B, its counterpart in ∂ + B is given by r(ℓ) -1 . By (3.30), the coherent state associated with r(ℓ) -1 is labelled by

As discussed in (2.12), we are considering superposition of various coherent states resembling the same geometry. Following the convention introduced in (2.12), the values for g for other signs of s 1 and s 2 . Consider the transformation ēi a → -ē i a in ∂ -B. This transformation does not alter the densitized triad field E a i in ∂ -B, and consequently, it does not change À ℓn and p ℓ . Furthermore, the spin connection Γ i a and the extrinsic curvature K ab remain unaffected. As a result, ¹ ℓ and Ξ ℓ remains constant under this transformation. However, due to the relationship sgn(det ēi a ) = -sgn(det -ē i a ), ϕ ℓ will be transformed to

This leads to the following transformation for the coherent state label ḡℓ ,

(3.33) Similarly, under the transformation r * (ē i a ) → -r * (ē i a ), gr(ℓ) -1 for all r(ℓ) transformed as

.34) As a summary, the list of coherent state labels {g (ss) ℓ } for links lying entirely in either ∂ + B or ∂ -B is given by {g

where we defined g(+) r(ℓ) -1 = gr(ℓ) -1 and ḡ(+) ℓ = ḡℓ .

D. 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 (3.35). 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, ℓ = l ∪ l, where l = r( l) -1 . In considering the values of g (ss) ℓ for different signs of s, s, we correspond to intrinsic geometries given by (sẽ i a , sē i a ) = (sr * (ē i a ), sē i a ) (refer to (2.12) for the convention). For each of the segments l and l, there exist associated coherent states È t g(s) l =n(À lsl )e (-ipl+ϕ s l ) n(À ltl ), ḡ(s) l =n(Àl sl )e (-ipl+ϕ s l ) n(Àl tl ),

where it is considered that the signs s and s only affect the value of the angle ϕ, as indicated by (3.29), (3.30), (3.33), and (3.34). Due to he time reversal symmetry, we have

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

(3.38) where h ℓ = hlh l is the holonomy for the entire link. Noticing the change of parameter t → 2t, we find that the state (3.38) is approximated by |g (ss) ℓ ð (refer to (2.6)) with

(3.39) Indeed, rather than the state (3.38), we prefer to take the boundary state associated with ℓ as È t g (ss) ℓ . Even though both of the two state can be approximated by |g 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 [76][77][78][79]. This method involves obtaining the complex critical point for each set {g (ss) ℓ } ℓ¢∂B with s, s = ±1. Since the procedure for each {g (ss) ℓ } ℓ is the same, we will use g ℓ to represent either of g (ss) ℓ in the following. By definition, the parameter t is equal to ´»ℏ divided by some unit of area. In our computation, the unit, denoted by ∆, will be chosen as some typical area in the triangulation. Then, in the computation of the spin foam action (2.17), the parameter t becomes

Then, the effective action for a boundary face can be expressed in terms of ∥ ⃗ A∥ as

where (3.12) is applied, and all of the remaining terms are absorbed into

This allows the action for a boundary face to take the form t -1 S f . It is noteworthy that the substitution j f → 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

where the Poisson resummation formula is employed to approximate the summation by an integral (see [77] for more details). Thank to (4.4), we are allowed to discuss the asymptotic behavior of A as t -1 → ∞. This involves the complex critical point method.

A. Complex critical point method

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

over some N -dimensional configuration space C with some parameter t. The behavior of A t as t → 0 is guided by the critical points of S, where a critical point is a point

It is worth noting that, in general, the equations ∂ 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 C C of C and satisfies ∂ z S = 0, where S needs to be analytically continued to C C ∋ z.

In our computation, we will introduce some external parameters Ç ∈ [Ç o , Ç 1 ] as in [77], to rewrite the integral as

Here S(x, Ç) is assumed to be an analytical function with respect to Ç. Additionally, Ç should be selected in a way that S(x, Ç 1 ) precisely describes our model of interest. Simultaneously, for Ç = Ç o , S(x, Ç o ) may represent a distinct model, but it possesses a real critical point 1 . Let Z(Ç) be the solution to

Here, the real critical point is defined as the solution to both ∂ x S(x, Ç o ) = 0 and Re(S(x, Ç o )) = 0. As Ç deviates away from Ç o , Z(Ç) moves away from C to C C analytically. Consequently, Z(Ç 1 ) ∈ C C will be a complex critical point contribute the integral in our initial model (see Fig. 8 for the illustration). With the complex critical point, the small-t asymptotic expansion for the integral in A t can be established as

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

where C is a positive constant [80]. This condition implies that in (4.5), the purely imaginary phase can only occur at the real critical point, where Im(Z) = 0 and Ç = Ç o . When Ç deviates from Ç o , causing Im(Z) finite and Re(S) to become negative, the result in (4.5) is exponentially suppressed as t approaches 0. Now let us return to the integral (4.4) to apply this method. As discussed in Sec. II D, the spin foam action is a function on the space Γ ∋ (j f , g ve , z vf ). However, 1 As we will see, in the concrete example, χ will be related to boundary data.

Re(Z )

FIG. 8. The real and complex critical points xo and Z(χ). S(z, χ) is analytic extended from the real axis to the complex neighborhood illustrated by the red disk.

the integral in (4.4) is merely performed over the quotient space Γ/ ∼, where the gauge orbit given by (2.18) 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

where ∂ ′ denotes the derivative only along the direction transversal the gauge orbit.

As illustrated in Fig. 8, the process of obtaining a solution is divided into two steps. The first step is to find a real critical point, x o . In the second step, x o is perturbed along the flow parameterized by Ç until reaching Ç 1 , resulting in a solution to (4.8) and (4.9). In the spin foam model, a real critical point can be constructed from flat geometry. To construct such a critical point, it is necessary to modify the boundary data of our model of interesting to ensure that the modified boundary data permit a flat bulk geometry. The critical point x o can then be constructed from this flat geometry. To better understand the role played by x o , the following observation is helpful. Consider selecting the boundary coherent state corresponding to the boundary data of flat geometry to obtain the spin foam effective action, denoted by S. In this case, x o becomes the solution to the critical equations (4.8) and (4.9), with S replaced by S. It is important to note that the critical equation (4.8) associated with S remains the same as the original one associated with S, while (4.9) differs between S and S. Therefore, x o is a solution to (4.8) but not to (4.9). To obtain a solution to both (4.8) and (4.9), we need to perturb x o . This resulting solution is complex. This observation also helps guide the choice of the parameter Ç. By definition, Ç is chosen to connect the concerning model with a model possessing a real critical point. Given the above understanding, we need to establish a relationship be-tween the parameter Ç and the boundary data, so that S(x, Ç o ) ≡ S.

To maintain the flow of the discussion, the subsequent subsections will disregard the parameter Ç and focus on obtaining x o from geometry, as well as searching for complex critical point from x o . When the parameter Ç is considered, the procedure remains similar. We still need to obtain x o from geometry initially, but search for the critical point Z(Ç) associated with S(x, Ç) from x 0 for all Ç ∈ [Ç o , Ç 1 ] 2 . The exact definition of Ç will be provided in Sec. V B. From now on, x o will be referred to as the real critical point, while it does not fully satisfy the critical equations (4.8) and (4.9).

B. 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 ∆ as the unit. Under this unit, the parameter t is equivalent to the Planck area ´»ℏ.

As discussed below (4.8) and (4.9), the real critical point is a solution only to (4.8). The absence of (4.9) allows us to disregard the boundary conditions arising from the coherent state data p ℓ , À ℓ , and ϕ ℓ . In other words, with any boundary data that permits a flat geometry, one can construct a corresponding real critical point x o , which will serve as a solution to (4.8) but not to (4.9). However, since the goal of introducing x o is to solve both (4.8) and (4.9), we will select the boundary data corresponding to x o as closely as possible to that of our model, ensuring that x o is an approximate solution to (4.9). Specifically, the boundary data for x o will have the same intrinsic geometry as our model, while differing in extrinsic geometry. This means that the area vectors for each link ℓ will remain as p ℓ and À ℓ , but the dihedral angles associated with each link will differ from ϕ ℓ . By adjusting the dihedral angles, we can achieve a flat geometry.

Given a set of boundary data {(p ℓ , À ℓ )} ∀ℓ∈∂B , the typical procedure to find the real 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 ℓ , À ℓ )} ∀ℓ∈∂B is flat. In the case at hand, the answer is negative. As we will see, it is possible to construct a flat geometry in the region B -that satisfies the boundary data on ∂ -B. 2 In the practical computation, we search for the critical point Z(χo + δχ) associated with S(x, χo + δχ) with δχ j 1 from x 0 , then search for the critical point Z(χo + 2δχ) associated with S(x, χo + 2δχ) from Z(χo + δχ), and so forth, until we obtain the critical point associated with S(x, χ 1 ), i.e., the solution to both (4.8) and (4.9).

However, extending this flat geometry into the region B + results in a mismatch with the boundary conditions on ∂ + B. Now, a tension arises: we are faced with two options. We can either construct a flat geometry throughout the entire region B, which will not match the intrinsic boundary geometry, or we can construct a geometry that matches the intrinsic boundary geometry but is not flat throughout the entire region. In this work, we choose the latter option. Specifically, we will construct a flat geometry in the B -region that matches the intrinsic boundary geometry of ∂ -B, and then reflect this geometry into the B + region, due to the time reversal symmetry, to obtain a geometry for the entire B region. However, this resulting geometry is only partially flat, as there will be a deficit angle present on the transition surface T . In Sec. III A, to construct the discrete geometry approximating the intrinsic geometry of ∂ -B, we introduced the points i(p a ), i(p a ′ ) and i(p a ′′ ) 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 p a , p a ′ and p a ′′ . We thus replicate the connections between the points p a , p a ′ and p a ′′ , to connect the points i(p a ), i(p a ′ ) and i(p a ′′ ). This process results in a flat Regge geometry in the region B -. Utilizing the geometric interpretation of the real spin-foam critical points [75,81,82], 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 ℓ and À ℓ , these data are clearly compatible with the boundary data on ∂ -B. Nevertheless, these data are specifically defined on the 2-complex in B -and solely solve the equations in (4.8) associated with j fi , g ve , and z vf within B -. They even do not solve the equations ∂ j f i S = 0 for the faces f i across the transition surface T .

Given a vertex v, an edge e and a face f contained in B -, the time reversal symmetry r : B -→ B + defines the corresponding objects r(v), r(e), 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 r(f ) , g r(v)r(e) , z r(v)r(f ) ) = (j f , g ve , z vf ) .

(4.10)

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 ∂B and solve the equations in (4.8) 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.

C. Complex critical point

The absence of (4.9) 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 Ξ ℓ computed in Sec. III B. Introducing the dihedral angles Ξ ℓ 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 [83], we cannot obtain a real solution to (4.8) and (4.9), 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 (4.8) 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 ( jf ,g ve , zvf ). 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 veo(v) as

Here, µ vo(e)e , along with analogous parameters introduced subsequently, is determined by det g vo(e)e = 1. This fixes the SU(2) gauge transformation in (2.18). Note that (4.12) needs to be applied to all internal edges, including those e o (v) that were chosen for gauge fixing the SL(2, C) transformation. For each e o (v) ≡ ẽ, its corresponding v o (ẽ) 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

Furthermore, the spinors z vf are parametrized as

which fixes the rescaling transformation of z vf in (2.18). Finally j f is parametrized as

Substituting the above parametrization into the action (2.15), 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 ∂ z S = 0 and ∂ z S = 0 for z taking the complexification of x i ve , y i ve , x vf or x f , where z denote the complex conjugate of z.

A technique issue appearing here is the incorporation of the value Ξ ℓ provided in Sec. III B into our actual computation. This issue is noteworthy due to the lack of explanation regarding the definition of the twist angle ¹ ℓ . According to the expression of the action given in (II D), S takes the form

where ϕ ℓ contains the twist angle ¹ ℓ and the dihedral angle Ξ ℓ . We introduce φℓ as the coherent label given by

.17) φℓ ensures ∂ j f b S takes vanishing value at x o for all boundary face f b . As we discussed in (C 6), φℓ f is equal to the twist angle plus a term proportional to the dihedral angle. This expression can be viewed as the discrete counterpart of the Ashtekar connection, A i a = K i a + µΓ i a

, where K i a corresponds to the dihedral angle and Γ i a to the twist angle [73,84]. Thus, the twist angles depend only on the intrinsic geometry which remains consistent at both the real and the complex critical points (see [84] for more details). As a consequence, the twist angle in φℓ is identical to that in ϕ ℓ . Combining this observation with (3.24), we arrive at the expression

where Ξℓ = -∂ ´φ ℓ denotes the dihedral angle compatible with the flat Regge geometry associated with x o . Substituting the result into (4.17), we get

For a given real critical point x o = ( jf ,g ve , zvf ), the parity transformation

yields another real critical point denoted by xo . As shown in [75], under the parity transformation, Ξℓ is transformed as

Ξℓ → -Ξℓ . (4.21)

In the numerical procedure for searching the complex critical point, the initial point should be chosen from either x o or xo such that the corresponding Ξℓ is closer to Ξ ℓ .

V. NUMERICAL INVESTIGATION

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

In this case, the parameter ³ takes the value

´»ℏ.

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.

A. Fix the location of ∂B Even though the free parameters as shown above are fixed, we still have the ambiguity on choosing the size of the B region in the entire spacetime. This needs us to fix the location of ∂B. As outlined in the beginning of Sec. III, this requires us to set the parameters r a with a = 1, 2, 3 and ¹ i with i = 1, 2 for ∂ -B. Subsequently, the boundary ∂ + B can be determined accordingly, taking into account the time reversal symmetry.

In our computation, we choose

(5.1)

Here, r 1 is chosen as the minimal radial where the collapsing dust ball bounces in the quantum Oppenheimer-Snyder model. r 2 < r 3 are choose as two radials nearby the horizon. Consequently, the resulting B region spans from the minimal radius up to just outside the horizon, effectively capturing the region expected to contain the quantum gravity regime. Additionally, r 2 and r 3 are positioned close to the horizon, ensuring that the B region remains compact and computationally manageable.

We select the values of ¹ 1 and ¹ 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 ¹ 1 and ¹ 2 . In the current work, we choose the following values as an example

When we embed the region B -into Minkowski spacetime using the map i, ¹ (c) 1

and

are the critical values to ensure the following conditions (see App. E for detailed derivations):

(1) All tetrahedra in the transition surface i(T ) are spacelike;

(2) The boundaries i(∂ -B 1 ) and i(∂ -B 2 ) are spacelike;

(3) The time coordinates t 2 and t 3 of the spheres i(S 2 ) and i(S 3 ) respectively satisfy t 2 < t 3 , so that i(S 3 ) occurs in the future of i(S 2 ).

Once the boundary ∂B is fixed, the unit ∆ is chosen as the minimal area in the triangulation of ∂B. In our case, it is the area of the triangle p 1 p 2 p 3 . We have

)

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

C. Property of the real critical point

As mentioned earlier, we need to find the complex critical point for each set {g (ss) ℓ } ℓ¢∂B with s, , s = ±1. However, the following considerations suggest that we only need to focus on the case where s = s.

At first, according to the discussion for (4.6) and (4.7), 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 (5.8), where the coefficient of the quadratic term is proportional to the sum of (ϕ ℓ f b -φℓ f b ) 2 for all boundary faces f b . Consequently, instead of considering all possible {g (ss) ℓ } ℓ¢∂B with s, s = ±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 φℓ defined in (4.17) satisfies φℓ = -φr(ℓ) -1

(5.12)

and, for each link ℓ c at the corners, φℓc = 0.

(5.13)

As a consequence, among the boundary conditions corresponding to {g (ss) ℓ } ℓ¢∂B with s, , s = ±1, those closer to the boundary conditions given by x o are the ones satisfying

(5. Indeed, given a link ℓ c at the corner, for s = -s, (3.40) yields |ϕ ℓc | = 2|Ξ ℓc |, relating |ϕ ℓ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 ∂B ensures its local intrinsic flatness, the extrinsic curvature nearby S 1 is expected to be very large. That implies a large value of |ϕ ℓc | for ℓ c passing through S 1 . Namely, ϕ ℓc for ℓ c passing through will deviate significantly from φℓc = 0. Consequently, this leads to a highly negative value of the real part of the action at the complex critical point for s = -s, as indicated in (4.7).

D. Final results

As discussed in the last subsection, our attention should be solely directed towards the complex critical point for the boundary data {g (++) ℓ } ℓ¢∂B and {g (--) ℓ } ℓ¢∂B . The detailed numerical results are available in [58]. We denote the values of Sres associated with these boundary data as S(+) res and S(-) res respectively. Our numerical investigation reveals that both S(+) res and S(-)

res are real values, as illustrated in Fig. 9. According to the numerical results, for Ç = 1, we have S(+) res Ç=1 = -0.0458193513442056, S(-)

res Ç=1 = -0.0458193513442275.

(5.15)

It is interesting that there exist a tiny difference between S(+) res and S(-) res . According to the numerical results in Fig. 9 with the residuals satisfying |ϵ ± (Ç)| < 2 × 10 -12 .

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. Namely, the orientation of 4-simplices in B -could differ from one 4-simplex to another.

By computing the deficit angle ϑ pap b ′ p c ′′ hinged at the internal triangle p a p b ′ p c ′′ ¢ B -, we obtain

(5.17)

Let Θ eve ′ be the dihedral angle between the 2 tetrahedra e and e ′ in the 4-simplex v. Then, we have

where the sum is taken over all 4-simplices sharing p a p b ′ p b ′′ . Equation (5.17

Furthermore, the fact that x o is a real critical point implies

where V v represents the volume of the 4-simplex v. Consequently, sgn(V v ) varies among the 4-simplexes sharing the triangles p a p b ′ p b ′′ . For example, we consider the case where a = 2 and b = 1. The triangle p 2 p 1 ′ p 1 ′′ is shared by four 4-simplices

We find

Indeed, for each triangle p a p b ′ p b ′′ , its boundary always contains a unique 4-simplex v 1 given by v 1 = p a p c p d p b ′ p b ′′ , where c and d are neither a nor b. The numerical computation confirms that the sign of V v1 is opposite to the signs of all other 4-simplices sharing the triangle.

The dihedral angles occurring on the transition surface T are associated with the triangles in the triangulation of 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 p 1 p 2 p 2 ′′ . 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 p 1 p 3 ′′ p 4 ′′ . In the current numerical computation, the dihedral angle of the first type takes the value

and the second type takes

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 [50]. To construct the spin foam amplitude, we first introduce a systematic procedure to yield a bulk triangulation from the boundary triangulation in the Bregion. 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 the 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 ¶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 C, the amplitude is dominated by the terms allowing the change of orientation from the past boundary ∂ -B to the future boundary ∂ + B. This result indicates that the black-to-white hole transition should be accompanied by the quantum tunneling process of the change of orientation from ∂ -B to ∂ + B. It is noted that the change of orientation was also discovered in the loop quantum reduced model [40]. 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 (++) ℓ } ℓ¢∂B and {g (--) ℓ } ℓ¢∂B as described in Sec. III C. According to the numerical results presented in Sec. V D, 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.

In this work, we apply the conjecture introduced in previous research (see, e.g., [52]). In these studies, the region of classical spacetime where QG effects become relevant is divided into three sub-regions, with the B region-our focus in this work-being one of them. The QG effects in these three sub-regions are addressed using different theories. Specifically, different from the B region, the other two regions are handled using effective loopquantum-cosmology-like dynamics. Although the results from loop-quantum-cosmology-like dynamics have gained broad acceptance, a more consistent approach would be to study QG effects across all regions within a unified framework. This is a direction we plan to explore in future work. Indeed, there have already been attempts to apply the spin foam model to cosmological settings [79].

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 (4.6), 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 ¹ 1 and ¹ 2 . We have successfully obtained the results, demonstrating the feasibility of studying the BH-to-WH transition numerically using the complex critical point method. As the next step, these calculations need to be repeated for various parameter sets. Although the results are expected to be exponentially suppressed, their relative values will provide deeper insights into the BHto-WH transition within the framework of the spin foam model. Addressing these issues will be part of our future work. The primary obstacle in this study is computational power. For the current results, obtaining the two data points in (5.15) required approximately 20 days using Mathematica code. To improve computational efficiency, potential strategies such as implementing the code in other programming languages or utilizing more powerful computational resources should be considered.

where the action S f is

with Z vef = g ve z vf and ∥Z∥ = ïZ, Zð for Z ∈ C 2 .

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

In the spin foam model, the amplitude reads

The action S, in the EPRL model, takes the form (see, e.g., [74,76] and Appendix A)

where S f and S ′ f are functions depending merely on the SL(2, C) element g and the spinor z. Applying the Poisson summation formula for the amplitude [77], one has

where the focus is on the phase in the exponential:

To study the asymptotic behavior of A for t → 0, we rescale x ′ by introducing x = tx ′ so that we get

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

The critical point is the solution to the equation of motion ¶S eff = 0. The properties of the critical point have been well studied in cases where the 2-complex contains no boundary (see, e.g., [74,75]). We can apply the results obtained from these studies, with the only exception being the treatment of the boundary faces carrying the coherent state (A2). Therefore, our attention should be directed towards the effective action associated with a boundary face f

Then, (B6) is a result of the replacement j f → x f /t. Treating j f as an integration variable, we need to study the equation ∂ j f S = 0, indicating that both its real and imaginary part are 0. Starting with its real part, we get

We impose the following decomposition for g ℓ

where n e ℓ , n t ℓ ∈ SU(2), p, ϕ ∈ R. Then, we could define new variables

and rewrite (B8) as

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 È t g takes [67, 85] ∥È t g ∥ = 2 √ Ãe t/8 |p|e |p| 2 2t t 3/4 sinh(|p|) . (B12) Therefore, for effective action corresponding to the transition amplitude ïA f |È t g ð/∥È t g ∥, the condition of vanishing its real part results in internal edges e Re ln ïZ seef , Z teef ð ∥Z seef , ∥∥Z teef ∥ + Re ln ï Zse n enf , e (ip-ϕ)Ä3 Zte 0 e0f ð ∥ Zse n enf ∥∥ Zte 0 e0f ∥ -|p| 2 4j f t -t 4 j f = 0 (B13) Comparing this equation with (B11), we get -|p| 2 4j f t -t 4 j f = -t 2 j f , (B14) leading to the area matching condition j f = |p| t . (B15) Taking the area matching condition into account, one can demonstrate e -tj f Zse n enf ∥ Zse n enf ∥ , e (ip-ϕ)Ä3 Zse 0 e0f ∥ Zse 0 e0f ∥ 2 = sin ´0 2 sin ´n 2 e p-|p| 2 + cos ´0 2 cos ´n 2 e -|p|-p 2 2 e -|p| 1 2 sin(´0) sin(´n)[1cos(µ 0µ nϕ)]

Here, the parametrization is applied as follows

(B17)

The result (B16), combined with the inequality

∥Z se ef ∥ 2 ∥Z te ef ∥ 2 f 1, results in the solution of (B11),

(1, 0), p > 0, e i ³ 0 +µ 0 2 (0, 1), p < 0, Zse n snf Zse n snf = e i ³n -µn

(B18)

When (B18) is satisfied, we have

= e i(³0-³n)-sgn(p)(µ0-µn-ϕ) (B19)

where P(z) := z/|z| for z ∈ C. In addition, combining (B10) and (B18), we have

where Z(p) = (1, 0) T for p > 0 and Z(p) = (0, 1) T for p < 0.

The above discussion solves the equation ∂ j f S eff = 0 for a boundary face f . Let us summarize the results as follows: For a boundary face f carrying the coherent state (A2) labelled by g ℓ = n s ℓ e (-ip+ϕ)Ä3 n -1 t ℓ , the area matching condition j f = |p|/t leads to the following solution of ∂ j f S eff = 0,

0 f À s ℓ , Z se n enf ∥Z se n enf ∥ = e i³ en f À t ℓ , and, internal e∈∂f ´ln ∥Z seef ∥ ∥Z teef ∥ + internal edges e Im ln ïZ seef , Z teef ð ∥Z seef ∥∥Z teef ∥ + ´ln ∥Z se n enf ∥ ∥Z te 0 e0f ∥ + Im ln ïÀ s ℓ , Z te 0 e0f ð ∥À s ℓ ∥∥Z te 0 e0f ∥ + Im ln ïZ te n enf , À t ℓ ð ∥Z te n enf ∥∥À t ℓ ∥ + sgn(p) 2 ϕ = 0,

for ³ ef ∈ R and k ∈ Z, where we define

The solutions to the equations ∂ z S = 0 and ∂ 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 [74,75].

It is worth noting that the first three equations in (B21) can be recovered through the usuall model, where the boundary state is selected as

denote the value of the action taken at the critical points in that conventional model. Then, the last equation in (B21) can be expressed as

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

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

where in the last line we introduce the following convention: when e = e 0 , the term Z seef becomes À s ℓ , and when e = e n , the term Z teef becomes À t ℓ .

Appendix C: Bridge Continuous Geometry with Coherent States

Classical phase space

As shown in [75], the value of the spin foam action at the critical point takes the form

where ϵ = ± is a global sign on the two complex, sgn(V

with » = 8ÃG and ´being the Barbero-Immirzi parameter, where

L + 1 2 ϵ M N KL , and Ω KL Ã¹ is the curvature form of the sl(2, C) connection compatible with e I µ . Note that the integral in (C2) depends on the orientation of the coordinate system, unlike the case where the action takes |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 (C2). Then a minus sign should be added if a left hand coordinate system is chosen to do the integration.

Let M be diffeomorphic to R × Σ with Σ dnoting the spatial manifold. Performing the canonical analysis on the action (C2), we get the canonical pairs (A i a , E b j ), where A i a is the usual Ashtekar connection but E b j becomes

As a consequence, the dual of the flux, i.e., E b j ε bcd , is a 2-form. This is in contrast to the case where the flux is √ qe b j = |e|e b j , resulting in its dual to be a pseudo 2-form.

The Poisson bracket between the canonical pairs is

Due to the definition (C3) of E, the ¶-function in the equation is the one such that ¶ 3 (y 0 , y)d 3 y = 1 for a right hand coordinate system y a , and ¶ 3 (y 0 , y)d 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 i a (x), E b j (y)f b j (y)d 3 y} = »´f b j (x), where f b j is a smeared function.

2. From the continuous to the discrete phase spaces

The holonomy h ℓ along an oriented link ℓ : t → ℓ(t) with t ∈ [0, 1] is the solution to

For a link ℓ, let S ℓ be a 2-surface transversally intersecting ℓ at at p = ℓ(1/2). The covariant fluxes are defined as

where Ä s ℓ (x) is a path starting from the source of ℓ, traveling along ℓ until reaching p, and then running in S ℓ to reach x ∈ S ℓ , and Ä t ℓ (x) is a similar path but starting from the target of ℓ. In order to provide a clearer understanding of the definition in (C7), let us clarify how to do the integration. To begin with, a right hand coordinate system {x 0 , Ã 1 , Ã 2 } is chosen in a region containing the surface S ℓ . It is required that (Ã 1 , Ã 2 ) could be chosen as a coordinate on S ℓ . With this coordinate, n c is defined as the conormal with direction compatible with the orientation of S ℓ , satisfying n c (∂/∂x 0 ) c = ±1. Without loss of generality, the orientation of S ℓ is chosen to be compatible with ℓ in what follows. In the general case where the orientation of S ℓ is allowed to be incompatible with ℓ, one only needs the replacements p i s,ℓ → sgn( la n a )p i s,ℓ and p i t,ℓ → sgn( la n a )p i t,ℓ to ensure the validity of all subsequent conclusions. According to (C7), one has

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

where Ä k := -iÃ k /2 with Ã k being the Pauli matrix. In the quantum theory, a wave function È(A) is given by a cylindrical function

where µ is a graph embedded in Σ comprising links ℓ ∈ E(µ) and È µ is a complex function on SU(2) |E(µ)| with |E(µ)| denoting the number of elements in E(µ). The action of the flux operators on the wave functions are

Coherent state

In this section, we focus on a graph of a single link ℓ for our calculation. Given g ℓ ∈ SL(2, C) which is parametrized by

the coherent state labeled by g ℓ , as shown in (A2), is given by

Using these notions, we can decompose g ℓ as

where p ℓ ∈ R is given by |p ℓ | = ∥⃗ p ℓ ∥, and ⃗ p ′ ℓ as well as ϕ ℓ is defined by

It can be verified that

The geometry interpretation for the parameters ⃗ p ℓ , ⃗ p ′ ℓ and ϕ ℓ can be deduced from the expectation values of operators. We have [23]

In terms of the expection values, g ℓ is expressed as

4. Compute g ℓ from the intrinsic and extrinsic geometry

The Ashtekar connection can be written in terms of the spin connection Γ i a and the extrinsic curvature for K i a as

Given a link ℓ : [0, 1] ∋ t → ℓ t ∈ Σ, we consider G ℓt as the parallel transport of the spin connection, i.e., solution to

we could obtain the equation for V ℓt ,

Let the link ℓ be a geodesic with respect to the intrinsic geometry given by the triad e

Let n s and n s be given by (C14) for ⃗ v being -la 0 e i a (ℓ 0 )/∥ l0 ∥ and -la 1 e i a (ℓ 1 )/∥ l1 ∥ respectively. They satisfy

Taking advantage of n s and n t and applying (C25), we get that G ℓ takes the following form,

for some real number ¹ ℓ .

In order to calculate V ℓt , we assume that the extrinsic curvature K ab takes the form,

along the link ℓ, for a field À on ℓ. Indeed, (C28) is assumed because the extrinsic curvature in the Regge calculus takes the same form with a distributional field À. It is easy to verify that

where n a (t) = q ab (t) lb t /∥ lt ∥ with q ab = ¶ ij e i a e j b being the spatial metric. The assumption (C28) leads to

As a consequence, we get

This equation, together with (C23), results in V ℓ ≡ V ℓ1 as follows

where we have change the variable of integration to the curve length s. Combining the equations (C22), (C27) and (C32), we get

This expression for h ℓ ensures h ℓ e i a (ℓ 1 ) la

which shares a similar form as (C8), but has a different sign on the right hand side. Now let us come to the fluxes p i t,ℓ and p i s,ℓ associated with a surface S ℓ intersecting ℓ at x o = ℓ to . We introduce h ℓ[0,to] as the holonomy along the segment [ℓ 0 , x o ] ¢ ℓ, and h ℓ[to,1] as the one along [x o , ℓ 1 ] ¢ ℓ. By definition of the fluxes in (C7), we have

where Ä xo (Ã) is the path starting at x o and running in S ℓ till Ã ∈ S ℓ .

To align with our discussion in (C 3), we introduce p ℓ as we did in (C15). Specifically, p ℓ ∈ R satisfies

Let us choose S ℓ appropriately to make the following relation holds

Indeed, (C37) is inspired by the facts that for a general S ℓ , one has

and that one can choose a S ℓ ensuring

Then, we get

and

leading to

where (C33) is applied and 1ℓ is given by

In this expression, 1ℓ has the geometric meaning of the twist angle (see, e.g., [73] for the definition), and ℓ À(s)ds is related to the extrinsic curvature. Due to (C29), the integration for À(s) can be alternatively rewritten as where q is the determinant of the metric q ab . If the cylinder is choose such that the area of A s are the same for all s, (C47) becomes

where ∥A∥ = As √ qn a s ϵ abc dx b dx c denotes the area of arbitrary A s .

Geometric meaning of ℓ ξ(s)ds

Embed our spatial manifold Σ into a spacetime (M, g µ¿ ). Let n µ be the unit conormal field of Σ, i.e., n µ n µ = -1. We have

where q ¿ µ = ¶ ¿ µ + n µ n ¿ is the projection that projects tensors to their spatial part. Given a link ℓ lying in the spatial manifold Σ, we consider g µ ¿ ∈GL(3) as the parallel transport of the Levi-Civita connection Γ µ Ã¿ along ℓ. Namely, g µ ¿ (ℓ t ) is the solution to

This equation together with (C49) leads to

Let e µ I be the vielbein field in M which takes e µ 0 = -n µ and e µ i = ¶ µ a E a i / | det E| for i = 1, 2, 3 on the hypersurface Σ. The spin connection associated to e µ I is given by

The holonomy p I J of the spin connection É I µJ is given by

Along the link ℓ, we have

This equation motivates us to contract e ¿ I and p -1 with the both sides of (C51), resulting in

It is useful to recognize the geometric meaning of the right hand side. By definition, n ¿ (ℓ t )(g(ℓ t ) -1 ) ¿ ³ = ñ³ (t) is the vector obtained by parallel transporting n µ (ℓ t ) to the starting point ℓ 0 of ℓ. Then, ñ³ (t)e ³ J (ℓ 0 ) is the inner product of ñ³ (t) with the veilbein at ℓ 0 , giving the information of the direction of ñ(t) with respect to e ³ J (ℓ 0 ). Due to e ³ 0 (ℓ 0 ) = n ³ (ℓ 0 ), one has |ñ ³ (t)e ³ 0 (ℓ 0 )| = cosh the dihedral angle between ñ³ (t) and n ³ (ℓ 0 ) .

In the formula (C55), the holonomy p -1 is the one of the 4-D spin connection. To interpret the geometry meaning of ℓ À(s)ds, the 3-D connection and extrinsic curvature should be incorporated. leading to

where (C66) is applied to get the last equality. Combining the result (C68) with (C55) results in

where we have considered the initial condition that both sides takes value 1 2 at ℓ 0 . This equation implies

where we introduced the convention ñ³

t)e i ³ (ℓ 0 )Ã i ∥ñ ³ (t)e i ³ (ℓ 0 )∥ u(t) (C72) for some element u(t) ∈ SU(2). To relate p 1 to V ℓt given in (C23), we need the relation between Kµ lµ t and the extrinsic curvature K i a la t lµ t Kµ = 1 2 lµ t É l0 µ Ã l = i la t K l a Ä l . (C73) Since we have identified p 2 with the holonomy G ℓt , the equation for p 1 now becomes ∂ t p 1 (ℓ t ) = -ip 1 (ℓ t )G ℓt K a (ℓ t ) la t G -1 ℓt . (C74) Focusing on the special case considered in (C28) and applying (C31), we could get p 1 (ℓ t ) = expℓ[0,t] À(s)ds la 0 e i a (ℓ 0 ) ∥ l0 ∥ Ã i 2 , (C75) where ℓ[0, t] ¢ ℓ denotes the segment of ℓ with the parameter taking values in [0, t]. Comparing this result with (C72), we get u(t) = 1 and ℓ[0,t] À(s)ds = arccosh(-n ³ (ℓ 0 )ñ ³ (t)), la 0 e i a (ℓ 0 ) ∥ l0 ∥ = ñ³ (t)e i ³ (ℓ 0 ) ∥ñ ³ (t)e i ³ (ℓ 0 )∥ , or ℓ[0,t]

À(s)ds = -arccosh(-n ³ (ℓ 0 )ñ ³ (t)), la 0 e i a (ℓ 0 ) ∥ l0 ∥ = -ñ³ (t)e i ³ (ℓ 0 ) ∥ñ ³ (t)e i ³ (ℓ 0 )∥ .

Compatibility between the spin foam action and the coherent state

Consider a 2-complex K and the spin foam model on it. Given a non-degenerate critical point, one could construct a flat simplical geometry dual to the 2-complex [75]. In this dual simplical geometry, each 4-simplex is convex. Given a boundary triangle f surrounded by tetrahedra Ä 1 , • • • , Ä n in the dual simplicial geometry, the spin foam action associated with f takes the following value at the critical point [75],

where ϵ = ± and sgn(V 4 ) = ± are some sign factor, and Θ k,k+1 are the dihedral angle between Ä k and Ä k+1 . Since Ä k and Ä k+1 share a 4-simplex denoted by T k,k+1 , the dihedral angle between them can be well-defined as the Lorentzian angle between the outwards normals of Ä k and Ä k+1 . In addition, since the dual geometry is flat, we could also calculate the dihedral angle, denoted by Θ f , between the tetrahedral Ä 1 and Ä n , by treating Ä 1 ∪ Ä n as the boundary of the simplicial complex T 1,2 ∪ T 2,3 ∪ • • • ∪ T n-1,n . Let us investigate the relation between the term k Θ k,k+1 in (C78) and the dihedral angle Θ f between Ä 1 and Ä n . Since the dual geometry is flat, the normals of the tetrahedra are coplanar. The plane containing all of the normals are referred to as the normal plane. The normal plane intersects the triangle f at a point, denoted by O. When this normal plane intersects the tetrahedron Ä k , it generates a line originating from O. In Fig. 10, we illustrate a general situation of the intersection between the normal plane and the simplicial complex. Here, we assume that the intersection fills the region beneath Ä 1 ∪ Ä k .

The direction of the line corresponding to the tetrahedron Ä i is specified by the coordinates (t i , x i ), given by (t i , x i ) = (sinh(¹ i ), cosh(¹ i )), i f m (sinh(¹ i ),cosh(¹ i )), i > m (C79)

Thus the future pointing unit normal vector of Ä i is v i = (cosh(¹ i ), sinh(¹ i )), i f m, (cosh(¹ i ),sinh(¹ i )), i > m.

(C80)

For the tetrahedra Ä k and Ä k+1 contained in T k,k+1 , their dihedral angle Θ k,k+1 is given by

leading to

Since the 4-simplex T m,m+1 shared by Ä m and Ä m+1 is convex, we have ¹ m + ¹ m+1 < 0. Therefore, (C83) becomes

Note that ¹ 1 + ¹ n f 0 implies that Ä 1 ∪ Ä n is convex, and ¹ 1 + ¹ n > 0, implies that Ä 1 ∪ Ä n is not convex. We thus get

In addition, Ä 1 ∪ Ä n , as the boundary of the simplicial complex, carries the extrinsic curvature of the form given in (C28), where À(t) becomes a distribution and the curve ℓ is considered to be the straight line perpendicular to the triangle f . In accordance with (C77), taking into account the geometric interpretation of Θ ℓ within that expression, we get The spin foam data in B -, which satisfies the boundary condition on ∂ -B, is replicated to obtain the data in B + . Consequently, the extended data in the entire region B remains compatible with the boundary condition on ∂B. Moreover, the replication ensures that the extended data solve the equations associated with all v, e, f except for those located on the transition surface T . These v, e, f that are not located on the transition surface T will be referred to as the interior objects. According to this discussion, we are allowed to focus only on the equations that involve objects crossing the transition surface T .

At first, ∂ j f S = 0 in (4.8) implies Re(S f ) = 0, where Re(S f ) denotes the real part of S f . As shown in literature [74], Re(S f ) = 0 is equivalent to

for some real number ³ f e , implying that Z vef is parallel to Z v ′ ef . Exploiting this result, we can simplify the last two equations in (4.8) to (see, e.g., [74] for the derivation)

Here, the first equation involves two edges e and e ′ satisfying e ∩ e ′ = v and e, e ′ ∈ ∂f . In the last equation, ϵ ef (v) is define by

Since the equations in (D2) do not involve objects crossing T , their validity can be directly derived from the fact that the extended data solve the equations associated with the interior objects. Thus, we only need to prove the validity of ∂ j f S = 0 for all faces f crossing T . By definition (2.16), we have

e∈∂f ln ïZ seef , Z teef ð 2 ïZ seef , Z seef ðïZ teef , Z teef ð + iµ ln ïZ seef , Z seef ð ïZ teef , Z teef ð . (D4) The edges e ∈ ∂f comprise of those entirely lie in B ± and these crossing T . For edges that entirely lie in B -, we apply (D1) to get ln ïZ seef , Z teef ð 2 ïZ seef , Z seef ðïZ teef , Z teef ð = ln e i³ f e ln ïZ seef , Z seef ð ïZ teef , Z teef ð = ln ∥Z seef ∥ 2 ∥Z teef ∥ 2 . (D5) For the edge e ′ = r(e) that entirely lie in B + , the compatibility of the orientation of the face f leads to s e ′ = r(t e ) and t e ′ = r(t e ). This fact, together with (4.10) and (D5), results in ln ïZ s e ′ e ′ f , Z t e ′ e ′ f ð 2 ïZ s e ′ e ′ f , Z s e ′ e ′ f ðïZ t e ′ e ′ f , Z t e ′ e ′ f ð = ln e -i³ f e ln ïZ s e ′ e ′ f , Z s e ′ e ′ f ð ïZ t e ′ e ′ f , Z t e ′ e ′ f ð = ln ∥Z teef ∥ ∥Z seef ∥ . Substituting (D7) and (D8) into (D4), we finally get ∂ j f S = 0 for f crossing T .

According to this two equations, ∆ 2 approaches ∞ as ϵ = 0, implying that when ϵ = 0 one of the adjacent tetrahedra with dihedral angle ∆ 2 becomes null. Moreover, the boundary ∂ -B 1 and ∂ -B 2 are expected to be spacelike, leading to

, and

Taking advantage of (E1) and (E4) and considering t 1 = 0 in the construction, we finally obtain

ϵ (r 3 -3r 1 ), t 2 = 1 3 (r 2 -3r 1 ) - ¶ (r 3 -3r 1 ) (E5) with ϵ and ¶ satisfying 0 < ¶ < ϵ.

(E6)

Since S 3 has to be in the future of S 2 , we need t 3 > t 2 , which leads to ϵ - ¶ < r 3r 2 3(r 3 -3r 1 )

. (E7)

Here, we use the fact that r 3 > 3r 1 for a large black hole.

According to (3.20), the values of a 1 and a 2 can be written in terms of ¶ and ϵ as

The values of ¹ 1 and ¹ 2 , according to (3.10) and (3.11), can also be related to ¶ and ϵ. With this relation and the inequalities (E6) and (E7), we finally get 2 f (r 1 ) -1 < cosh(¹ 1 ) < 2 f (r 1 ) 1 -r 2 -3r 1 3(r 2r 1 )

where we used f (r 3 ), f (r 1 ) > 0.