Improved Effective Dynamics of Loop-Quantum-Gravity Black Hole and Nariai Limit

We propose a new model of the spherical symmetric quantum black hole in the reduced phase space formulation. We deparametrize gravity by coupling to the Gaussian dust which provides the material coordinates. The foliation by dust coordinates covers both the interior and exterior of the black hole. After the spherical symmetry reduction, our model is a 1+1 dimensional field theory containing infinitely many degrees of freedom. The effective dynamics of the quantum black hole is generated by an improved physical Hamiltonian ${\bf H}_\Delta$ where the holonomy correction is implemented by the $\bar{\mu}$-scheme regularization with a Planckian area scale $\Delta$. The effective dynamics recovers the semiclassical Schwarzschild geometry at low curvature regime and resolves the black hole singularity with Planckian curvature. Our model predicts that the evolution of the black hole at late time reaches the charged Nariai geometry ${\rm dS}_2\times S^2$ with Planckian radii. The Nariai geometry is stable under linear perturbations but may be unstable by nonperturbative quantum effects. Our model suggests the existence of quantum tunneling of the Nariai geometry and a scenario of black-hole-to-white-hole transition. During the transition, the linear perturbations exhibit chaotic dynamics with Lyapunov exponent $\lambda=2\pi T_{\rm dS}\sim \Delta^{-1/2}$ relating to the Hawking temperature $T_{\rm dS}$ of ${\rm dS}_2$. In addition, the Nariai geometry in our model provides an interesting example of Wheeler's bag of gold and contains infinitely many infrared soft modes with zero energy density. These infrared modes span a Hilbert space carrying a representation of 1-dimensional spatial diffeomorphisms (or the Witt/Virasoro algebra) which are conserved charges of the effective dynamics by ${\bf H}_\Delta$.


Introduction
The research on quantum black holes has recently made important progress in the framework of loop quantum gravity (LQG) (see e.g. [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17], see also [18] for a recent review). LQG, as a background-independent and non-perturbative approach of quantum gravity, has the advantage for studying non-perturbative quantum effects in strong gravitational fields such as inside the black hole or the beginning of the Universe. In particular, LQG leads to the success of resolving both black hole and big-bang singularities (see e.g. [19,20] for cosmology). In both cases, the classical curvature singularities are replaced by the non-singular bounces where their curvatures are finite and Planckian.
These developments of LQG black hole are based on symmetry reduced models. Instead of the full quantum theory of gravity, these models of black hole quantize spherical symmetric gravitational degrees of freedom (DOFs). They fall into two categories: the models in the first category e.g. [1-5, 10, 11, 13, 15, 21, 22] quantizes the black hole interior and exterior separately using the homogeneous Kantowski-Sachs slices, and they only quantize a finite number of DOFs due to the spherical symmetry and homogeneity. The models in the second category e.g. [6,7,14,16,17,23] only reduce the DOFs by spherical symmetry, so these models are 1 + 1 dimensional field theories which still contain infinitely many DOFs, and they can treat both the black hole interior and exterior in a unified manner. Our approach in this work belongs to the second category.
The effective dynamics is one of important tools for studying black holes in LQG. The effective Hamiltonian of the black hole modifies the classical Hamiltonian in the canonical spherical symmetric gravity by including holonomy corrections to the Ashtekar-Barbero connection A, motivated by the LQG quantization. It is shown that in loop quantum cosmology (LQC), the effective dynamics is able to accurately capture key features of the quantum dynamics [24]. It is expected that the effective dynamics should also be valid for black hole, or at least be a crucial first step toward understanding full quantum effects of black hole. The effective dynamics from all the existing models lead to the black hole singularity resolution and bounce, but their detailed behaviors are different and depend on their schemes of the holonomy corrections. The choice of schemes are involved in these models due to the regularization/discretization of the Hamiltonian in which the curvature of A is discretized by the loop holonomy around a plaquette.
Here are a few examples of schemes used in these models: among models in the first category, the earliest models [1,2] follow the strategy similar to the μ 0 -scheme in LQC, and use plaquettes with constant fiducial scale. Their resulting behaviors of the bounce depend on the fiducial scale of plaquette. The model in [3] follows the similar strategy of the improvedμ-scheme of LQC, where the scale of the plaquette is dynamical (see also [4]). This model removes the dependence on the fiducial scale, and interestingly relates the final state of black hole to the Nariai geometry dS 2 × S 2 (similar in [4], see also e.g. [25] for earlier studies of the Nariai solution). But the model [3] has two problems (1) it gives large quantum effect near the event horizon where the spacetime curvature is small (this may be the consequence from choosing the Kantowski-Sachs foliation whose spatial slice approaching null near the horizon), and (2) the area of two-sphere becomes too small in the evolution, making the model self-inconsistent with theμ-type regularization (see [11] for the discussion). The more recent model [11] applies a new regularization scheme. In this scheme, the fiducial scale is a conserved Dirac observable and remains constant along the trajectory of evolution, but it may change value along different trajectories. The models [6,7,16,17] in the second category applies theμ-type regularization, complemented by certain gauge fixing and choice of lapse function. The foliations of these models are not Kantowski-Sachs but cover both the black hole interior and exterior.
In this work, we propose a new model of the spherical symmetric LQG black hole and study the effective dynamics. The model is embedded in the framework of the reduced phase formulation. We deparametrize gravity by coupling it to the Gaussian dust [26,27]. The dust fields provide a material reference frame of the time and space. The dynamics in the reduced phase space is governed by a physical Hamiltonian, which in our case is identical to the Hamiltonian constraint with unit lapse. When reducing to the sector of spherical symmetric DOFs, we obtain an 1 + 1 dimensional Hamiltonian field theory describing the spherical symmetric gravity-dust system. Classical solutions of this theory are Lemaitre-Tolman-Bondi spacetimes (classically, our model closely relates to the earlier spherical symmetric reduced phase space model [28], see also [16,[29][30][31][32] for recent works on coupling dust to black hole). We prefer the reduced phase space formulation because the dynamics is generated by the Hamiltonian which can make sense the unitarity when we pass to the quantum theory.
Our model applies theμ-scheme regularization (following [6]) to include the holonomy correction to the physical Hamiltonian. The improved effective Hamiltonian H Δ depends on the Planckian area scale Δ ∼ P which usually is set to be the minimal LQG area gap. The improved effective dynamics is given by solving the Hamiltonian equations of motion (EOMs) of H Δ . The EOMs are a set of partial differential equations (PDEs) which we call the effective EOMs. The corrections of O(Δ) are understood as quantum corrections to the classical EOMs. As a reason of choosing theμ-scheme regularization, it has the nice properties that the improved effective dynamics from H Δ has infinitely many conserved charges corresponding to spatial diffeomorphisms, which play an interesting role in our discussion.
An important advantage of our model is that it includes infinitely many DOFs and describes the effective dynamics of both interior and exterior of the black hole in one set of effective EOMs. Indeed, in the case of zero dust density and classical limit Δ → 0, the solution of effective EOMs is the Schwarzschild geometry in Lemaitre coordinate, which covers both the interior and exterior of the black hole. Classically the spatial slice of Lemaitre coordinate starts from the spatial infinity, crosses the event horizon, and ends at the singularity (the slice is further extended to another infinity by the singularity resolution in our model). Our model based on the reduced phase formulation is different from other ones in the second category. H Δ corresponds to the unit lapse function, and does not use the areal gauge fixing, as other two differences from [7,16,17].
When solving the effective EOMs of H Δ , we mainly focus on the situation with small dust density ρ. We firstly derive non-perturbatively the vacuum solution in the case of negligible ρ, then turn on nontrivial ρ by including linear perturbations. The vacuum solution is obtained by the numerical method with high precision. Some features of this solution are highlighted below, while the details is given in section 5.
(a) The solution satisfies the semiclassical boundary condition, i.e. reduces to classical Schwarzschild geometry in Lemaitre coordinate near spatial infinity. The quantum correction are negligible in the low curvature regime. (b) The black hole singularity is resolved and replaced by the non-singular bounce of the spatial volume element. The Kretschmann scalar K ∼ Δ −2 at the bounce is Planckian as Δ ∼ 2 P . The Lemaitre spatial slice, classically ending at the singularity, is now further extended to another infinity. (c) After the bounce, the evolution stabilizes asymptotically to the charged Nariai limit dS 2 × S 2 with different dS and S 2 radii. It is interesting that we obtain the Nariai limit similar to the result in [3], although our model is very different from theirs (the model in [3] belongs to the first category while ours belongs to the second category). (d) Although we have the similar result as in [3], our model is free of its problems: thanks to the foliation used in our approach, the geometry near the event horizon is semiclassical with negligible quantum correction. The area of two-sphere never becomes smaller than Δ in the evolution, consistent with theμ-scheme regularization.
The charged Nariai limit dS 2 × S 2 in our model is due to quantum effect. Both radii of dS and S 2 are of O( √ Δ). The appearance of the Nariai geometry is an interesting feature of the model, given that extensive studies in 90s suggesting the relation between the Nariai geometry and the black hole remanent e.g. [25,[33][34][35]. As a difference from early results on the Nariai limit, in our model the Nariai geometry dS 2 × S 2 is stable when we turn on linear perturbations. This result is similar to [3,36]. But we find evidence that the Nariai geometry is unstable by taking into account the non-perturbative quantum effect. The quantum tunneling effect may send dS 2 × S 2 to dS 2 × S 2 with opposite time orientation. Then following the effective EOMs, dS 2 × S 2 decays to the Schwarzschild white hole spacetime with complete future timelike and null infinities. The solution evolving from dS 2 × S 2 to the white hole is the time reversal of the vacuum solution discussed above. The entire picture taking into account the black hole evaporation and the quantum tunneling proposes a scenario of black-hole-to-white-hole transition similar to [8,9,37]. The detailed discussions of the black-hole-to-white-hole transition and evidences of quantum tunneling is given in sections 6-8. There are another two interesting aspects that are analyzed heuristically with the effective approach in this paper, and are discussed as parts of the outlook (see section 9): firstly, when we turn on linear perturbations on dS 2 × S 2 , the perturbations exhibit chaotic dynamics where we can extract the Lyapunov exponent λ = 2πT dS relating to the Hawking temperature T dS at the cosmological horizon of dS 2 . Moreover λ ∼ Δ −1/2 suggests that this chaos is due to the quantum gravity effect. The relation between λ and T dS resembles the known relation in the black hole butterfly effect [38,39].
Secondly, we suggest that the Nariai limit dS 2 × S 2 is an example of Wheeler's bag of gold. The foliation corresponds to the inflationary coordinate of dS 2 and gives large space volume behind the event horizon of small area at the late time of the evaporation. When we turn on perturbations, we find that all perturbations become infinitely many infrared modes with zero energy density in dS 2 × S 2 . It is also the reason why dS 2 × S 2 geometry is perturbatively stable. This result seems to suggests that dS 2 × S 2 should have quantum degeneracy, and all the infrared modes should span a Hilbert space H dS 2 ×S 2 . In the reduced phase space formulation, the spatial diffeomorphisms are the symmetries of the physical Hamiltonian H Δ and give infinitely many conserved charges. Then H dS 2 ×S 2 is the representation space of the spatial diffeomorphism group Diff(S 1 ) on the space of dS 2 . The Lie algebra of Diff(S 1 ) is the Witt algebra, or the Virasoro algebra if the central extension is considered.
Here is the organization of this paper: section 2 discusses the preliminaries including the reduced phase space formulation with Gaussian dust and the spherical symmetric reduction. Section 3 constructs the improved effective Hamiltonian by theμ-scheme regularization. Section 4 discusses the effective EOMs. Section 5 discusses the numerical black hole solution of the effective EOMs and the Nariai limit. Section 6 takes into account the black hole evaporation and motivates the extension of the effective spacetime. Section 7 discusses the scenario of the black-hole-to-white-hole transition. Section 8 discusses the evidence of quantum tunneling from dS 2 × S 2 to dS 2 × S 2 . Section 9 discusses some interesting future perspectives and open questions.

Deparametrized gravity with Gaussian dust
The reduced phase space formulation couples gravity to clock fields at classical level. In this paper, we mainly focus on the scenario of gravity coupled to Gaussian dust [26,27]. The action is given by where S GR is the host action of gravity [40] S GR e μ I , where the tetrad e μ I determines the four-metric by g μν = η IJ e μ I e ν J , and Ω IJ μν is the curvature of the so(1, 3) connection ω IJ μ . β is the Barbero-Immirzi parameter. S GD is the action of the Gaussian dust: where T, S j=1,2,3 are clock fields and defines time and space coordinates in the dust reference frame. ρ, W j are Lagrange multipliers. The energy-momentum tensor of the Gaussian dust is which indicates that ρ is the energy density and W μ relates to the heat-flow [26]. We assume M R × Σ and make Legendre transform of dust variables: where q αβ (α, β = 1, 2, 3) is the three-metric and L n denotes the Lie derivative along the normal to the hypersurface Σ. The constraint analysis [26,27] results in Hamiltonian and diffeomorphism constraints C tot , C tot α which are first-class constraints, and eight second-class constraints z, z j , ζ 1 , ζ 2 , s, K: where T ,α ≡ ∂ α T. Solving second-class constraints gives by a choice of sign in the ratio between W j , P j . These relations simplifies C tot , C tot α to equivalent forms: where C, C a are pure gravity Hamiltonian and diffeomorphism constraints from S GR . S GR leads to gravity canonical variables A a α (y), E α a (y), where A a α (y) is the Ashtekar-Barbero connection and E α a (y) = √ det qe α a (y) is the densitized triad. a = 1, 2, 3 is the Lie algebra index of su (2). Dirac observables are constructed relationally by parametrizing (A, E) with values of dust fields T(y) ≡ t, S j (y) ≡ σ j , i.e. A a j (σ, t) = A a j (y)| T(y)≡t, S j (y)≡σ j and E j a (σ, t) = E j a (y)| T(y)≡t, S j (y)≡σ j , where σ, t are physical space and time coordinates in the dust reference frame. Here j = 1, 2, 3 is the dust coordinate index (e.g. A j (y) = A α (y)∂y α /∂σ j ). A a j (σ, t), E j a (σ, t) depending only on values of dust fields are independent of gauge choices of coordinates y. They are proven to be invariant (on the constraint surface) under gauge transformations generated by diffeomorphism and Hamiltonian constraints [41][42][43]. Moreover A a j (σ, t), E j a (σ, t) satisfy the standard Poisson bracket in the dust frame: where β is the Barbero-Immirzi parameter and κ = 16πG. A a j (σ, t), E j a (σ, t) are the conjugate pair in the reduced phase space P red .
The evolution in physical time t is generated by the physical Hamiltonian H 0 given by integrating h on the constant T(y) = t slice S. The constant T slice S is coordinated by the value of dust scalars S j = σ j thus is referred to as the dust space [27,43]. T ,α = 0 on S leads to (2.14) H 0 formally coincides with smearing the gravity Hamiltonian C with the unit lapse, while here C(σ) is in terms of Dirac observables A a j (σ), E j a (σ) and σ j=1,2,3 are dust coordinates on S. The dynamics is governed by for all functions f on the reduced phase space P red . The recent result in [44] leads to an understanding of dusts or other clock fields from the LQG point of view. The altitude is that when we quantize P red to obtain the physical Hilbert space and defineĤ as the quantization of H 0 , the quantum theory should be taken as the fundamental theory and starting point of discussions. Although the quantum theory is formally obtained by quantizing the classical theory, the classical theory is not fundamental but emergent from the fundamental quantum theory. From this point of view, the classical gravity-dust theory is not fundamental, but is the low-energy effective theory emergent from the quantum theory via the semiclassical approximation, as demonstrated in [44]. The quantum theory is defined byĤ and the physical Hilbert space. H Δ to be constructed in section 3 is expected to describe the quantum effective theory ofĤ, such that the gravity-dust theory discussed in this section is the low-energy effective theory.

Spherical symmetric reduction
We assume the dust space S R × S 2 and define the spherical coordinate σ = (x, θ, φ). When the reduced phase space P red is further reduced to the phase space P sph of spherical symmetric field configurations, the symplectic structure reduces to (2.16) where in the second step we fix the gauge such that E j a (σ) is diagonal while A a j (σ) is not (details of this gauge fixing and solving Gauss constraint can be found in [6,7,45]).
Here K x is 1/2 of the extrinsic curvature along x. The canonical conjugate pairs in P sph are [7] The physical Hamiltonian reduced to P sph is expressed as The time evolution of any observable f on P sph is given by d f /dt = { f , H 0 } where the Poisson bracket is from (2.16). Comparing δH 0 = S d 3 σδC(σ) to the standard pure gravity Hamiltonian δH GR = S d 3 σ(NδC + N j δC j ). The dust coordinates correspond to the lapse function N = 1 and zero shift. Any solution E x (x), E ϕ (x) of EOMs can construct a 4d metric by expressing in the dust frame (t, x, θ, ϕ) General solutions of the EOMs by H 0 give the Lemaître-Tolman-Bondi (LTB) metric [28] Λ with arbitrary real functions E(x), F (x). When E(x) > 0, it relates to the energy per unit mass of the dust particles. F (x) > 0 relates to the gravitational mass within the sphere at radius x. Solutions can be classified into three cases E(x) > 0, E(x) = 0, or E(x) < 0: where β(x) is an arbitrary function. t = β(x) is the singularity of the metric. The solution at E = 0 and constant F = R s = 2GM (β(x) = x) is the Schwarzschild metric in Lemaître coordinates.
It is necessary to discuss the boundary condition and boundary term in the physical Hamiltonian, since the dust space is noncompact. But we postpone this discussion to section 3.
The time evolution by H 0 has infinitely many conserved charges V(N) from spatial diffeomorphisms, 30) for all N(x) (vanishing at boundary). Moreover since the Poisson bracket {H 0 , C(x)} is proportional to C x (x), C(x) becomes infinitely many additional conserved charges when the initial value of the dynamics satisfies C x (x) = 0.

Improved Hamiltonian
Recall relations between K x , K ϕ and Ashtekar-Barbero connection: A 1 1 (σ) = 2βK x (x), A 2 2 (σ) = βK ϕ (x), A 3 3 (σ) = βK ϕ (x) sin(θ), we following [6,16,17] to modify K ϕ , K x (x) by including a deformation parameter Δ ∼ l 2 P = Gh which may be chosen as the minimal area gap in LQG: Clearly as Δ → 0, (3.1) reduces to K x . To motivate these modifications, let us construct SU(2) holonomies of A 1 1 , A 2 2 and A 3 3 along edges in x, θ, and ϕ directions 3 : Firstly let e be an edge along the x-axis toward the positive direction. We define the holonomy h Δ (A 1 1 ) along e as below, and change parametrization of e from x to u such that 3 A a j contains nonzero off-diagonal components (see (2.19)). The discussion here constructs holonomies only from the diagonal components of A a j . We define a different connection fieldÃ a j , such thatÃ a j = A a j at the diagonals whereas A a j = 0 at the off-diagonals. The holonomies h Δ are defined fromÃ a j . The holonomies h Δ only capture the components of A a j relating to K x , K ϕ while forgetting the off-diagonal components of A a j which are determined by E j a , because the triad should be regularized differently in the Hamiltonian. This regularization is not new and has been appeared in the literature [6,46].
where σ a=1,2,3 are Pauli matrices. If we set dx is fixed to √ Δ. Therefore h Δ (A 1 1 ) along the fixed-length edge e in the x-direction can be written as where we assume the fields are approximately constant along e whose length is √ Δ ∼ l P . In the result of h Δ (A 1 1 ), E x , E ϕ , K x are evaluated at the starting point of e. We define the holonomy h Δ (A 2 2 ) to be along an edge e toward the θ-direction. Again by changing parametrization of e from θ to s such that s ∈ [0, 1] → e(v) ∈ S, and noticing that If we let dθ ds = √ ΔR −1 = √ Δ/ |E x | ≡μ θ , the length of e is fixed to be √ Δ: Therefore the holonomy along a fixed-length edge e in the θ-direction can be written as where E x , E ϕ , K x are evaluated at the starting point of e. The holonomy h Δ (A 3 3 ) along ϕ direction can be derived similarly. We let dϕ ds = These fixed-edge-length holonomies are analogs of dressed holonomies used inμ-scheme in LQC [20,47]. Therefore the modification (3.1) is called the 'μ-scheme' of the spherical symmetric LQG. Note that the holonomies only capture components relating to K x , K ϕ in the Ashtekar-Barbero connection, but do not capture the spin-connection This choice of treatment is often referred to as the K-holonomy regularization, and has often been used in LQC and black holes [7,16,46,48].
The modification (3.1) can be obtained from regularizing the curvature where is the elementary plaquette with fixed area Δ [6,46]. The plaquette is in a cubic lattice adapted to the dust coordinate (x, θ, ϕ). Let coordinate lengths of edges of areμ x ,μ θ ,μ ϕ / sin(θ) (the edge along ϕ has constant θ) as defined above, they indeed give the fixed area to every plaquette : h Δ ( ) is made by four holonomies along edges with fixed length √ Δ: Then the improved physical Hamiltonian H 0 can be obtained by inserting the above h Δ ( ) in the regularized Euclidean : where e( jk ) denotes the area element on i j (see [47]). The terms independent of K in equation (2.23) are contributions from the spin-connection compatible to E j a , and are kept unchange. Comparing the above result to the terms depending on K in equation (2.23) justifies the replacement (3.1).
We define the improved Hamiltonian H Δ to be H 0 modified by the replacement (3.1): Entries of the above sine functions are restricted into a single period (−π, π]. H Δ is the same as (3.10). The dynamics generated by H Δ is referred to as the improved effective dynamics.

In deriving EOMs from H
and integrations by part result in boundary terms respectively We assume that E x , E ϕ behaves asymptotically the same as in the Lemaître coordinates of the Schwarzschild spacetime as x → ∞: . This boundary condition corresponds to the asymptotically flatness formulated in the dust coordinate, given that E x , E ϕ reduce to the Schwarzschild spacetime at infinity. We add the following boundary Hamiltonian to H ∞ cancel the boundary terms in (3.13) respectively, with the boundary condition (3.14). However, H ∞ vanishes at x → ∞ by (3.14), so we end up with zero boundary Hamiltonian at x → ∞. Alternatively, instead of the Schwarzschild boundary condition (3.14) and (3.15), we may make an infrared cut-off of the dust space at bdy = {x = L 1} and impose the Dirichlet boundary condition δE x | bdy = 0. In this case, we have to add the boundary term to the physical Hamiltonian up to a term proportional to δE x which vanishes by the Dirichlet boundary condition.
On the other side x → −∞, we impose the Neumann boundary condition so that both terms in (3.13) vanish at the other asymptotic boundary x → −∞. The reason to choose this boundary condition at x → −∞ will become clear when we analyze solutions in section 5.2.
As an advantage of the improved Hamiltonian H Δ , V(N x ) are still conserved by H Δ , i.e. for all N x (x) vanishing at boundary, where K x , K ϕ are not modified in C x . The computation of the bracket between V(N) and H Δ is given in [6,49]. If we consider the standard formulation of pure gravity in terms of constraints, and understand C Δ (x) as the improved Hamiltonian constraint, nicely resembles the Poisson bracket between classical Hamiltonian and diffeomorphism constraints [6]. However, the Poisson bracket between two Hamiltonian does not vanish when the diffeomorphism constraint is satisfied. The physical implication is that the resulting dynamics may depend on foliation of the spacetime. This issue is relieved in the reduced phase space formulation where the dust clock fields provide a specific foliation of the spacetime. C Δ (x) and C x (x) are not understood as constraints. Both Hamiltonian and diffeomorphism constraints have been resolved before we arrive and the classical counterpart only indicates that C Δ (x) are not anymore conserved in the dynamics when the initial value satisfies V(N x ) = 0, although its classical version C(x) is conserved. This may indicate certain quantum symmetry anomaly (rather not gauge anomaly) for special initial states. When the initial condition does not satisfy V(N x ) = 0, there is no symmetry anomaly since the classical C(x) is not anymore a conserved charge.
Note that C Δ (x) only takes into account the holonomy correction while neglecting the inverse triad correction. It was argued in [50] that at least for LQC, the effect from the inverse triad correction may be negligible in the effective dynamics on physical grounds.

Effective equations of motion
The signs of E x , E ϕ has to be fixed in order to derive the Hamiltonian EOMs from H Δ . Here we focus on two cases: both When both E x , E ϕ > 0, the Hamiltonian EOMs from H Δ are given below: H Δ and above EOMs reduces to H 0 and corresponding EOMs in the limit Δ → 0. We refer to the dynamics determined by above EOMs with nonzero Δ as the improved effective dynamics in analogy with LQC and various models of black holes [6,16,17,20,47].
When E x < 0, E ϕ > 0, H Δ gives a different set of EOMs. We express the coordinates to be (t,x,θ,φ) and write Ex < 0, Eφ > 0 to distinguish from the above case of both E x , E ϕ > 0. The EOMs in this case is given by Two set of EOMs (4.1)-(4.4) and (4.5)-(4.8) can be related by the spacetime inversioñ with the fields transforming as follows This transformation can be used as a solution generating map. Given a solution with e.g.
is a solution to EOMs in terms ofx,t. E x flips sign in (4.12) as x → −x can be understood from the definition E j a = √ det(q)e j a . In order to obtain special solutions, we impose the following ansatz to simplify the nonlinear PDEs (4.1)-(4.4) (or (4.5)-(4.8)), This ansatz is motivated by the Lemaître coordinate of the Schwarzschild spacetime. The existence of the solution satisfying the ansatz (see section 5) indicates the consistency between slicing the spacetime with z and the dust foliation. Here z parametrizes the spatial slice when fixing t, while parametrizing the time evolution if x is fixed. In the case Δ → 0, solutions from the ansatz corresponds to equations (2.26)-(2.28) with constant F , E and β = x. The ansatz reduces (4.1)-(4.4) to nonlinear ordinary differential equations (ODEs) of E x (z), E ϕ (z), K x (z), K ϕ (z). The resulting ODEs are 1st order in E ϕ (z), K x (z), K ϕ (z) and 2nd order in E x (z) (resulting from ∂ 2 x E x in equation (4.1)). But the 2nd order derivative in E x (z) can be eliminated by using Moreover for the convenience of solving equations numerically, we make a change of variable: As a result, the EOMs (4.1)-(4.4) reduce to a standard form of 1st order ODE: Similar change of variables and reduction can be applied to (4.5)-(4.8) and result in d dz The explicit expressions of (4.16) and (4.17) contain long formulae and are given in appendix A and [49]. Equation (4.17) can transform to equation (4.16) byz = −z and deduced from (4.9)-(4.13).

Strategies
When z 1, we are far from the classical singularity, K x , K ϕ are small thus (4.1)-(4.4) (or equation (4.16)) reduce to the classical EOMs from H 0 . We focus on a solution which reduces to the Schwarzschild spacetime when far away from the singularity, by imposing (2.27) as the initial condition at z = z 0 1, i.e.
With the above boundary condition, we can obtain a family of numerically solutions satisfying the ansatz (4.14) to EOMs (4.1)-(4.4) using Mathematica or Julia (for higher precision). The solutions are labelled by different values of parameters R s , Δ, z 0 . These solutions resolve the Schwarzschild black hole singularity with regular spacetime with finite but large curvature (see figures 1 and 2).
The strategy of finding numerical solutions to equation (4.16) is the following: since both E x and E ϕ are large in the semiclassical regime z 1, and the initial condition is at z 0 1, it is more efficient to change variables: and express EOMs in terms of Explicit expressions of (5.5) are given in [49]. v 1 , v 2 are not large in z 1, so less numerical errors are produced at the early stage of numerical evolution in z. In the z-evolution across z = 0 and toward −∞, E x is suppressed and stabilized at a nonzero constant value, while E ϕ grows exponentially as −z goes large (see figure 3).
In the region where E ϕ is exponentially large, we can expand EOMs (4.16) in 1/E ϕ , and neglect O(1/E ϕ ) to simplify the EOMs: where v 2 = ln E ϕ only appears in the equation (5.9). In practice, the numerical z-evolution from z 0 to z < 0 follows the EOMs (5.5) until certain z mid < 0 where E ϕ is too large. Then the solution E x , E ϕ , K 1 , K 2 of (5.5) evaluated at z mid serves as the initial condition for equations (5.6)-(5.9), which can be further evolved to arbitrarily large −z. The approximation in equations (5.6)-(5.9) by neglecting O(1/E ϕ ) is consistent because the solution E ϕ keeps growing exponentially as z → −∞, while all other quantities are bounded by O (1). A full solution of the EOMs is given by connecting two solutions of (5.5) and (5.6)-(5.9) at z mid .
Note that the similar approximation with Eφ 1 can be applied to (4.5)-(4.8) and leads to

Properties of solutions
The numerical solutions exhibit following properties: the spacetime curvature is finite on the entire range of z (from z 0 1 to negative z with arbitrarily large |z|). The classical singularity at z = 0 is resolved. The entire spacetime is nonsingular and has large but finite curvature at   Figure 1 plots the Kretschmann invariant K = R μνρσ R μνρσ of the solution. It demonstrates two groups of local maxima of K located respectively in the neighborhood N 0 of z = 0 ( figure 1(b)) and in a neighborhood N < of z < 0 ( figure 1(c)). The oscillatory K in these two neighborhoods indicates strong quantum fluctuations in these regimes. We denote by K max,0 and K max,< the maximal K in N 0 and N < respectively, and test their dependences with respect to Δ and horizon radius R s (we fix β = 1). The numerics demonstrate that both K max,0 and K max,< are proportional to Δ −2 (see figure 4) The behavior has qualitative similarity with results in [11,17]. The behavior of Kretschmann scalar K ∼ Δ −2 motivates us to understand Δ ∼ 2 P such that the singularity resolution and bounce of spatial volume (figure 2) happen at the Planckian curvature. Both K max,0 and K max,< are Planckian curvatures. In models of LQC and LQG black holes, Δ is chosen to be the minimal nonzero eigenvalue of the LQG area operator. Here K max,0 > K max,< , and k 0 , k < have mild dependence on R s and Δ (see figures 4 and 5). The R s and Δ dependences in k 0 , k < are subleading corrections. Asymptotically for large negative z, K approaches to be z-independent constant whose dependence on Δ is still ∼1/Δ 2 . We come back to this asymptotic behavior shortly. In the regime z > 0 and far away from N 0 , the solution is semiclassical and reduces to the Schwarzschild spacetime in Lemaître coordinates. The quantum effect is negligible in this regime. It is clear from EOMs (4.1)-(4.4) that as far as K x , K ϕ do not blow up, the classical Schwarzschild geometry is approximately a solution to (4.1)-(4.4) up to corrections of O( √ Δ). In particular, the numerical solution indicates that the quantum correction at the event horizon z = 2 3 R s is negligible. z = 2 3 R s is a marginal trapped surface with Θ k = 0 and Θ l < 0 where Θ k and Θ l are outward and inward null expansions (see figure 6).
The other asymptotic regime is in the opposite side where z < 0 and −z is large. In this regime E x approaches a constant, denoted by r 2 0 , and E ϕ grows exponentially: r 0 , α 0 , α 1 can be obtained numerically (see figure 7). E x (z) approaching to constant as z → −∞ (x → −∞ when fixing t) fulfills the boundary condition (3.18) discussed earlier.
Recall that the spacetime metric is given by (2.24), the simple asymptotic behavior (5.15) indicates that at z → −∞ is the metric looks like dS 2 × S 2 : a product of two-dimensional de Sitter (dS) spacetime and two-sphere, as indicated in figure 1. The asymptotical regime as z → −∞ still has Planckian curvature. α 1 depends on both Δ and R s (see figure 10). Although the asymptotic dS 2 × S 2 geometry is independent of R s or the mass of black hole, the dust coordinate x in the geometry depends on R s since α 1 depends on R s . dS 2 × S 2 with r 0 = α 0 is also known as the charged Nariai geometry. It can be obtained as the near horizon limit of the (near) extremal Reissner-Nordstrom-de Sitter (RN-dS) spacetime where the cosmological horizon and event horizon coincide [25,33]. The relation between LQG black hole and the Nariai geometry has been proposed in earlier studies of effective dynamics using the Kantowski-Sachs foliation of the black hole interior [3,36]. However as is pointed out in [11] that the analysis in [3,36] suffers from two problems: (1) theirμ-scheme model of black hole interior produce large quantum effect near the event horizon which is of low curvature, and (2) the area of S 2 is even smaller than the minimal area gap Δ at certain stage of the time evolution, inconsistent with theμ-scheme treatment of holonomies (S 2 has no room for theμ-scheme holonomy since it has to be around an area of Δ). Our analysis does not have these problems: firstly the quantum effect is negligible in the low curvature regime including the event horizon as discussed above. Secondly, numerical results indicate that the area of S 2 is always larger than Δ during the evolution (see figure 11).
The recent work [45] develops the polymer quantization of the phase space P sph and the Hamiltonian H Δ , and derive a path integral formula for the quantum dynamics of the model.   The effective dynamics (4.1)-(4.4) emerges from the stationary phase approximation of the path integral. It is shown in [45] that the effective dynamics discussed above is a good approximation of the quantum dynamics with the small quantum fluctuation, when the black hole is massive. But for the small black hole (e.g. the quantum black hole close to the end of evaporation), the quantum correction to the effective dynamics may not be negligible, so the Nariai geometry may be broken down by the quantum effective in this case (we will come back to this point in section 8). The path integral associates dS 2 × S 2 a quantum states |r 0 , α 0 ; α 1 such that when z → −∞: x , E ϕ have small (large) quantum fluctuations given by |r 0 , α 0 ; α 1 for the massive (small) black hole. It is interesting that these states |r 0 , α 0 ; α 1 depend on the additional parameter α 1 which indicates that dS 2 × S 2 is not a single state, but has infinite degeneracy from the quantum point of view. Although different values of α 1 correspond to the same spacetime geometry, and are related by diffeomorphisms, here they indeed correspond to different physical states since we start with the reduced phase space formulation where gravity is deparametrized by dust fields. We would like to emphasize that our system is not pure gravity but gravity coupled to dust. Different α 1 indeed correspond to the same spacetime geometry. But different α 1 correspond to different values of dust fields. It is because different α 1 corresponds to different coordinate systems (t, x) while the coordinate system is defined by values of dust fields. For constants α 1 = α 1 , the dust coordinate x changes as Different α 1 's correspond to inequivalent solutions of EOMs, so we expect that they correspond to different physical states. As another way to understand this: due to the deparametrization, the diffeomorphism in the dust space (the space of values of dust fields) is not anymore gauge symmetry but becomes the global symmetry (see e.g. [51]). These diffeomorphisms change physical states, because they change the coordinate system which is made by values of the dust fields.
All states |r 0 , α 0 ; α 1 corresponding to dS 2 × S 2 should span a Hilbert space H dS 2 ×S 2 . We come back to discussing more details of H dS 2 ×S 2 in section 9.  t, x, θ, ϕ) coordinate. S (black curve) is a typical spatial slice with constant t. Dashed curves are another spatial slice in the far past. The green line illustrate the event horizon, which bounds the black hole region. The gray triangular region has strong quantum fluctuation and has Plankian curvature. It enclose the classical singularity at z = x − t = 0 near its past boundary. At the future boundary of the patch, the asymptotic geometry is dS 2 × S 2 with Planckian radii.
Viewing dS 2 × S 2 to be the asymptotic geometry, we draw the Penrose diagram of the effective spacetime in figure 12. The Schwarzschild event horizon (green line) only receive negligible quantum correction (see figure 6 and the discussion there). The resulting spacetime has a complete future infinity since dS 2 × S 2 is complete. The dust time t can extend to t → ∞ in the spacetime as the inflationary coordinate in dS 2 . The Penrose diagram is similar to the one obtained in [52], although the authors there find asymptotically dS 4 in the black hole interior.
We compute the Einstein tensor of the solution and define the quantum effective stress-energy tensor T (eff) μν by R μν − 1 2 g μν R = 8πGT (eff) μν . We find T (eff) μν violate the average null energy condition, i.e. T (eff) UU dU and T (eff) VV dV are negative (see figure 13 for the null rays parametrized by U, V). In concrete, when parameters are e.g. z 0 = 3 × 10 8 , The main contributions to T (eff) UU dU, T (eff) VV dV are from the regions with local maxima of K. Here T (eff) μν does not correspond to any physical matter (the dust density is approximately zero in this solution), but rather is an effective account of the LQG effect in the black hole.
The numerical errors can be tested by inserting numerical solutions back into the EOMs. We find the EOMs are satisfied by solutions up to numerical errors which are bounded by ∼10 −21 with Julia and by ∼10 −8 with Mathematica [49].

Perturbation and stability
In this subsection, we exam the stability of the asymptotic geometry dS 2 × S 2 by turning on some perturbations. We still assume perturbations to satisfy the spherical symmetry, so we are going to linearize EOMs (4.1)-(4.4) at the asymptotic background dS 2 × S 2 . In practice, we make the change of variable from (E and insert in the EOMs the perturbation ansatz: where initial perturbations f 1 , . . . , f 4 are arbitrary functions of x. As t → ∞, asymptotically p 1 , c 2 , c 2 exponentially damp off while is controlled by the initial perturbation. We conclude that the asymptotic geometry dS 2 × S 2 is stable with linear perturbations. The perturbation of p 2 in (5.29) modifies α 1 in (5.15) and effectively defines a transformation in H dS 2 ×S 2 from |r 0 , α 0 ; α 1 to |r 0 , α 0 ; α 1 . We are going to come back to this point in section 9. Although dS 2 × S 2 is stable with linear perturbations, as we are going to see in section 8, this dS 2 × S 2 geometry may be unstable by non-perturbative quantum effect, and have nontrivial quantum transit by tunneling effect. We can extend the study of perturbations on the entire effective black hole spacetime. Since the background spacetime metric components only depend on z, we can rewrite the perturbations in terms of ξ = −z = t − x and x and make Fourier transformations of perturbations along x, e.g.

Picture of black hole evaporation
The quantum correction at the event horizon z = 2 3 R s is negligible in the black hole solution of effective EOMs. The geometry near and outside the horizon almost has no difference from the classical Schwarzschild spacetime. By turning on quantum field perturbations, the Hawking's original derivation of black hole evaporation happening near the horizon carries over to the black hole spacetime obtained here. The quantum correction from nonzero Δ to Hawking's derivation is negligible. The back-reaction from Hawking radiation reduces the black hole mass and causes the horizon to shrink. Then the event horizon should be replaced by the T-DH. A T-DH is a three-dimensional time-like submanifold foliated by two-dimensional surfaces h with two-sphere topology, so that at each leaf h, Θ k = 0 and Θ l < 0 where Θ k and Θ l are expansions of outward and inward null normals of h [53]. Moreover, for the semiclassical spacetime outside and far from the black hole, the future null infinity I + is extended until the 'last ray': the last Hawking particle radiated from the black hole before the evaporation stops. The picture of quantum effective black hole spacetime is illustrated in figure 15. The black hole evaporation results in the existence of the classical asymptotic flat regime which we call the region (I). All future null rays from points in the region (I) does not intersect with the black hole horizon, in contrast to the spacetime in figure 12 where the inward future null ray always cross the horizon.
Here we assume the evaporation process is sufficiently slow such that at every instance the spacetime can be approximately described by the solution of effective EOMs with a fixed R s . Foliating the spacetime with solutions with different R s approximates the dynamical black hole spacetime. This is known as the adiabatic approximation of the black hole evaporation.
As an advantage of solutions following the ansatz (4.14), when we can set a constant t = t 0 , describe the geometry on a constant t = t 0 spatial slices S. The two-sphere h ⊂ S given by x − t 0 = 2 3 R s is a marginal trapped surface with Θ k = 0 and Θ l < 0 where Θ k and Θ l are outward and inward null expansions (see figure 6). We may implement two different numerical solutions of different parameters R s (t 0 ), R s (t 1 ) at different spatial slices S t 0 , S t 1 at two instances t = t 0 , t 1 . The marginal trapped surfaces h(t 0 ) and h(t 1 ) are at x − t 0 = 2 3 R s (t 0 ) and x − t 1 = 2 3 R s (t 1 ). On S t 0 and S t 1 respectively. If the evaporation is slow enough, the dynamical black hole can be approximated by a large number of spatial slices S t carrying different solutions with different horizon radii R s (t) which is monotonically decrease as t growing. The T-DH is foliated by the set of marginal We run numerical experiments of solving effective EOMs (equation (5.5) by the ansatz (4.14)) with smaller and smaller R s (R s is implemented by the initial condition (5.1) and (5.2)). From the results, we find that the asymptotic dS 2 × S 2 geometry is invariant under changing R s , although details of curvature fluctuations and bounces between semiclassical Schwarzschild Figure 15. The picture of quantum effective black hole by taking into account the backreaction from Hawking radiation. The blue curve is the timelike trapping dynamical horizon (T-DH) of the evaporating black hole. The blue diamond illustrates the regime where the horizon area is comparable to Δ. This spacetime diagram is incomplete since the future of the last ray is missing. The completion is given in figure 19. and dS 2 × S 2 can change. In particular two neighborhoods N < and N 0 (of two groups of local maxima of K) become closer as R s decreasing, and finally merge when the horizon area 4πR 2 s is comparable to Δ (see figure 16). When N < and N 0 merge, the bounce near z = 0 looks like a domain wall separating the semiclassical (low curvature) Schwarzschild spacetime and the quantum (high curvature) dS 2 × S 2 spacetime.
Remarkably, when R s is small such that 4πR 2 s is comparable to Δ, e.g. in the case that R s = 0.19 in figure 16, the T-DH disappears and is replaced by a spacelike transition surface with both Θ k = 0 and Θ l = 0 (see figure 17). It is followed by a marginal anti-trapped null surface with Θ k = 0, Θ l > 0 on the left (at smaller x), then followed by anti-trapped and trapped regions (due to small fluctuations of geometry after the bounce, as seen in figure 16) before approaching to dS 2 × S 2 . Since the T-DH disappears, Hawking's derivation of black hole evaporation fails to valid in this regime. We expect the quantum gravity effect to be strong in the region of transition surface since the spacetime curvature at the surface is Planckian. The last ray of Hawking radiation should happen before this instance.
Note that the small R s does not break 4πE x > Δ, 2πE ϕ > Δ (see figure 18) as long as 4πR 2 s > Δ, 5 so our results are self-consistent with the starting point ofμ-scheme regularization. 5 The situation 4πR 2 s Δ is not tested since the evolution in this situation reaches the singularity of the ODEs. But 4πR 2 s Δ is inconsistent with theμ-scheme regularization. The region where the Hawking radiation stops is the blue diamond in figure 15. This region is a neighborhood of the transition surface in figure 17, and has the strong quantum gravity effect [18,54] since the curvature is Planckian. This region should also be strongly dynamical. Our analysis of the effective dynamics and approximation using foliations of solutions with R s (t) should be not adequate although it provides a preliminary picture of this phase. The rigorous analysis of this region should apply the full theory of LQG. The study on this aspect is beyond the scope of the present paper.

Black hole to white hole transition
The black hole spacetime in figure 15 is incomplete due to the existence of the region (I). The future null infinity I + should be extended beyond the last ray of Hawking radiation (see [5] for an earlier study of extension in 2d dilaton black hole). We are going to use the effective EOMs of H Δ to derive the extension. We make the following assumptions of region (I) in order to obtain boundary conditions for the effective EOMs: (a) The quantum dynamics in region (I) near I + can be well approximated by semiclassical spacetime geometry. The dynamical effect is weak. The spacetime near I + is Schwarzschild with different horizon radii, i.e. the adiabatic approximation is assumed to hold. R (0) s = 2M 0 is the initial Schwarzschild radius, and R (r) s = 2M r is the remnant black hole radius before Hawking radiation stops. Figure 19(a), as the analog of figure 7.2 in [55], is the conformal diagram for the black hole evaporation: in a neighborhood of the contour labelled by M(t) (or M r ), the exterior spacetime is well-approximated by the Schwarzschild geometry with R (t) s = 2M(t) (or R (r) s = 2M r ).  (b) For all spatial slices region (I), their asymptotic geometries (metrics and extrinsic curvatures) near the spatial infinity are classical and asymptotically flat. Their geometries are continuous extensions from geometries in the past. Namely there exists a slice St 0 in region (I) such that the asymptotic spacetime geometry of the neighborhood of St 0 near the infinity should be the same as the asymptotic geometry of the neighborhood of S t 0 in the close past of the blue diamond (see figure 19(a)). In particular, their horizon radii are approximately the same.
We find above assumptions are physically reasonable and should be an excellent approximation of the full quantum dynamics. from the Schwarzschild spacetime: we assume that the dynamics near I + (enclosed by the light blue dashed circle) can be approximated by Schwarzschild geometries with different radii. The dashed contour is the spatial slice S t . Below the contour labelled by M 0 the spacetime exterior to the black hole is well-approximated by the Schwarzschild geometry with the initial radius R (0) s = 2M 0 . In a neighborhood of the contour labelled by M(t) (or M r ), the exterior spacetime is well-approximated by the Schwarzschild geometry with R (t) s = 2M(t) (or R (r) s = 2M r ). The spacetime geometry around St 0 near spatial infinity is the same as the geometry around the slice S t 0 in the close past of the blue diamond; (b) an extended effective spacetime geometry beyond the last ray suggests the black-hole-to-white-hole transition. There is a dS 2 × S 2 geometry (denoted by dS 2 × S 2 ) behind the white hole event horizon (green line). The gray region between dS 2 × S 2 and dS 2 × S 2 contains a quantum tunneling. The entire spacetime includes two coordinate patches (t, x, θ, ϕ) and (t,x,θ,φ). Now we focus on the asymptotic geometry of region (I) near the spatial infinity i 0 of St 0 : the neighborhood of St 0 near infinity has the semiclassical Schwarzschild geometry with R (r) s = 2M r obtained from extending the geometry from the past, i.e. from S t 0 . We denote R s ≡ R (r) s for simplicity in the rest of the section. For simplicity of the present model, we do not consider to vary the Schwarzschild radius in the future of St 0 . We may draw infinitely many spatial slices St in region (I) such that the neighborhoods of all these slices carry the same semiclassical Schwarzschild spacetime geometry near the infinity.
The foliation (t, x, θ, ϕ) cannot be extend to region (I) due to the strong quantum dynamical effect in the blue diamond region, thus a new foliation by (t,x,θ,φ) is necessary for studying the dynamics of the spacetime including St 0 and its future. We impose the following boundary conditions near i 0 corresponding toz =x −t → −∞ for slices in region (I): where R s equals 2M r . We are going to apply the boundary conditions (7.1) and (7.2) to the effective EOMs of H Δ . The resulting solution extends the effective spacetime beyond the last ray. The boundary conditions (7.1) and (7.2) are obtained as follows: in region (I) near spatial infinity, the foliation St is obtained from S t by the time reflection symmetry T → −T of the Schwarzschild-Kruskal geometry (see figure 20 for illustration of spatial slices S t and St when they are in the classical Schwarzschild spacetime). We define diffeomorphism R maps S t to St, and relates the coordinate basis by 6 for any point p ∈ St. In terms of coordinates, we obtain that as functions, , which give (7.1). Since R leaves the 4d metric invariant, the four-metric in (t,x,θ,φ) is given by The boundary conditions (7.2) are given by solving the classical EOMs of H 0 (in the coordinates (t,x,θ,φ)) with Ex(t,x), Eφ(t,x) given by (7.1), or equivalently, they can be obtained by computing the extrinsic curvature of St. 7 The extension of black hole exterior geometry in region (I) can be viewed as the geometry of white hole exterior, as illustrated by figure 20. This aspect is similar to [54].
In order to extend the geometry beyond the last ray from region (I) to the causal future of black hole interior, we apply the effective EOMs of the improved Hamiltonian H Δ and implement the boundary conditions (7.1) and (7.2). The effective EOMs are given by equations (4.5)-(4.8) since we have E x < 0 by the boundary condition (due to equation (7.3)).
Recall that the EOMs (4.5)-(4.8) in (t,x,θ,φ) relate to the EOMs (4.1)-(4.4) in (t, x, θ, ϕ) by identifyingx = −x,t = −t (thus z →z = −z) and (4.10)-(4.13) which consistently match the relation between the boundary conditions (7.1) and (7.2) of Ex, Eφ, Kx, Kφ atz → −∞ and (5.1) and (5.2) of E x , E ϕ , K x , K ϕ at z → ∞. We again apply the ansatz Ex(t,x) = Ex(z), Eφ(t,x) = Eφ(z), Kx(t,x) = Kx(z), Kφ(t,x) = Kφ(z) to reduce the PDEs to ODEs 6 We denote the standard Schwarzschild coordinate by (τ , r, θ, φ). Its transformation to Lemaitre coordinate is dt = dτ + rs r 1 − rs r −1 dr, dx = dτ +  (equation (4.17)), so that the solution is uniquely determined by the boundary conditions (7.1) and (7.2). As a result, the solution Ex, Eφ, Kx, Kφ in (t,x,θ,φ) is the spacetime inversion of the solution E x , E ϕ , K x , K ϕ in (t, x, θ, ϕ) by the mapping (4.10)-(4.13). All properties of solutions discussed in section 5.2 is carried over to solutions in (t,x,θ,φ) by (4.10)-(4.13). In particular, the solution gives asymptotically dS 2 × S 2 geometry asz → ∞ (t → −∞ with fixedx or x → ∞ with fixedt) and constantK 1 , (7.13) are constant. The extended spacetime is illustrated by figure 19(b). α 0 , r 0 are independent of R s and the same as in (5.17), while α 1 depends on R s . When we compare to the dS 2 × S 2 geometry obtained in section 5.2 and attempt to glue this two version of dS 2 × S 2 , we find that, if we identify their spatial slices, they have the same E j a (minus sign in Ex is due tox = −x) when they are obtained from the same R s , but they have opposite K a j (no minus sign in K 1 is again due tox = −x), so classically they cannot be glued. By this reason, we denote this new dS 2 × S 2 geometry by dS 2 × S 2 . The time orientation of dS 2 × S 2 is opposite to dS 2 × S 2 . The discussion in section 8 suggests that the transition from dS 2 × S 2 to dS 2 × S 2 may be due to the quantum tunneling.
The solution Ex, Eφ, Kx, Kφ extends the spacetime geometry beyond the last ray. When fixingt =t 0 , the solution Ex(x −t 0 ), Eφ(x −t 0 ), Kx(x −t 0 ), Kφ(x −t 0 ) extends the internal and external geometries of St 0 from region (I) to the future of the blue diamond. Given the numerical solution, we find atz =x −t 0 = − 2 3 R s the marginal anti-trapped surface where Θ l = 0 and Θ k > 0, whilez > − 2 3 R s is the anti-trapped region with both Θ k , Θ l > 0 (see figure 21). Θ l and Θ k are inward and outward null expansion. If R s = 2M r is fixed and we do not consider the loss of white hole mass,z = − 2 3 R s is a null surface being the white hole event horizon. The white hole event horizon may become an spacelike anti-trapping dynamical horizon if the dynamics further reduces M r .

Evidence of quantum tunneling
The discussion in the last section leaves a question about how dS 2 × S 2 and dS 2 × S 2 can be glued to make the black hole interior transit to the white hole interior. We notice that H Δ (recall (3.12)) has a symmetry of reflection K x → −K x and K ϕ → −K ϕ which precisely is the relation between dS 2 × S 2 and dS 2 × S 2 . It is well-known that the symmetry of the Hamiltonian does not imply the symmetry of the solution of EOMs. Indeed from the perspective of the effective dynamics, the reflection symmetry is broken spontaneously at either solution dS 2 × S 2 or dS 2 × S 2 .
We compute the improved Hamiltonian density at dS 2 × S 2 where we have implemented that ∂ x E x = ∂ 2 x E x = 0. ρ is also the density of Gaussian dust by equations (2.10) and (2.11). When fixing E x = r 2 0 , ρ as a function of K 1 , K 2 is plotted in figure 22, where we find that ρ is a double-well effective potential. The values of constant K 1 , K 2 for dS 2 × S 2 and dS 2 × S 2 are located at two zeros of ρ respectively. They are two points √ The right: the cross section of the left in the space of K 1 , K 2 located respectively in two potential wells. Both dS 2 × S 2 and dS 2 × S 2 are not the ground state of ρ since they are not exactly at the minima of ρ, although they are close to the minima. The minima of ρ corresponds to α 0 → ∞ (by computing ∂ K 1 ρ = ∂ K 2 ρ = 0 and applying equation (4.4)), which is not possible to approach by the effective dynamics. Figure 22 suggests an analog with the double-well potential in quantum mechanics, which is the standard example of demonstrating quantum tunneling from one potential well to the other. The low energy states of the double-well model are linear combinations of wave packets located respectively in two wells.
When we understand dS 2 × S 2 and dS 2 × S 2 as quantum states (or wave packets), they should have relatively large quantum fluctuation for the small remnant black hole at the end of the evaporation. This is suggested by the analysis in [45]. In particular if their K x fluctuation is large, the analogy with the double-well potential suggests that there should be a quantum tunneling effect transiting from dS 2 × S 2 to dS 2 × S 2 . We conjecture that the quantum state at the gray region in figure 19(b) should be a sum over dS 2 × S 2 and dS 2 × S 2 , namely the state |r 0 , α 0 ; α 1 in equation (5.19) should be |r 0 , α 0 ; α 1 = 1 √ 2 |r 0 , α 0 ; α 1 dS 2 ×S 2 + |r 0 , α 0 ; α 1 dS 2 ×,S 2 , (8.2) which is well-posed since both dS 2 × S 2 and dS 2 × S 2 share the same values of α 0 , α 1 , r 0 . It is |r 0 , α 0 ; α 1 that exists as the final state of the black hole and the initial state of the white hole. This state with two terms of opposite orientations looks similar to the largej asymptotics of the spinfoam amplitude (see e.g. [56]). The verification of our conjecture goes beyond the effective theory, and relies on studying more carefully the quantization of H Δ and the dynamics in this regime. We expect the effective dynamics is insufficient to approximate the quantum dynamics at the end of the black hole evaporation. The quantum theory of H Δ is just being developed very recently [45]. Interestingly [45] shows that although the effective dynamics developed above is a good approximation Figure 23. (a) Gluing t → ∞ of dS 2 × S 2 tot → −∞ of dS 2 × S 2 is understood as the quantum tunneling (the S 2 factor is suppressed in this figure); (b) both dS 2 × S 2 and dS 2 × S 2 are analytic continued to the Euclidean S 2 × S 2 . The S 2 × S 2 from dS 2 × S 2 has τ/α 0 ∈ [−π/2, π/2] and ψ ∈ [0, 2π) (black half-circle is the τ -obit), while the S 2 × S 2 from dS 2 × S 2 hasτ/α 0 ∈ [−π/2, π/2] (blue dashed half-circle is theτ -obit). The transition from S 2 × S 2 to S 2 × S 2 is at the south pole where the geometry is smooth.
to the quantum dynamics for massive black holes, it needs to be corrected by quantum effects for small black holes at the end of the evaporation, consistent with our above expectation.
The geometry of the Euclidean metric is S 2 × S 2 whose radii are α 0 and r 0 . The coordinate transformation cos(τ/α 0 ) = sin(χ) with τ/α 0 ∈ [− π 2 , π 2 ] makes ds 2 E = α 2 0 dχ 2 + sin 2 (χ)dψ 2 + r 2 0 dθ 2 + sin 2 (θ)dϕ 2 . In the black hole interior, dS 2 × S 2 and dS 2 × S 2 cannot be glued classically because the future boundary t → ∞ slice of the former and the past boundaryt → −∞ slice of the latter cannot be glued smoothly. However when we analytic continue to the Euclidean signature and denote their analytic continuation to be S 2 × S 2 and S 2 × S 2 , transiting from S 2 × S 2 to S 2 × S 2 is a coordinate transformatioñ τ /α 0 = −τ/α 0 ,ψ = −ψ in the first S 2 . The 'future boundary' of S 2 × S 2 at τ /α 0 = π/2 and 'past boundary' of S 2 × S 2 atτ /α 0 = −π/2 are glued at the south pole where the geometry is smooth (see figure 23). The above argument is based on the effective theory. Rigorously speaking, the quantum tunneling discussed in this section is a conjecture, or a proposal. Notice that dS 2 × S 2 has a complete future infinity. We cannot rigorously rule out the possibility that the black-hole-towhite-hole transition does not exist in this model, and the spacetime is completed by dS 2 × S 2 . But what we claim here is that if there exists the black-hole-to-white-hole transition from the Schwarzschild spacetime in our model, it must happen in the region where the effective dynamics is insufficient to describe the dynamics, thus a more careful quantum analysis would be needed. Such a region does exist at the end of the evaporation in our model.
A careful analysis of the quantum dynamics from H Δ is currently undergoing which is expected to provide more details of the expected quantum tunneling.

Conclusion, discussion, and outlook
This paper proposes a new model of the spherical symmetric quantum black hole in the reduced phase space formulation. We deparametrize gravity by coupling to the Gaussian dust. The effective dynamics of the quantum black hole recovers the semiclassical Schwarzschild geometry at low curvature regime and resolves the black hole singularity with Planckian curvature. Our model predicts that the evolution of the black hole at late time reaches the charged Nariai geometry dS 2 × S 2 with Planckian radii ∼ √ Δ. Our model also suggests the existence of quantum tunneling of the Nariai geometry and a scenario of black-hole-to-white-hole transition.
The analysis of this work opens new windows of developments in three phases of quantum black hole dynamics. These three phases are (1) from black hole to the Nariai limit, (2) near the Nariai limit, and (3) from the Nariai limit to white hole. The future analysis of these three phases needs the upgrade from the present effective dynamics to a quantum operator formulation, which is a research undergoing.
As an advantage of the reduced phase space formulation, a proper quantization of the physical Hamiltonian H Δ generates manifestly unitary time evolution. In the phase (1) from black hole to the Nariai limit, it is interesting to investigate the thermalization predicted from H Δ , for instance, questions like whether we can find local observables that thermalize after the formation of the black hole, and if they relates to the Hawking radiation. A standard formalism of addressing these question is the eigenstate thermalization hypothesis (ETH), whose purpose is to explain how (local) thermal equilibrium can be achieved by quantum evolutions from initially far-from-equilibrium states. A wide variety of many-body systems are shown to satisfy ETH, suggesting thermalization should be a generic feature for interacting quantum system (see e.g. [57]). We expect that the H Δ should lead to thermalization of certain local observables. Moreover, thermalization often combines the quantum chaos and information scrambling [58] which are other perspectives to be investigated. In addition, H Δ may be related to the recent studies in [59] on equilibrated pure states after long-time unitary evolution, in order to understand if the long-time unitary evolution of H Δ can be approximately typical and lead equilibrated pure states as outputs. The entanglement entropy of the final state might give an explain of the replica wormholes and page curves following the line of [59].
In the phase (2), it is important to carry out careful analysis for the expected quantum tunneling from the viewpoint of the quantum H Δ , as mentioned earlier.
The numerics shows that the background geometry has approximately vanishing dust density ρ throughout the evolution, given that the initial condition (5.1) and (5.2) corresponds to the vacuum Schwarzschild spacetime (there exists small numerical error, and ρ is bounded by ∼10 −6 throughout the evolution). But perturbations can make the dust density nonvanishing in principle.
Even though perturbations can turn on the dust density, ρ vanishes asymptotically in dS 2 × S 2 no matter if perturbations are turned on or not. Indeed let us consider the PDEs (4.1)-(4.4) and the solution with constant K 1 , K 2 , E x ≡ r 2 0 (these cannot be changed by perturbations). ∂ t E x = 0 and equation (4.3) leads to (a) cos 2β The option (b) is dropped since it reduces equation (4.2) to 1/E x = 0. The option (a) and equation (4.2) gives As a result, the dust stress-energy tensor always vanishes asymptotically In the dynamics on the reduced phase space, ρ = C Δ /(|E ϕ | |E x |) is the physical energy density since H Δ = dx C Δ is the physical Hamiltonian. In terms of states |r 0 , α 0 ; α 1 ∈ H dS 2 ×S 2 , ρ = 0 is understood as expectation values r 0 , α 0 ; α 1 |ρ|r 0 , α 0 ; α 1 = 0, which means that all states |r 0 , α 0 ; α 1 in H dS 2 ×S 2 are infrared soft modes. In particular they have no backreaction to dS 2 × S 2 . The black hole interior containing infinitely many infrared states are anticipated in existing studies of quantum black holes (see e.g. [18] for a summary). ρ → 0 relates to the exponentially large spatial volume in dS 2 × S 2 as t → ∞ (recall equation (5.16)). But the spatial slice has a very narrow throat, since S 2 area is small in the middle see figures 11 and 18 and the horizon radius R s becomes small at late time. dS 2 × S 2 gives an example of Wheeler's bag of gold, since it has the small horizon area but infinitely large volume (see figure 24).
The diffeomorphisms in x-space leaves dS 2 × S 2 invariant but changes α 1 . In quantum notation, the diffeomorphisms generated by V(N) = dx N(x)C(x) define operators acting on infrared states e iεV(N) |r 0 , α 0 ; α 1 = |r 0 , α 0 ; α 1 , where ε ∈ R is an infinitesimal parameter. α 1 (x) comes from the coordinate transformation x → x = x + εN(x) and dη = e −α 1 −α −1 0 x dx = e −α 1 −α −1 0 x dx . In the reduced phase space formulation, the diffeomorphisms are not gauge redundancy but symmetries of the theory. We find the Hilbert space H dS 2 ×S 2 of infrared states is a representation space of the group of onedimensional diffeomorphisms. The x-space is S 1 in dS 2 × S 2 as t → ∞. Therefore H dS 2 ×S 2 carries a representation of Diff(S 1 ) or equivalently H dS 2 ×S 2 carries a representation of Witt algebra: [L m , L n ] = (m − n)L m+n , L n = −i e inθ ∂ ∂θ , (9.7) or Virasoro algebra if we generally allow nontrivial central extension. L n as infinitesimal diffeomorphisms give infinitely many conserved charges (recall equation (3.19)). The above discussion of the infrared modes the Hilbert space H dS 2 ×S 2 are heuristic, and should be further analyzed in more rigorous manner. It is interesting to describe H dS 2 ×S 2 in terms of representations of Diff(S 1 ), and also understand the states as resulting from the unitary evolution of H Δ . We may also looking for their relation with the equilibrated states. In addition, H dS 2 ×S 2 might be embedded in the language of quantum error correcting code, as the code subspace similar as in [60]. It might also relate to states of baby universes [61]. There are debates about whether the modes in the Wheeler's bag of gold form infinitely dimensional Hilbert space, or these infrared modes inside the black hole are linearly dependent (the recent progresses from the AdS/CFT suggests that the Hilbert space may be actually one-dimensional, see e.g. [62]). The unitary dynamics of H Δ should help to clarify the dimension of H dS 2 ×S 2 of the infrared modes.
Rigorously speaking, the existence of the phase (3) relies on the precise description of the phase (2). We expect that the Nariai geometry dS 2 × S 2 is not stable at the quantum level. Its decay and relation to the white hole requires to be further analyzed in detail. It is also interesting to investigate the dynamics of the infrared modes in dS 2 × S 2 toward the white hole and the asymptotically flat regime. As discussed earlier, we expect this dynamics to be highly chaotic. The evidence is provided below.
The quantity C(t) is often called the out-of-time-order correlator, and considered as diagnosis of chaos in quantum systems [39,63,64]. Expectation values of commutators should relate to Poisson brackets in the effective dynamics [14,65], while the corresponding Poisson bracket gives the dependence of the final perturbation δEx(t) on small changes in the initial perturbation δEx(t 0 ). The exponential grows in equation (9.20) is the expected behavior of C(t) in a quantum chaotic system in early time (before the Ehrenfest timet s −t 0 ∼ 1 λ ln 1 Gh ). Our argument of equation (9.20) follows from the usual argument about relating the out-of-time-order correlator to Lyapunov exponent (see e.g. the introduction of [39]).
The dS temperature (the Hawking temperature at the cosmological horizon) relates to α 0 by T dS = 1 2πα 0 [66]. We relate the Lyapunov exponent to the dS temperature by λ 2πT dS . (9.21) This relation resembles the AdS/CFT black hole butterfly effect where the Lyapunov exponent λ CFT of the boundary CFT relates to the black hole Hawking temperature by λ CFT 2πT bh [38,39]. The exponential growth lasts only in a small period near dS 2 × S 2 . The time evolution in (t,x,θ,φ) is just the spacetime inversion of the evolution in (t, x, θ, ϕ) by (9.10)-(9.13). The evolution of perturbations is the same as in figure 14 but evolving in reverse direction.
The above argument is based on the effective approximation. Then it is interesting to understand the chaos in the full quantum theory instead of the effective theory. The detailed analysis of the black-hole-to-white-hole transition should shed light on the resolution of information paradox, given that our discussion is based on the unitary evolution of H Δ . Our model may also related to the argument proposed in [67] where information paradox may be resolved using microscopic structure of the quantum geometry at the Planck scale.
Another important aspect is the strong quantum dynamical regime in the blue diamond in figure 19. The analysis may be carried by the quantization of H Δ , or may require the full theory of LQG (see [9] for a discussion based on spinfoams). We should apply the full theory of LQG to black holes, preferably using the new path integral formulation similar to the recent works on cosmology [44,68,69]. deduced from (4.9)-(4.13).