Evaporation and information puzzle for 2D nonsingular asymptotically flat black holes

We investigate the thermodynamics and the classical and semiclassical dynamics of two-dimensional ($2\text{D}$), asymptotically flat, nonsingular dilatonic black holes. They are characterized by a de Sitter core, allowing for the smearing of the classical singularity, and by the presence of two horizons with a related extremal configuration. For concreteness, we focus on a $2\text{D}$ version of the Hayward black hole. We find a second order thermodynamic phase transition, separating large unstable black holes from stable configurations close to extremality. We first describe the black-hole evaporation process using a quasistatic approximation and we show that it ends in the extremal configuration in an infinite amount of time. We go beyond the quasistatic approximation by numerically integrating the field equations for $2\text{D}$ dilaton gravity coupled to $N$ massless scalar fields, describing the radiation. We find that the inclusion of large backreaction effects ($N \gg 1$) allows for an end-point extremal configuration after a finite evaporation time. Finally, we evaluate the entanglement entropy (EE) of the radiation in the quasistatic approximation and construct the relative Page curve. We find that the EE initially grows, reaches a maximum and then goes down towards zero, in agreement with previous results in the literature. Despite the breakdown of the semiclassical approximation prevents the description of the evaporation process near extremality, we have a clear indication that the end point of the evaporation is a regular, extremal state with vanishing EE of the radiation. This means that the nonunitary evolution, which commonly characterizes the evaporation of singular black holes, could be traced back to the presence of the singularity.


Introduction
Despite the huge recent progress on the observational side, achieved through gravitational-wave detection [1,2] and imaging [3,4], black holes are still a source of challenges for theoretical fundamental physics.The usual black-hole solutions of general relativity (GR) harbour, shielded behind event horizons, spacetime singularities, i.e. regions where the classical and semiclassical descriptions break down [5,6].As discovered more than 50 years ago, black holes behave as thermodynamic systems, whose microscopic description remains, however, still mysterious [7][8][9][10][11][12][13][14].They emit thermal radiation, but the description of the information flow during the evaporation has led to the information paradox, which most embodies the apparent incompatibility between quantum mechanics and GR [15][16][17].
A possible solution to the information puzzle, which has been pursued in the literature, is linking it to the singularity problem [18][19][20][21][22][23].The presence of a spacetime singularity makes the very notion of a global quantum state for matter fields in the black-hole background ill-defined.The loss of unitarity in the evolution of quantum states could be, therefore, traced back to the bad definition of the latter.The main objection to this argument is that the "unitarity crisis" shows up also for macroscopic black holes, i.e. those with masses hierachically larger than the Planck mass m p 10 19 GeV, that is at energy scales where the singularity cannot play any role.Also the possibility to shift the solution of the problem to the late stages of the evaporation, for instance through the formation of Planck-scale remnants, seems untenable owing to difficulty of storing/recovering the huge amount of information at these small scales [24,25].
The recent reformulation and proposal for a solution of the information puzzle [26][27][28][29] seem to bring further arguments against a close relationship between the singularity and the information problem.This is because this novel approach is focused only on reconstructing the correlations between early and late Hawking radiation and, thus, refers mainly to near-horizon physics.
There is, however, an important feature of black-hole solutions, which could change drastically the debate about the relationship between the singularity problem and the information paradox.The most commonly used spacetime setup is that of a black-hole solution with a single event horizon.Black holes with two (an inner and an outer) horizons introduce a new ingredient, which drastically changes the rules of the game.First of all, these black holes typically admit a ground state (GS) represented by an extremal configuration, in which the inner and outer horizons merge in a single one.Moreover, the radius of the extremal black hole could be hierachically larger than the Planck scale [23,30].In the near-horizon region and in the extremal regime, the geometry factorizes as a two-dimensional (2D) anti de Sitter (AdS 2 ) spacetime times a 2D sphere of constant radius.This opens the way to the intriguing possibility that the information issue could be solved in the final stages of the evaporation process using properties of AdS 2 quantum gravity, e.g., by reconstructing correlations between the two disconnected parts of AdS 2 spacetime [31] or by the topological properties of the fragmented GS [32].Moreover, there is some evidence that, for black holes with two-horizons and without a central singularity [33][34][35], the evaporation process could be unitary.Specifically, the presence of the inner Cauchy horizon could act as a trapping region for high energy modes, which could be responsible for the release of information at late times, when the two horizons are about to merge [36,37].
The most natural candidates for testing these ideas are four-dimensional (4D) nonsingular black holes with a de Sitter (dS) core [23,33,[38][39][40][41][42][43][44].They appear as static solutions of Einstein's equations sourced by an anisotropic fluid.The corresponding spacetime is asymptotically flat (AF) and at great distances is indistinguishable from the Schwarzschild solution, whereas the singularity at the origin of the radial coordinate is regularized due to inner dS behavior.The latter also produces an additional hair , which could have interesting observational signatures in the geodesics motion of massless and massive particles, quasinormal modes spectrum and gravitational waves (see, e.g., Refs.[23,30,[45][46][47][48] and references therein).Another consequence of the dS core is the presence, depending on the value of , of two horizons and an extremal solution.From the thermodynamic point of view, these models are characterized by a second order phase transition: the spectrum has a branch of large unstable configurations and a stable branch of near-extremal solutions.
There are two main obstructions that prevent the direct use of such 4D models to address the information paradox.Firstly, we do not have a microscopic model describing the sources of the solutions.We can just give a coarse-grained description in terms of an anisotropic fluid, with equation of state p = −ρ, and a given profile for the energy density ρ.Secondly, there is the usual difficulty of describing semiclassical dynamics, including backreaction effects on the geometry, of quantum Hawking radiation in the 4D classical black-hole background.
In this paper, we show that both issues can be addressed by considering 2D dilaton gravity models of AF, nonsingular black holes with a dS core.As we shall show, these models can be formulated at Lagrangian level and describe, in a simplified setting, the S-wave sector (radial modes) of their 4D nonsingular counterparts.This will allow us to retain the qualitative features of the higher-dimensional models, keeping, however, under control their dynamics.For concreteness, our investigations will be focused on a particular, but quite relevant, case, namely the 2D Hayward black hole.
We will be able to capture the main thermodynamic features of 4D regular models and, at the same time, to describe their evaporation process and to solve the classical and semiclassical dynamics, including the backreaction of Hawking radiation on the geometry.Having under control the latter will allow us to partially answer some important questions regarding the end point of the evaporation process, the time evolution of the entanglement entropy and the shape of the related Page curve [16,49].The main limitation of this approach is obviously represented by the limit of validity of the semiclassical approximation.Nonetheless, our results, together with some known features of AdS 2 quantum gravity, will allow us to have clear indications about the fate of information during the evaporation of nonsingular black holes with a dS core.
The structure of this paper is as follows.
In Section 2, we review some general properties of 2D dilaton gravity models and we present our class of 2D, nonsingular, AF solutions with a dS core.
In Section 3, we introduce the prototype-model we will use throughout our paper, namely the 2D Hayward black hole.
Section 4 is devoted to the investigation of the thermodynamic properties of our 2D models.
In Section 5, we discuss black-hole evaporation using a quasistatic approximation.
The coupling with conformal matter, in the form of N massless scalars, is introduced in Section 6.We consider, in particular, classical solutions corresponding to a shock wave.
In Section 7, we discuss the evaporation process by quantizing matter in the classical gravitational background and by including backreaction effects.The semiclassical dynamics cannot be solved analytically, so we resort to numerical integration.
The entanglement entropy of the Hawking radiation is computed in Section 8 and its Page curve is presented.
Finally, we draw our conclusions in Section 9. Some details of the calculations concerning the absence of divergences for the stress-energy tensor in the GS and the boundary conditions used for the numerical integration of the field equations are presented in Appendix A and Appendix B, respectively.

Two-dimensional regular dilatonic black holes
The simplest black-hole models can be constructed in a 2D spacetime.However, the pure 2D Einstein-Hilbert action is a topological invariant and a metric theory of gravity has to be built by coupling the Ricci scalar with a scalar field φ, the dilaton.2D dilaton gravity is generally described by the action (see Ref. [50] for a review; for a generalization, see Ref. [51]) where D, K and V are functions of the dilaton, representing, respectively, the coupling with the Ricci scalar R (the dimensionless inverse Newton constant), the kinetic term of the scalar field and the potential.Using a Weyl transformation of the metric together with a field redefinition it is always possible to set K = 0 and to recast the lagrangian into the simpler form [52,53] This choice of the conformal (Weyl) frame is typically used when dealing with asymptotically AdS black holes, while for AF configurations (which is the focus of the present paper), a conformal frame with K = 0 is generally considered more appropriate [54][55][56][57][58].This is particularly true when 2D dilatonic black holes are used to describe the S-wave sector of higher dimensional models (but see, e.g., Refs.[59][60][61][62][63]).However, as we shall see below, the description of also AF black holes is much simpler in the conformal frame (4) than in a frame with K = 0.Moreover, the lagrangian ( 4) is fully characterized by the dilaton potential V (φ); this allows for a simple classification of regular black-hole models in terms of the properties of V .For these reasons, in the following, we will work in this frame, although this choice will introduce some difficulties concerning the physical interpretation of the parameters we use to describe black holes.The equations of motion stemming from Eq. ( 4), in the absence of matter fields, read

Linear dilaton solution
Let us now first consider static solutions of Eq. (5b).In this case, the dilaton can be used as a spacelike "radial" coordinate of the 2D spacetime, φ = λr, where λ is a constant, with dimensions of the inverse of a length, characterizing the potential V .This parametrization of the dilaton is particularly useful when the 2D theory is used to describe the S-wave sector of 4D black holes.In this case, the dilaton is proportional to the radius of the transverse two-sphere.Notice also that the dilaton represents the inverse of the 2D Newton constant.This means that the region r λ −1 is in a strong coupling regime, whereas r λ −1 is a weak-coupling region.To be consistent, both interpretations of φ require to limit the range of variation of the radial coordinate to r ∈ [0, ∞).
Eq. ( 6) allows us to write the most general static solution of Eqs.(5a) and (5b) as the linear dilaton solution (LDS) where c 1 is a dimensionless integration constant, which can be written in terms of the covariant mass M, which can be defined for a generic 2D dilaton gravity theory [64].For a static spacetime, M is the conserved charge associated with the Killing vector χ µ = F 0 µν ∂ ν φ, generating time translations 1 .F 0 is a constant, which is fixed by the normalization of χ µ .As we will see in the following, to make contact with four dimensional models, a convenient choice2 is F 0 = −1/λ.
In our Weyl frame, M reads Choosing appropriately the form of V allows to generate different solutions.In particular, we focus on dilaton gravity models allowing for AF nonsingular black holes.
A useful information for classifying different classes of models can be obtained from the existence of solutions characterized by a constant dilaton, the so-called constant dilaton vacua (CDV).Owing to the r-dependent parametrization of the dilaton (6), these solutions cannot be obtained as particular LDSs given by Eq. ( 7) and must be discussed separately.
Notice that the condition for having a flat CDV is rather strong.It requires φ 0 to be both a zero and an extremum of V .
Asymptotic flatness, however, also implies that we always have V (∞) = dV /dφ| ∞ = 0.In this case, formally, we can consider φ = ∞, which corresponds to a decoupled configuration (the 2D Newton constant vanishes), as a flat CDV.
In the first case above, instead, we can define From this, using Eq.(5a), we have which describes two-dimensional dS (dS 2 ) spacetime, with an associated dS length L dS .
In the second case, we define instead Equation (5a) now yields which describes the AdS 2 spacetime, with an associated AdS length L AdS .

General class of 2D nonsingular, asymptotically-flat black holes with a de Sitter core
We are interested in 2D dilatonic black holes, which mimic the behavior of 4D regular black holes.Therefore we have to choose the form of the potential V such that: (a) the spacetime has one or at most two horizons; (b) the spacetime curvature remains everywhere finite, in particular the usual 4D curvature singularity at r = 0 is removed; (c) the spacetime is AF, with a Schwarzschild subleading behavior at r → ∞, i.e., f ∼ 1 − c/r, with c a constant.
Condition (a) requires f (r) to have at least one zero and to stay positive outside the (outer) horizon, i.e., Condition (b) can be reformulated as a condition on the first derivative of V .In two dimensions, the Ricci scalar is the only curvature invariant.Therefore, bounding it is sufficient to generate regular spacetime metrics.According to Eq. (5a), this translates into requiring regularity of dV /dφ.Regularity at r = 0 can be achieved in different ways.The simplest and more physical one, which has also been used for 4D models [23,33,[38][39][40][41][42][43], is to impose dV /dφ| 0 to be finite and According to the general discussion of Section 2.2, this implies that our model must allow for a dS 2 CDV at φ = 0, given by Eq. ( 10), which will therefore describe the inner core of our black-hole solutions.Using Eq. ( 10), one can easily find the form of the potential in the φ ∼ 0 region is, Condition (c), i.e., asymptotic flatness and a Schwarzschild subleading behavior, can be implemented by fixing the asymptotic behavior for φ → ∞ Figure 1: Qualitative behavior of the potential, characterizing the broad class of the regular models satisfying the conditions (1), (2), (3).We restricted ourselves to the case of a single maximum and a single minimum.The horizontal and vertical lines correspond to V = 0 and φ = 0, respectively.Equation ( 14) implies the dilaton potential to be zero at φ = 0, and to become negative and decrease near φ ∼ 0. However, it has to grow again, cross the φ-axis at a finite value φ = φ 1 , develop at least one minimum and one maximum to guarantee the positive fall of V at asymptotic infinity, implied by Eq. ( 16).The model, therefore, must allow for three different CDV solutions at φ = 0, φ = φ 1 and φ = ∞, describing, respectively, a dS 2 , AdS 2 and flat spacetimes.
In principle, the potential could show any number of oscillations, but for simplicity in the following we restrict ourselves to potentials with a single maximum and a single minimum.The qualitative behavior of the general form of our potential V is shown in Fig. 1.Quantitatively, the potential will depend on some dimensional parameters.The most natural, minimal choice is a potential depending on two parameters.Since φ is dimensionless, the parameter λ, introduced above, is needed to give the right dimensions to V .At least a second parameter, however, which we will call , is needed if we want to express the CDV φ 0 in terms of parameters of the model.The values of such parameters will impact on the behavior of the metric function.Indeed, since V = λ df /dr from Eq. ( 7), we see that the presence of two zeros for V implies the existence of a minimum for f (r).Depending on the value of , this minimum can be above, below or exactly at the r-axis, producing horizonless, two-horizon or extremal configurations, respectively.The parameter λ, instead, does not affect the presence of extrema in f (r), nor their location in the radial direction.Notice that the AdS 2 CDV describes the near-horizon behavior of the extremal black hole [67,68].
One can easily construct dilatonic potentials V (φ) behaving as in Fig. 1.Basically, for every sphericallysymmetric, regular 4D black hole, characterized by a single metric function f , one can easily construct the corresponding 2D dilaton gravity theory by solving Eq. (7), determining in this way the form of V .For instance, notable models are those which can be obtained from the Hayward black hole [33], Gaussian-core black hole [40], the Fan-Wang model [43] or the Bardeen solution [38].
For the sake of concreteness, in the following we will focus on a 2D dilaton gravity model reproducing the Hayward black hole.However, all the considerations of the next sections can be extended to the general class of models described in this section.

Two-dimensional Hayward black hole
One of the simplest cases of potentials behaving as shown in Fig. 1 is given by where is a parameter with dimensions of length.The potential has a zero at φ = 0, which gives the dS CDV with a related dS length (9), given in this case by L 2 dS = λ 3 , and goes to zero for φ → ∞.The other zero is at φ λ which gives the AdS 2 CDV and, as we shall see below, describes extremal black holes in the near-horizon regime.The associated AdS length (see Eq. ( 11)) is L 2 AdS = 3λ 3 .With the potential (17), solving Eq. ( 7) yields which interpolates between the Schwarzschild spacetime at great distances and the dS one at r ∼ 0, modulo a rescaling of the coordinates t and r by the constant quantity 2M/λ.This peculiar behavior, in which the mass term in the line element dominates at great distances, was analyzed in 2D very recently in Ref. [69] and termed "mass-dominated" dilaton gravity.The 4D Hayward black hole [33] is described by the metric element ds 2 4 = −f H (r) dt 2 +f −1 H (r) dr 2 +r 2 dΩ 2 , with the metric function given by where m is the 4D ADM mass.One can easily check that the (constant) Weyl rescaling of the 2D metric, together with a rescaling of the time coordinate brings the metric into the form This transformation leaves the 2D dilaton gravity action invariant up to a constant factor, which does not alter the equations of motion.
If we now write the covariant mass of the 2D solution in terms the mass m of the 4D black-hole solution (20) we get The 4D metric element of the Hayward black hole can be simply written in terms of the 2D one (22) and the dilaton as ds 2 4 = ds 2 2 + (φ/λ) 2 dΩ 2 .The peculiar relation (23) can be seen as a consequence of both the specific conformal frame chosen here (see the discussion at the beginning of Section 2), which is particularly suited for asymptotically AdS spacetimes, and of the normalization of the Killing vector F 0 adopted in Section 2.1.As we shall see, the minus sign in F 0 implies that an asymptotic observer in 2D spacetime measures the energy of the system with the opposite sign with respect to the asymptotic observer in the 4D spacetime.Hence, when the 2D mass becomes bigger, the corresponding 4D mass decreases and viceversa, which is reflected in the inverse relation (23).We will further confirm this below, when studying the thermodynamic properties of the 2D model.
In the remainder of the paper, we will consider the metric function f in the form (19).It has a minimum at If f (r min ) < 0, the metric has two horizons, solutions of f (r) = 0, while if f (r min ) > 0, it has no horizons.If f (r min ) = 0, the two horizons merge, become degenerate and the configuration becomes extremal, with an event horizon located at r min ≡ r ext .Using the latter and setting f (r ext ) = 0 yields the critical value of at extremality Thus, for < c the black hole has two horizons; for = c the two horizon merge in a single one; whereas for > c the spacetime has no horizons.Interestingly, the value of r ext in Eq. ( 24) is the same at which the potential (17) changes sign (see Eq. ( 18)).Indeed, as it is usually the case for two-horizon models [23,70,71], the extremal, near-horizon metric is that of AdS 2 spacetime.Our 2D black-hole solutions can be considered as thermodynamic systems, characterized by a Hawking tempertaure T H , an internal energy E and an entropy S. The Hawking temperature is given by the standard formula where r H is the radius of the (outer) event horizon.The temperature becomes zero both at extremality, i.e., for r H = r ext =3 √ 2 , and for r H → ∞, whereas it reaches a maximum at The qualitative behavior of T H is plotted in Fig. 2. The internal energy is usually identified with the black-hole mass.However, in our case, we have chosen a negative normalization of the Killing vector generating time translation, opposite to the usual positive one.Consistency with this normalization requires E = −M.Using Eq. ( 19), we can express M as a function of the outer event horizon radius A negative internal energy is somehow unusual for black holes, but it is perfectly consistent with their description (and normalizations) as thermodynamic systems.Indeed, we will confirm this below by using the euclidean action approach and proving the consistency of the first law of thermodynamics.The energy E, thus, is always negative and goes to its maximum value E = 0 as r H → ∞, whereas it reaches its minimum value for the extremal black hole, whose radius is given by Eq. ( 24) 3 .Moreover, for r H ≥ r ext , as expected, E(r H ) is a monotonic increasing function.
Let us now calculate the entropy of the 2D Hayward black hole using the Euclidean action formalism, which allows us to calculate the partition function Z of the thermodynamic ensemble in terms of the Euclidean action I, i.e., I = − ln Z.By a Wick rotation of the time t, the action of the lagrangian (4) becomes the euclidean bulk action I b .It has been shown that, in order to have a well-defined variational principle for "mass-dominated" dilaton gravity theories, the action must be supported by (one half) the usual Gibbons-Hawking-York (GHY) boundary term I GHY plus an additional one, containing the normal derivative of the dilaton [69].
The action reads We enclose the system into a hypersurface at constant r = r ∞ , where we define an induced, one-dimensional, metric h, whose extrinsic curvature is described by K µν (K is its trace).τ is the euclidean time, which is periodic with period equal to the inverse of the temperature T −1 H ≡ β.All quantities on the boundary will be evaluated at the cutoff r = r ∞ and then we will let r ∞ go to infinity.
Moreover, one could add a purely topological, Einstein-Hilbert term of the action I topo ∝ φ 0 d 2 x √ −gR, with φ 0 a constant, which only changes the value of the entropy by an additional constant S 0 , depending on φ 0 .This constant value can be identified as the entropy of the extremal configuration.
Let us now evaluate the boundary term on the LDS.The induced metric with euclidean signature reads h µν = h 00 = f .The extrinsic curvature is defined in terms of the normal vector to the hypersurface n µ as where the normal vector reads, in this case, n µ = f −1/2 δ r µ .Therefore, when evaluated on the solution for the dilaton, we have I GHY = − λ 4 β r f | r=r∞ , which vanishes in the limit r ∞ → ∞ when f is given by Eq. (19).The remaining boundary term, instead, gives Usually, one has to add a counterterm to the boundary action, needed to regularize divergences arising in the limit r ∞ → ∞.This counterterm is written in terms of the extrinsic curvature of the boundary embedded in flat spacetime.In our case, such term is not needed because there are no divergences.Moreover, there is no contribution from flat spacetime (f (r) = 1) since Eq. ( 29) gives S flat = 2 dt dr φ , which is zero for the LDS (7).
We now evaluate the bulk action.Using Eq. ( 7), we have where we used f (r Combining Eqs. ( 30) and ( 31) yields I = −βM − 2πλr H = − ln Z, where Z is the partition function.The internal energy and entropy, thus, read Equation (32a) confirms Eq. ( 28), as expected.
The black-hole entropy, instead, scales as the dilaton, i.e., as the inverse 2D Newton constant, evaluated at the horizon.This is the usual formula for the entropy of 2D black holes [72,73] and represents the extension to two spacetime dimensions of the usual area law in higher dimensions.This is also quite evident when the 2D black hole is derived from the dimensional reduction of the 4D one, with the dilaton playing the role of the radius of the transverse S 2 sphere, as it is here the case.
Contrary to standard 4D regular black-hole solutions [23], here the entropy naturally follows the area law.This is because our 2D solutions do not require external matter sources and, therefore, are not coupled to a stress-energy tensor, which in the 4D case describes an anistropic fluid.In general, for regular models, the latter is characterized by a density that depends on the ADM mass of the model, which introduces extra (bulk) terms altering the area-scaling of the entropy [23,74].
Adding the contribution S 0 of the topological action leads to Because of the minus sign in Eq. ( 28), we need to check the consistency of our derivation with the first principle of thermodynamics.By differentiating E and S, given respectively by Eqs. ( 28) and (32b), with respect to r H , and using Eq. ( 26) for T H , we can easily check that the identity dE = T H dS is satisfied.
The previous results allow us to compute the energy difference between two configurations, characterized by the two values φ 1 and φ 2 of the dilaton, in terms of the integral of the dilaton potential.This can be done in all generality by integrating the first law, considering Eqs. ( 26) and (32b) Let us assume φ 2 < φ 1 and that φ 2 represents a black-hole configuration.Then, for φ 2 , the condition (13) holds, whether or not φ 1 represents a black hole.Therefore, ∆E 1,2 > 0 and the black-hole energy increases monotonically, as already seen by analysing Eq. ( 28).The consequence of this is that the configuration retaining the least internal energy will be the one with the least dilaton, i.e., the extremal configuration.

Thermodynamic stability and second order phase transition
Let us now investigate the thermodynamic stability of our regular black-hole solutions.This can be done by studying the specific heat and the free energy.

Specific Heat
The specific heat of our solutions is given by In Section 4.1, we showed that the internal energy of black-hole configurations is always increasing with φ H . Therefore, the sign of C H is determined by the sign of dV /dφ H .As already discussed in Section 2, requiring the potential to satisfy the minimal requirements listed in Section 2.3 implies V to be necessarily nonmonotonic.Moreover, imposing a dS-like behavior in the interior constrains the potential to have another zero at φ H = φ 0 , which further restricts the interval to φ H ∈ [φ 0 , ∞) where we have at least an extremum.The specific heat, thus, shows a single maximum in this interval, located at φ peak (see Fig. 1).Here, dV /dφ H changes sign (from positive to negative) as the potential falls monotonically at infinity.Therefore, we have This correspods to a branch of thermodynamically stable configurations; • For φ H > φ peak , dV /dφ H < 0 and therefore C H < 0. This is, instead, the branch of thermodynamically unstable configurations; • For φ H = φ peak , dV /dφ H = 0 and C H → ∞, which signals the onset of a second order phase transition.This general discussion also applies to the 2D Hayward black hole, described by the potential (17).Using Eqs. ( 26) and ( 28), Eq. (35) reads which diverges at the peak temperature (27).In terms of the event horizon radius we have • An unstable branch of large black holes (r H r ext ), with negative specific heat; • A stable branch of black holes close to extremality (r H r ext ) with positive specific heat.Here, we highlighted the two configurations corresponding to φ 1 and φ 2 , respectively in the stable and unstable branch and their common temperature T = T (φ 1 ) = T (φ 2 ).The areas representing (φ 2 − φ 1 )T and the right-hand side of Eq. ( 38), are coloured in blue and red, respectively.
C H goes to zero at extremality and, near r H ∼ r ext , it behaves as which scales linearly with r H , as the specific heat of AdS 2 black holes [75].This is fully consistent with the AdS 2 behavior of the 4D LDS Hayward metric in the near-extremal, near-horizon regime.However, differently from its four dimensional counterpart, for which the AdS 2 spacetime represents only an approximate solution, the AdS 2 appears also as an exact CDV solution for the 2D model.

Free energy
The analysis of the specific heat allows us to distinguish between stable and unstable branches, but it is not sufficient to select the energetically preferred configurations.To this end, we need to analyze the difference in the free energy between different configurations sharing the same temperature.We will consider LDS on different thermodynamic branches, with the same temperature, but different dilaton (and, therefore, different horizon radius).In our analysis, we will not consider the CDV solutions, but only LDS.The inclusion of the CDV will be discussed in the next subsection.
We consider the situation depicted in Fig. 4, where we present a qualitative plot of the temperature T (φ) of our 2D black holes.We focus on two different configurations, φ 1 and φ 2 , in the stable and unstable branches, respectively, but with the same temperature T (φ 1 ) = T (φ 2 ) ≡ T .We now evaluate the free energy difference ∆F 2,1 ≡ F (φ 2 ) − F (φ 1 ) = ∆E 2,1 − T ∆S 2,1 .Using Eqs.(32b) and (34), we have From a geometric point of view, ∆F 2,1 /2π represents the area limited by the graph of T (φ) (solid blue line in Fig. 4) and the straight line T (dashed horizontal line in Fig. 4).Since T (φ) shows a single maximum located at φ = φ peak and φ 1 < φ peak < φ 2 , then the relation This implies ∆F 2,1 to be strictly positive, i.e., F (φ 2 ) > F (φ 1 ).Therefore, generic configurations in the stable branch are energetically preferred with respect to configurations in the stable branch.These results confirm those obtained in Section 4.2.1 and show that configurations in the stable branch retain the least free energy and are, thus, thermodynamically favoured.

Including the constant dilaton vacuum
So far, in our thermodynamic considerations, we have considered only the LDSs.We have already seen in Section 2, however, that our dilaton gravity model allows, in its spectrum, a solution with constant dilaton describing an AdS 2 spacetime, the CDV.This solution represents a GS of the theory, with different asymptotics not only with respect to the "excited" LDS, but also to the extremal one.In fact, the LDS, including the extremal one, are metrically AF and the dilaton depends linearly on r.Conversely, the CDV describes an AdS 2 spacetime and the dilaton is identically constant.Extremal LDSs, in the near-horizon approximation, are described by an AdS 2 spacetime, endowed, however, with a linear dilaton, the so called linear dilaton vacuum (LDV) [68].
The situation described above is very common for 4D charged black holes [76], for which we have both an extremal, AF solution, described in the near-horizon approximation by AdS 2 with a linear dilaton, and an AdS 2 × S 2 solution (our CDV).
Working in the context of 2D dilaton gravity, it has been shown that the AdS 2 CDV does not admit finite energy excitations, i.e., it is separated from the AdS 2 LDV by a mass gap [77].Moreover, there is the additional difficulty that the two spacetimes have different asymptotics (linearly varying versus constant dilaton).The latter point makes it conceptually problematic to compare the free energies of the two configurations and to assess which one is thermodynamically favoured.These difficulties can be circumvented, and one can show, computing the free energy, that the AdS 2 CDV is energetically preferred with respect to the AdS 2 LDV [68].Here, we will use a similar procedure for the solutions of the 2D Hayward model and compare the free energy of the AdS 2 CDV with that of the extremal LDS.
Let us first note that we can formally consider zero mass, thermal excitations of the CDV using a Rindler-like coordinate transformation, which generates a horizon with a related temperature We can now evaluate the bulk euclidean action of the CDV (39), considering that the dilaton is constant φ = φ CDV and the potential is zero when evaluated in the CDV.We have4 The free energy F CDV = −T ln Z reads We can now compute the difference in the free energy between a generic black-hole configuration of the LDS and the CDV.Using the result of Eq. (32b) for the entropy of the LDS and the equation Since, for every black-hole solution, φ H > φ CDV , we have F BH < F CDV , and thus the black-hole configuration is always thermodynamically preferred.However, this is true whenever T H = 0.At T H = 0, we do not have thermal contributions to the free energy anymore and ∆F reduces to the difference of the masses contributions.On the other hand, at T H = 0 the semiclassical approximation is broken, as signalized by the generation of the mass gap (see Section 5.3) [32,68,77], and thus we cannot rely on the euclidean action approach anymore.This does not allow to define a proper mass for the CDV, consistently with the fact that pure AdS 2 spacetime does not admit finite energy excitations.One could argue, following the argument of Ref. [68], that at T H = 0 the only contribution to ∆F BH−CDV comes from the mass difference, which, due the absence of finite energy excitations, should diverge, which makes the CDV energetically preferred with respect to the extremal LDV.

Black-hole evaporation in the quasistatic approximation
In the following, we will describe the evaporation process of our regular 2D black hole working in the quasistatic approximation and in the semiclassical regime, in which the mass is slowly varying with time so that it can be considered almost constant for each individual evaporation step.In this way, the backreaction of the geometry due to the radiation is not taken into account in a fully dynamic way, but it is described in a very simplified, rough manner.The dynamic character of this backreaction will be, instead, fully taken into account in the next section, where we will consider the coupling of gravity to the matter fields describing Hawking radiation.We expect our quasistatic approximation to hold for black holes very far from extremality and to break down in the near-extremal regime, where the semiclassical approximation is not capable of describing the dynamics.
Since our black holes behave as black bodies with a Planckian thermal spectrum, we use the Stefan-Boltzmann (SB) law to describe the time variation of the internal energy, which in arbitrary d + 1 dimensions reads [78] dE dt where σ is the SB constant, A d−1 is the (d − 1)-dimensional emitting surface, Γ(x) and ζ(s) are the gamma and Riemann zeta functions, respectively.In the present case, d = 1, and thus where we used the fact that E = −M.

Evaporation time
To compute the evaporation time, we express M as a function of the event horizon radius and use dM/dt = (dM/dr H )(dr H /dt) together with Eqs. ( 26), ( 28) and (44), to obtain Inverting and integrating yields the evaporation time required to pass from an initial configuration with event horizon radius r H,0 r ext to a final one with radius r H, final From this expression, it is already evident that, as T H → 0, i.e., as we approach extremality, ∆t → ∞.In other words, reaching the extremal configuration, in the quasistatic semiclassical approximation, requires an infinite time, consistently with the thermodynamic stability analysis of the previous subsections.
For the particular model under investigation described by Eq. ( 19), Eq. ( 46) reads In the extremal limit r H, final → 3 √ 2 = r ext we have a logarithmic divergence, as expected.This is consistent with the behavior of 4D regular models (see [79][80][81] and references therein).

Time variation of mass and entropy
Let us now compute how the mass of our 2D Hayward black hole evolves in time due to the emission of Hawking radiation, according to the SB law (44).This can be done by first inverting r H and T H , to express   44) (blue solid line), which shows the time evolution of the mass of an evaporating black hole as a function of time, in the quasistatic approximation.We see that at large times, the numerical solution asymptotes the extremal mass (orange dashed line).For the numerical integration, we set M(t = 0) = 0.1 λ as an initial condition.Figure (b): Time variation of the entropy of the black hole according to Eq. (33).At late times, the entropy goes to zero, after subtracting the topological term S 0 .In both figures, we set λ = = 1.
them as functions of M and then by numerically solving Eq. ( 44) written in the form dM dt = π 6 T 2 H (M) .As a boundary condition, we impose M(t = 0) = 0.1 λ.We also set λ = = 1.Using the expression for r H as a function of the mass obtained by inverting Eq. ( 28), for M = 0.1, we have r H 4.96, which is about twice as large as that pertaining to the temperature peak (27) (r H, peak 1.9 with = 1), which confirms that the initial state belongs to the unstable branch.We have also checked that the final results are independent of the value of the initial mass.
The result of the numerical integration is shown in Fig. 5(a).We see that, at large times, the mass asymptotes the extremal one.The mass increases during the evaporation due to the negative normalization of the Killing vector generating time translations.Indeed, the internal energy E = −M is decreasing, as it should be during the evaporation.
We can also derive the time evolution of the black-hole entropy.This can be done by simply combining the solution of Eq. ( 44) together with the function r H (M) (obtained by inverting Eq. ( 28)) and Eq. ( 33).The result is reported in Fig. 5(b) and confirms the fact that, as the solution asymptotes the extremal one, the entropy reduces to zero.This is true only if we subtract the contribution S 0 of Eq. (33).
These results will be extended in Section 7, where we will go beyond the quasistatic approximation and we will consider the full dynamics of the backreaction of Hawking radiation on the background geometry.We will show that the inclusion of the latter causes the evaporation process to take place in a finite time.

Approaching extremality and breakdown of the semiclassical approximation
The divergence of the evaporation time for black holes approaching extremality, found in the previous section, signalizes the breakdown of the semiclassical approximation.Let us now study in details how an excited configuration approaches the extremal limit, by solving Eq. ( 45) at leading order around extremality, i.e., around r H ∼ 3 √ 2 .Near extremality, at leading order, the temperature varies linearly with r H , which gives, after solving Eq. ( 45), the time-dependence of r H r H (t) where α 1 is an integration constant.Moreover, M(t) behaves, near extremality, as which reduces to 1 3 3 √ 2 = M ext only for t → ∞, confirming the numerical results of the previous subsection.The entropy, instead, approaches exponentially that of the extremal configuration at t → ∞, according to Eq. ( 49).
However, one must question the validity of the semiclassical approximation near extremality.The latter breaks down when the energy of Hawking quanta, which is of order T H , becomes comparable with the energy of the black-hole energy above extremality ∆E = |M − M ext | (see, e.g., Ref. [32]).The energy scale at which this breakdown occurs determines the mass gap separating the CDV from the continuous part of the spectrum (the LDS) [32,77].
This mass gap can be determined by expanding Eq. ( 28) near extremality, which yields From ∆E T H , we easily find the energy gap This result is consistent with those obtained for 4D black holes with two horizons merging into a single one, in particular for charged black holes.For instance, in the Reissner-Nordström case, the energy gap behaves as E gap ∝ Q −3 , where Q is the black-hole charge [32].As one can expect, here the role of Q is played by 5 .Figure 6 It is interesting to notice that the limit of validity of the semiclassical approximation sets also the limit of validity of the quasistatic one, which cannot be valid when the former is broken.In fact, the quasistatic approximation is valid in the initial stages, when the black hole is macroscopic, which essentially evaporates as a GR one.It remains also valid for most of the evaporation process because the evaporation time is much larger than the typical black-hole time scale 1/∆E (see also, e.g., Refs.[81,82]).Only at extremality, when ∆E goes to zero, the two time scales become comparable.
Moreover, it should be considered that the time at which the semiclassical approximation breaks down (which can be read from Fig. 6(a)) is close to the time at which the evaporation process reaches the maximum of the temperature (27), where we have the onset of the second order phase transition (see Fig. 6(b)).

Coupling to conformal matter
In the previous section we have described the black-hole evaporation process in the semiclassical and quasistatic approximations.Within these approximations, the backreaction effects of the geometry on the presence of Hawking radiation is completely encoded in the change in time of the black-hole mass M. As we have seen above, this may be a good approximation in the early stages of the evaporation, but it is   52).
As it can be seen, the breakdown of the semiclassical approximation happens after the evaporation process reaches the temperature peak.
expected to fail at later times.Another shortcoming of the approximation is that the backreaction is not fully dynamic, since it does not involve the full metric solution.Its role is simply encoded in the variation of the black-hole mass.
In this section, we will give an exact semiclassical description of the evaporation process by studying the coupling of our 2D model to quantum conformal matter, in the form of N massless scalar fields.This coupling is most easily analyzed in the conformal gauge, where the 2D metric reads where e 2ρ is the conformal factor of the metric.The transition from the metric in the Schwarzschild gauge (7) to that in the form Eq. ( 53) is realized by using the coordinates, x ± = t ± r * , where r * ≡ dr/f is the tortoise coordinate.The system of coordinates x ± does not cover the interior of the black hole, but only the region outside the outer horizon.Indeed, x + − x − → −∞ corresponds to the horizon (x + → −∞ gives the past horizon, while x − → ∞ future horizon), while x + − x − → ∞ corresponds to asymptotic infinity.The field equations stemming from Eq. ( 4) are now The solution reads where we defined J ≡ 1 λ 2 φ dψ V (ψ).

Coupling to matter: shock wave solution
We now couple 2D dilaton gravity to the N massless scalar fields describing conformal matter.The full action reads The stress-energy tensor of matter fields reads The field equations ( 54) now become The matter-coupled field equations above admit an exact solution if we consider an ingoing shock wave starting at x + = x + 0 and propagating in the x − direction, while no energy flux is present in the x + direction The minus sign in T ++ is again due to the normalization of the Killing vector of the metric.From Birkhoff's theorem, we can write the full solution by patching, on the infall line x + = x + 0 , the vacuum solution together with the one after the shock wave [62].
At the end of Section 4.1, we showed that, in the linear dilaton case, the GS of the theory, i.e., the state retaining the least internal energy, is the extremal black-hole configuration, characterized by a mass M ext given by Eq. (25).
The vacuum solution (before the shock wave) therefore is equivalent to Eq. ( 55) with x ≥ x + 0 Since T −− = 0, now the solution is (see Ref. [62]) where ) is a function needed to map the old coordinate x − of the observer in the GS solution into a new coordinate, which pertains to an observer in the excited solution.Its form is fixed by requiring continuity of the metric function across the shock wave at x + = x + 0 .This defines the function up to an integration constant, which is fixed by requiring the continuity of the dilaton across the shock wave.
As usual, and as we will see in more details in Section 7.2, the function F (x − ) generates the Hawking flux of particles, which can be described in terms of the change of the coordinate x − defined by F .In fact, in 2D, the flux of Hawking particles can be described in terms of the Schwarzian derivative of the function F (x − ) (see, e.g., Refs.[62,83]).
Due to the sign of the shock wave (59) and to the minus sign in Eq. ( 28), the mass M of the excited states is less than the extremal mass, as the shock wave increases the internal energy of the system.The physical picture we expect is therefore the following (see also Fig. 7).The shock wave increases the internal energy of the initial system, i.e., the GS, extremal configuration: the degenerate horizon splits into two apparent horizons.According to the thermodynamic analysis and the results of Section 5.2, we then expect that, when dynamically evolving the system, these two horizons will meet again at the end of evaporation, and merge to give again the extremal configuration.
Figure 7: Metric functions with varying mass M. The blue solid line refers to the extremal configuration with M = M ext .The orange dashed line, instead, refers to an excited state (after the shock wave) with M = 0.9 M ext , which presents two event horizons.
We will confirm this in the next section by keeping into account the dynamic contributions of backreaction effects of the radiation on the geometry and dynamically describing the evaporation process by numerically integrating the field equations.

Black hole evaporation and backreaction
The evaporation and its backreaction effects on the spacetime geometry, are studied by quantizing the conformal matter on the curved 2D background.An important consequence of the curvature of the spacetime is that the otherwise classically traceless stress-energy tensor acquires a nonzero trace, proportional to the Ricci scalar, which is the so-called conformal anomaly [84] T where N is the number of matter fields.This can be accounted for by adding, to the classical action, the nonlocal Polyakov term where 2 −1 is the scalar Green function.Using Eq. ( 62), we can derive the expectation value of T +− , which in the conformal gauge is entirely local The nonlocal effects stemming from the action ( 63) are, instead, encoded in the other components of the stress-energy tensor T ±± , which can be obtained from covariant conservation of the latter Here, t ± (x ± ) are integration functions, which depend on the boundary conditions and therefore encode the nonlocal effects of the Polyakov action (63).These functions are also sensitive to the choice of the coordinates.Indeed, under a conformal transformation of the coordinates x ± → y ± (x ± ), the conformal factor transforms as which, plugged into Eqs.( 64), (65a) and (65b), yields the anomalous transformation of the stress-energy tensor with the Schwarzian derivative where the prime indicates a derivative with respect to x.The form of the stress-energy tensor in the new coordinate system is preserved if the t ± 's trasform as Including the conformal anomaly, the field equations ( 58) become Eqs. ( 69) can be solved once suitable initial conditions to fix the functions t ± (x ± ) are imposed.These can be determined assuming the GS as the initial state.In conventional AF models, like the CGHS one [54], the GS is pure Minkowski spacetime.One can therefore define a global coordinate transformation in which the conformal metric ( 53) is manifestly flat, i.e., we can define a system of coordinates in which e 2ρ = constant = 1.One can then assume that there is no incoming radiation (except from the classical shock wave) and that there is no net outcoming flux, so that identically, which implies t ± = 0 on the GS in this system of coordinates.One can then transform back to the original coordinates and exploit the anomalous transformation (68) to obtain their final form in the new coordinates.
In the case under consideration, however, we saw that the GS does not correspond to Minkowski spacetime (which is only reached asymptotically), but it is given by the extremal configuration (60).This, of course, prevents from defining a global coordinate transformation which brings e 2ρ → constant.Despite this difficulty, we can still use Eq. ( 70) as a boundary condition, similarly to the CGHS model.
Once the boundary conditions on the GS have been imposed, the solution before the shock wave (x + < x + 0 ) is Eq. ( 60), the vacuum one, while after the shock wave (x + > x + 0 ) it is given by an evaporating black-hole solution.

Adding counterterms and fixing the boundary conditions
In order to preserve the physically motivated boundary condition (70), we can follow Refs.[62,65] and modify the usual Polyakov action ( 63) by adding the most general local covariant counterterms with no second order derivatives where A and B are functions of the scalar field.The presence of these new terms, of course, does not alter the classical limit N → 0. Also notice that the addition of new counterterms was already employed in, e.g., the CGHS model, in order to make the theory exactly solvable [85,86].In Ref. [62], the addition of the counterterms was necessary to prevent T µν GS from diverging for φ → ∞.In the present case, it can be shown that divergences are absent due to the peculiar properties of the potential outlined in Section 2.3 (see Appendix A).Nevertheless, adding counterterms is needed to implement the boundary condition (70) in a consistent way.With the new terms, the components of the stress-energy tensor ( 64), (65a) and (65b) become We now impose the boundary condition (70).Requiring t ± (x ± ) = 0 on the GS completely fixes the two functions A and B (see also Ref. [62]) where the subscript GS indicates that J is computed at extremality.

Equations (72a) and (72b) now read
Since the scalar field is a function of x + and x − , while we are treating J as a function of φ, it is convenient to rewrite all derivatives of J with respect to the coordinates as derivatives with respect to φ (to lighten the notation, we indicate derivation with respect to φ with a subscript ,φ ).With the new components of the stress-energy tensor (74) and also using Eq. ( 60), the field equations ( 69) become (76a)

Hawking flux and apparent horizon trajectory
We now derive the asymptotic form of the Hawking flux in our model.This can be done by studying the behavior of T µν at future null infinity x + → ∞.In this region, we are considering φ → ∞.We are, therefore, in the decoupling regime, where the gravitational coupling is weak, so that the effects of backreaction can be approximately neglected.This means that the solution in the region of interest (x + > x + 0 ) corresponds to the classical one (61).In this limit, J → 0, J ,φ → 0 and J ,φφ → 0. We have thus (77a) where now indicates differentiation with respect to x − .This result agrees with that of Ref. [62], as it is naturally expected.As it is noted there, this expression diverges once the (outer) event horizon is reached, due to the choice of coordinates adopted.One way to solve this problem is to redefine the x − coordinate as x− ≡ F (x − ) and exploit the anomalous transformation (67).This leads to a well-behaved expression at the horizon, which reads Using the form of F given by Eq. ( 61), we obtain, approaching the horizon The proportionality relation is the same as the SB law (44), which confirms the Planckian nature of the emitted spectrum.Here, however, the flux is modified with respect to standard singular models due to the specific form of the potential (17).
Since the outgoing flux of Hawking radiation is positive, we expect the (outer) apparent event horizon to recede.To see this, we closely follow the approach adopted in Ref. [87].The apparent horizon trajectory x− = x− (x + ) can be derived using the definition of apparent horizon, which satisfies ∂ + φ = 0.This implies from which follows Combining Eqs.(76a), (76b) and (76d) into the above yields The qualitative behavior of the trajectory of the apparent horizon is determined by the sign of the right hand side of the expression above.In order to assess the latter, we would need the full solution of Eqs.(76a) and (76d), which however can only be solved numerically.This will be done in the next section.Here, as a first test, we exploit the fact that the full solution of the field equations approaches the classical one in the asymptotic region φ → ∞, where the coupling with matter fields and backreaction effects become negligible.In this region, the solution is then given by Eq. ( 61).Using this solution into Eq.( 82), it can be shown that dx − /dx + is indeed positive.This implies a receding outer apparent horizon, as expected.Moreover, it confirms the qualitative picture we expect from evaporation, outlined in Section 5 and at the end of Section 6.1: the outer horizon recedes and approaches the horizon of the extremal GS.Notice that, due to the limitations caused by the adopted system of coordinates, we are able to describe the behavior of the outer horizon only.

Numerical results
We now numerically solve the equations of motion given in Eqs.(76).This will allow us to capture the full dynamics of the evaporation process, keeping into account also backreaction effects.A numerical study of evaporating 2D models was performed in the past for the CGHS [88][89][90] and other (singular or regular) 2D models [61,91].
To numerically integrate the equations of motion, we construct a spacetime lattice by means of a grid of null lines and we impose two different sets of boundary conditions, one along x + and one along x − : • At the shock wave, i.e., at x + = x + 0 , we require the solution to coincide with the GS, given by Eq. ( 60); • Above the shock wave, along I − , i.e., at x − → −∞, where backreaction effects are expected to be negligible, we require the solution to match Eq. ( 61), the classical one.Of course, we cannot numerically set a condition at infinity, so we choose a reasonably large negative value.Here we set x − ∞ = −220.This value is found to be the minimal one for which the numerical solution coincides, in the classical limit N → 0, with Eq. ( 61), for every value of x + and x − , within a reasonably small numerical error.For values of x − ∞ greater than −220, the numerical solution deviates from the expected analytical one, while for smaller values the results are the same as the one obtained with x − ∞ = −220, but with a much higher integration time.
Both the boundary conditions at x + = x + 0 and at x − → −∞ are given in implicit form.We therefore first need to integrate and invert the corresponding expressions (the details of this computation are reported in Appendix B).
To numerically integrate the field equations, we also need to select an appropriate integration interval along x + , from the shock wave at x + 0 up to a maximum value x + max .We chose the interval x + ∈ [x + 0 , 5] (where we set x + 0 = 1).We expect the general results to hold also for larger values of x + max .We chose it equal to 5 to have a reasonable integration time interval.For larger values of x + max , the time required to complete a computation increases considerably, given the high computational cost of the algorithm.
The interval on the x + -axis is then discretized into a number n steps of small intervals, with length ∆x = (x + max − x + 0 )/n steps .The number of steps was set equal to n steps = 1000.We checked that, for larger values of n steps , the results of the integration remain qualitatively the same, at the price of having, again, a much longer computational time.
Each point of the discretized x + -axis is labeled by an index i.We choose to discretize the derivatives in the x + direction accordingly, Notice that, with this choice, our algorithm converges to the solution only at first order in ∆x.However, since we are interested in the qualitative behavior of the solutions, Eqs. ( 83) and ( 84) represent a good approximation for the derivatives of φ and ρ.
Along x − at fixed x + , the field equations reduce to ordinary differential equations.Therefore, for each step in the x + direction, we numerically integrate the equations along x − by means of a 4th order Runge-Kutta algorithm.The outcome, thus, is a list of x − -profiles of φ and ρ for each point of the discretized interval on x + (see Fig. 8).
Figure 8: Schematic representation of the numerical algorithm adopted to numerically integrate Eqs.(76a) and (76d).We discretize the x + -axis, and for each interval on the latter, we numerically integrate the field equations along x − using a Runge-Kutta algorithm.
For all the cases considered here, the values of the parameters are set equal to λ = 1 and = 1.The mass of the evaporating solution is fixed equal to M = 0.1 (in these units).
As a first test, we have verified the accuracy of the integration algorithm in the absence of backreaction, i.e., for N = 0, by comparing the numerical solution with the analytical classical one (61).Overall, we find that the relative difference between the numerical and the analytical solutions is smaller than 1 % as long as we consider large negative values of x − , while it increases for x − → 0, staying however Figure 9: Relative difference between the numerical (obtained by setting N = 0) and the exact (61) solutions at x + = x + 0 and x + = 5.As it can be seen, the difference increases as we move towards small negative values of x − .The numerical integration has been carried out on an x + interval of [x + 0 , x + max ] = [1,5], while on the x − direction, in the interval [−220, 0].The parameters of the integration are the following: λ = = 1, M = 0.1 (in these units), n steps = 1000, so that ∆x = 4 • 10 −3 .We checked that, for higher values of n steps , the relative differences decrease, without however altering the qualitative final results and with a much higher computational time.
(see Fig. 9).We checked that increasing n steps leads to an improvement in the accuracy of the integration algorithm (the relative differences decrease), without, however, altering the qualitative final results and at the price of having a much longer computational time.Although a relative discrepancy of 20 % may seem quite important, one should consider that, in the presence of backreaction (see below), the extremal solution is reached at values of x − for which the relative discrepancy stays always below 5 %.As in the following we are not interested in the exact details of the evaporation process, but rather in its qualitative evolution and outcome, we will adopt n steps = 1000 anyway, favouring time efficiency over high precision.
After this preliminary test, we analyse three different cases of increasing N : N = 0, N = 24 and N = 2400, to study the backreaction effects in different regimes.
The x − -profiles of φ and ρ, together with their variations ∆φ = |φ(x − ) − φ ext (x − )|, ∆ρ = |ρ(x − ) − ρ ext (x − )| with respect to the extremal configurations, are computed numerically for several values of x + in the range x + ∈ [x + 0 , 5].At x + = x + 0 , the numerical solutions match exactly the extremal ones, as it should be according to the boundary condition imposed at the shock wave.For simplicity, in Figs. 10 to 12, we only show the plots for x + = 5 and N = 0, 24, 2400.The plots for the other values of x + in the range considered here have the same qualitative behavior.Moreover, given the increase in computational errors near x − ∼ 0, we performed the integration in the range x − ∈ [−220, 0] for convenience.We have checked that the results do not differ from those shown, even if we extend the x − axis to positive values: the convergence to the extremal solution either does not occur in the entire axis or always occurs in the x − < 0 region.
For N = 0, i.e., in the absence of backreaction, we see, as expected, that the black-hole solution remains different from the extremal one for every value of the coordinate x − .
For N = 24, namely when backreaction effects begin to become relevant, we see that, although ∆φ and ∆ρ remain different from zero for all values of x − , they begin to decrease towards zero after reaching a maximum.
For N 24, i.e., N ∼ 2400, when backreaction effects become stronger, we see that ∆φ and ∆ρ become always zero at some finite (negative) value of x − .In general, the larger N , i.e., the stronger backreaction effects, the faster the evaporating configuration reaches the extremal GS.As remarked above, the convergence to the extremal configuration occurs at values of x − for which relative numerical errors are less than 5 %.
It is very important to notice that the convergence of the excited, evaporating solution towards the extremal one is non-monotonic.As one can see clearly from the plots shown (but the same happens also for other values of N not shown here), ∆φ and ∆ρ stay almost flat in the region of large φ (corresponding to x − 0).Then, they reach a sharp maximum at relatively large values of x − before falling rapidly towards zero.This behavior cannot be traced back to backreaction effects, since it is present also in the N = 0 case.The sharp maximum seems to be related to the presence of the maximum in the potential V at (relatively) small values of the dilaton (see Fig. 1), thus to a self-interaction effect of the dilaton.On the other hand, this maximum in V is also responsible for both the presence of two horizons (instead of only one) and for the phase transition small/large black holes (see Sections 2.3 and 4).We will come back to this intriguing point in Section 9.
Summarizing, the numerical integration of eqs.(76) clearly shows that, differently from what obtained in the rough quasistatic description, the effect of the backreaction is to bring the excited, evaporating solution back to the extremal state after a finite time, when N is chosen to be sufficiently large, i.e., at least N ∼ O(10 2 − 10 3 ).

Entanglement Entropy and the Page curve
In this section, we compute the entanglement entropy (EE) of Hawking radiation, described here as a collection of N massless scalar fields, in the 2D nonsingular black-hole geometry.By assuming that the evaporation process is quasistatic, we also determine the time variation of the EE and construct the related Page curve.
The EE of the radiation can be computed by using Kruskal coordinates, covering the region outside the outer event horizon of the black hole, where κ is the surface gravity at the outer event horizon.In these coordinates, the conformal factor of the metric ( 53) can be written as The entanglement entropy of N massless scalar fields in two spacetime dimensions on a line can be evaluated by tracing out the degrees of freedom in a spacelike slice [x, y] connecting two points.The resulting expression is 6 (see, e.g., Refs.[27,28,92,93]) where d(x, y) is the geodesic distance between x and y.In principle, Eq. ( 87) is valid for a QFT on a flat spacetime [94], but it has been generalized to static curved spacetime [92], where d 2 (x, y) reads To compute the entanglement entropy, we construct a spacelike surface encompassing different regions of the black hole.In Fig. 13, this surface is Σ L ∪ I ∪ Σ R , where Σ L and Σ R are two hypersurfaces on the outside regions of the two copies of the black hole, where an observer collects Hawking radiation.They are the portion of the hypersurface, where the radiation degrees of freedom are defined.They are anchored to two timelike surfaces (dashed black lines in Fig. 13 (left wedge).The surface J defines, instead, the interior region of the black hole.The radiation quantum state over the whole hypersurface Σ L ∪ I ∪ Σ R is pure.When tracing out the interior degrees of freedom in I, we obtain the mixed state of the radiation described by the density matrix ρ rad , which can therefore be used to compute the entanglement entropy.This is reminiscent of the thermofield double state of the black hole [31,95]: the entanglement entropy takes into account the correlations between the two disjointed copies of the black hole (right and left wedges).Since we have radiation outside the black hole, there will be two copies of this thermal bath (the two regions Σ L and Σ R ).
In our case, Eq. ( 88) reads This is valid off-shell.On shell κ → κ H , we have βκ H 2 = π, and thus Finally, the EE of the matter fields is As stressed above, this is valid as long as we consider the static case.However, κ H varies due to the evaporation process.To get a qualitative picture of the behavior of the entropy in time, we can assume that the evaporation process happens in an adiabatic way, so that we can use a quasistatic approximation.The evaporation is thus again modelled in terms of a sequence of static states with decreasing mass.As we have seen explicitly in Section 5.3, the quasistatic approximation is reliable as long as the semiclassical one is valid.In a first approximation, therefore, we can use the time variation of the event horizon r H = r H (t), computed as a solution of the SB law (45), and plug it into the expression of the surface gravity The qualitative result (obtained neglecting the irrelevant constants N 6 ln [4f (b)]) is plotted in Fig. 14.As in singular black-hole models, initially the entanglement entropy of the radiation grows.However, this growth reaches a maximum at the "Page time" t P and then the entropy starts decreasing, due to the peculiar form of the surface gravity, which is related to the absence of a singularity.
It is interesting to note that t P physically coincides with the onset time of the second order phase transition (graphically, the intersection point between the solid blue and the horizontal dotted lines in Fig. 6(b)).This feature was found before for nonsingular black holes in Ref. [36], where it was noted that the presence of the dS core traps Hawking modes, which cause a decrease in entropy once freed from the trapping region.This indeed happens as we get closer to the extremal configuration, right after the onset of the second order phase transition, as the role of the inner horizon becomes increasingly important.This release of information could also be related to the peculiarities of the latter, which has negative surface gravity, causing an outburst of energy in the final stages of the evaporation [37], a process similar to the mass inflation.91), where we considered the time variation of the surface gravity, calculated using the solution of Eq. ( 45).The vertical, dashed green line corresponds to the maximum of the curve (the Page time t P ), while the vertical dashed orange one indicates the time t B at which the semiclassical approximation should break down.
The mechanism described above is qualitative similar to that taking place in the island proposal [26][27][28][29] (for an application to two-horizon models, both singular and regular, see, e.g., Refs.[96][97][98]; for an application to dS spacetime, see, e.g., Refs.[99,100]), where the resolution of the information paradox is traced back to a transition in the behavior of the entanglement entropy functional.At late times (right after the Page time), this functional starts receiving contributions from nontrivial configurations in the black hole interior (the "islands").This allows to correctly keep track of the entanglement structure of both the black hole and the radiation subsystems and to reconstruct the Page curve.Our approach provides, at qualitative level, a physical explanation for the appearance of such configurations in the black-hole interior.In our description, the islands correspond to the inner horizon, while the transition in the entropy functional is physically realized by the second order phase transition 7 .In this paper we have only considered nonsingular black holes with two horizons, characterized therefore by a timelike singularity.It is, therefore, currently unclear if and how the results of the present paper compare with the island rule applied to black holes possessing spacelike singularities (see, e.g., Refs.[101][102][103]) 8 .In general, we expect the single-horizon case to be qualitatively different from the two-horizon one.On the other hand, at least in some particular cases, the behavior of regular models with a single event horizon could be not so drastically different from that of two-horizon black holes (see, e.g., Refs.[58,59,79,104]).
The assumptions used so far are the validity of the quasistatic and the semiclassical approximations.As we have seen in Section 5.3, the semiclassical approximation (and hence also the quasistatic one) breaks down near extremality, when we reach the energy gap (52).This happens at the time corresponding to the vertical dashed orange line in Fig. 14.Therefore, we have to cut the Page curve when the semiclassical approximation breaks down.What happens beyond this point cannot be inferred from our semiclassical description of the dynamics.In particular, we cannot assess whether the decrease in EE continues until it becomes zero at extremality, as it would be expected for an evaporation process that leaves behind a quantum pure state.Results from AdS 2 quantum gravity indicate the occurrence of a quantum phase transition from the LDS vacuum to the AdS 2 CDV [68].Similarly to what happens in the case of extremal charged black holes, the near-extremal, near-horizon state of 4D nonsingular Hayward black holes, described by the AdS 2 × S 2 spacetime, could have a purely topological entropy content, explained in terms of AdS 2 fragmentation [32].

Conclusions
In this paper we have investigated the thermodynamics and the classical and semiclassical dynamics of 2D, AF, nonsingular dilatonic black holes with a dS core.The aim has been trying to understand both the end point of the black-hole evaporation and the related information flow.Our analytic and numerical results provide some evidence that the latter leads to a regular extremal configuration in a finite amount of time.This conclusion is supported both by a thermodynamic analysis, showing that extremal configurations are energetically preferred, and by the numerical integration of the semiclassical field equations, which allows to take into account the full dynamics of the backreaction.Concerning the information flow during the evaporation process, the Page curve we have constructed clearly shows the presence of a maximum at the Page time, followed by a descent, in which information is recovered from the hole.These features have a nice physical explanation in terms of the trapping/flow of Hawking modes in the region between the inner and outer horizons.Additionally, they may provide some physical insight on the mechanism underlying the recently advanced "island" proposal to address the information paradox.
On the other hand, the intrinsic limitation of our approach, the validity of the semiclassical approximation, prevents us from putting a definitive final word on the issue.Near extremality, the semiclassical approximation breaks down and a mass gap is generated, which is a well-known feature of AdS 2 quantum gravity.This gap separates the extremal configuration, endowed with a linear dilaton, from another vacuum of the theory, the AdS 2 CDV.There are strong indications that, once the former is reached during evaporation, a quantum phase transition occurs that causes a transition to the latter.A similar phase transition between the extremal LDS and the AdS 2 CDV seems to occur also in the case under consideration.If this could be independently confirmed, it would imply that the end point of the evaporation process is the full AdS 2 spacetime endowed with a constant dilaton, which is a perfectly regular configuration.This would be a clear confirmation that the evaporation process is unitary.
In terms of the 4D solution, which we are here modelling with the 2D dilaton gravity models, this corresponds to a (topology changing) phase transition between an extremal solution, which is an AF spacetime with an infinitely long throat, and an AdS 2 × S 2 spacetime, describing instead the near-horizon region.
Our investigation has also revealed some intriguing peculiar features of the semiclassical dynamics of regular black holes with a dS core.The approach to extremality of excited black-hole solutions is not monotonic, but presents a sharp maximum of ∆ρ and ∆φ near extremality, at relatively small values of the dilaton.The presence of this maximum can be explained in terms of the self-interaction of the dilaton, i.e., the presence of a maximum in the dilaton potential V (φ).This, in turn, is what determines the presence of two horizons, instead of a single one, the presence of a large/small-black-hole phase transition and also the decreasing tail in the Page curve for the EE of Hawking radiation.Since the behavior of the dilatonic potential is essentially determined by the presence of a dS core, all the previous features are a consequence of the absence of singularities, which is achieved exactly thanks to this dS behavior in the black-hole interior.This strongly supports the hypothesis that the nonunitary evolution of the evaporation process has to be traced back to the presence of a singularity in the black hole interior.

Acknowledgements
We thank D. Grumiller for valuable comments and for having pointed our attention to Ref. [69].
A No divergences in the ground state stress-energy tensor In Ref. [62], counterterms are added to the Polyakov action (as in Eq. ( 71)) to eliminate a divergence in the stress-energy tensor of the GS T µν ∼ φ 2a , which, for a > 0 diverges at asymptotic infinity φ → ∞.We show that in our case the rapid fall of the potential at infinity allows to avoid this divergent behavior.
We consider, as an example, T +− GS , but similar considerations also hold for the other components of T µν .Using Eq. ( 64), we have where ρ, computed at the GS, is given by Eq. ( 60) Therefore where J ,φ = λ −2 V and we used the vacuum solution Eq. (55b).Differentiating with respect to x + yields For φ → ∞, e 2ρ → constant, while V ,φ ∼ −φ −3 → 0. So we do not have divergences.

B Boundary conditions for numerical integration
B.1 Boundary condition at x + = x + 0 For our numerical integration, we set λ = 1 and = 1.
When the shock wave is turned on, the solution is set equal to the extremal configuration (60).Our procedure to implement this boundary condition is the following: 1. We first integrate the equation for the dilaton (expression on the right of Eq. ( 60)), leading to an implicit relation between the dilaton and the coordinates; 2. We numerically invert it to have explicitly φ = φ(x + , x − ); 3. We plug the result into the equation for e 2ρ .
The first point is achieved by solving the differential equation This integral can be done analytically.With M ext = 3 √ 2 /3 and = λ = 1 it reads r * ,ext = 2 5/3 ln We then numerically invert the result to get φ = φ(x + , x − ) and plug the result into Eq.( 60) to get ρ = ρ(x + , x − ), both evaluated at extremality.
Again, we numerically invert this expression, to get φ = φ(x + , x − ).In this case, however, the quantity we are (improperly) calling r * contains the function F (x − ), which has to be computed to fully obtain φ as a function of the coordinates.Starting from we note that, when evaluated at x + = x + 0 , it reads from which we get This expresssion of φ(x + 0 , x − ) can be used in Eq. ( 61) to compute J 0 .This gives us a differential equation in terms of F which reads This equation is solved numerically.The integration constant is chosen so that, once the solution of F is plugged into the dilaton solution (61), at x + = x + 0 , the dilaton is equal to the extremal one computed at x + 0 .This guarantees the continuity of the scalar field across the shock wave.Once F (x − ) is computed, we plug it into φ = φ(x + , x − ).With this and Eq. ( 61), we obtain also ρ = ρ(x + , x − ) above extremality.

4 Black-hole thermodynamics 4 . 1
Thermodynamic potentials and the first principle

Figure 2 :
Figure 2: Qualitative behavior of the temperature of the 2D black hole, according to Eq. (26).

Figure 3 :
Figure 3: Qualitative behavior of the specific heat of the 2D black hole.The vertical dashed line corresponds to the peak of the potential, φ peak .

Figure 4 :
Figure 4: Temperature T (φ) as a function of φ.Here, we highlighted the two configurations corresponding to φ 1 and φ 2 , respectively in the stable and unstable branch and their common temperature T = T (φ 1 ) = T (φ 2 ).The areas representing (φ 2 − φ 1 )T and the right-hand side of Eq. (38), are coloured in blue and red, respectively.

Figure 5 :
Figure 5: Figure (a): Numerical solution of Eq. (44) (blue solid line), which shows the time evolution of the mass of an evaporating black hole as a function of time, in the quasistatic approximation.We see that at large times, the numerical solution asymptotes the extremal mass (orange dashed line).For the numerical integration, we set M(t = 0) = 0.1 λ as an initial condition.Figure (b): Time variation of the entropy of the black hole according to Eq.(33).At late times, the entropy goes to zero, after subtracting the topological term S 0 .In both figures, we set λ = = 1.
(a) shows the time variation of ∆E = |M − M ext |: the intersection with the horizontal dashed line, corresponding to ∆E ∼ E gap , identifies the time at which the semiclassical approximation breaks down.

Figure 6 :
Figure 6: Figure (a): Plot of ∆E as a function of time (solid blue line).The horizontal dashed orange line corresponds to the value of the mass gap (52).The intersection point between the two curves gives the time where the semiclassical approximation breaks down.Figure (b): Mass evolution in time (solid blue line).The horizontal dashed orange line corresponds to the value of the mass of the extremal configuration.The horizontal dotted green line corresponds to the value of the mass at the temperature peak(27).The vertical dashed line, instead, signals the instant of time where |M − M ext | ∼ E gap , with E gap given by Eq. (52).As it can be seen, the breakdown of the semiclassical approximation happens after the evaporation process reaches the temperature peak.For both figures, we set M(t = 0) = 0.1 λ and = λ = 1.

Figure 10 :
Figure 10: Upper figures: Comparison between the numerical (solid blue line) and analytical extremal (dashed orange line) dilaton solutions (left figure), and difference between the two (right figure), as functions of x − .Lower figures: Comparison between the numerical (solid blue line) and analytical extremal (dashed orange line) metric solution (left figure), and difference between the two (right figure), as functions of x − .All figures are evaluated at x + = 5 and with N = 0, in units where λ = = 1.

ϕFigure 11 :Figure 12 :
Figure 11: Upper figures: Comparison between the numerical (solid blue line) and analytical extremal (dashed orange line) dilaton solutions (left figure), and difference between the two (right figure), as functions of x − .Lower figures: Comparison between the numerical (solid blue line) and analytical extremal (dashed orange line) metric solution (left figure), and difference between the two (right figure), as functions of x − .All figures are evaluated at x + = 5 and with N = 24, in units where λ = = 1.

Figure 13 :
Figure 13: Penrose diagram of the maximally extended spacetime of the nonextremal configuration.The two points b + and b − , belonging to the right and left wedges, respectively, are highlighted, and represents points anchored to two timelike curves (dashed black lines in the two wedges).The union between the three hypersurfaces Σ L ∪ I ∪ Σ R (the red and orange curves) is a spacelike surface and the state defined on it is pure.The radiation is defined on Σ L and Σ R and its state is mixed.

Figure 14 :
Figure14: Qualitative time variation of the entanglement entropy of matter fields, according to Eq. (91), where we considered the time variation of the surface gravity, calculated using the solution of Eq. (45).The vertical, dashed green line corresponds to the maximum of the curve (the Page time t P ), while the vertical dashed orange one indicates the time t B at which the semiclassical approximation should break down.