We provide a modesum prescription to directly compute the renormalized stressenergy tensor RSET for scalar fields in the Boulware vacuum. The method generalizes the recently developed extended coordinate method which was previously only applicable to HartleHawking states. We exhibit the accuracy and efficiency of the method by calculating the RSET in subextremal and extremal ReissnerNordstrm spacetimes. We find numerical evidence for the regularity of the RSET at the extremal horizon regardless of the field mass and its coupling. We employ our numerical results of the RSET to source the semiclassical Einstein equations demonstrating that if the RSET is considered as a static perturbation it will either deextremalize the black hole or convert it into a horizonless object.
['General Relativity and Quantum Cosmology', 'High Energy Physics - Theory', 'Quantum Physics']
The renormalized stressenergy tensor for scalar fields in the Boulware state with applications to extremal black holes
{'The renormalized stress-energy tensor for scalar fields in the Boulware state with applications to extremal black holes': 'Julio Arrechea ∗  \nIFPU, Institute for Fundamental Physics of the Universe, via Beirut 2, 34014 Trieste, Italy SISSA, International School for Advanced Studies, via Bonomea 265, 34136 Trieste, Italy and INFN Sezione di Trieste, via Valerio 2, 34127 Trieste, Italy', 'Cormac Breen †': 'School of Mathematics and Statistics, Technological University Dublin, Grangegorman, Dublin 7, Ireland', 'Adrian Ottewill ‡': 'School of Mathematics and Statistics, University College Dublin, Belfield, Dublin 4, Ireland  \nLorenzo Pisani § and Peter Taylor ¶  \nCenter for Astrophysics and Relativity, School of Mathematical Sciences, Dublin City University, Glasnevin, Dublin 9, Ireland (Dated: September 10, 2024)  \nWe provide a mode-sum prescription to directly compute the renormalized stress-energy tensor (RSET) for scalar fields in the Boulware vacuum. The method generalizes the recently developed extended coordinate method which was previously only applicable to Hartle-Hawking states. We exhibit the accuracy and efficiency of the method by calculating the RSET in sub-extremal and extremal Reissner-Nordstrom spacetimes. We find numerical evidence for the regularity of the RSET at the extremal horizon regardless of the field mass and its coupling. We employ our numerical results of the RSET to source the semi-classical Einstein equations, demonstrating that if the RSET is considered as a static perturbation, it will either de-extremalize the black hole, or convert it into a horizonless object.', 'I. INTRODUCTION': "In the semiclassical approximation to quantum gravity, one replaces the stress-energy tensor in the gravitational field equations with the corresponding expectation value for the quantized matter fields, including ostensibly quantized graviton fields considered as linear perturbations about a background geometry. The gravitational field equations in this framework are  \nG µν +Λ g µν + α 1 H (1) µν + α 2 H (2) µν = 8 π ⟨ ˆ T µν ⟩ ren , (1)  \nwhere g µν is the metric of spacetime, Λ is the cosmological constant, G µν is the Einstein tensor and ⟨ ˆ T µν ⟩ ren is the renormalized expectation value of the stress-energy tensor (RSET) of the quantum fields in some quantum state. The tensors H (1) µν , H (2) µν are geometrical and include terms that are quadratic in the curvature; their inclusion arises through the point-splitting regularization process that yields ⟨ ˆ T µν ⟩ ren . This regularization process allows for the calculation of the finite physical RSET with an infinite renormalization of the constants Λ, α 1 and α 2 . It is difficult to conceive how to practically solve Eq. (1) in  \na general self-consistent way since the right-hand side depends on the unknown solution g µν , and yet to compute the stress-energy tensor explicitly requires that we know the geometry of spacetime a priori . It is only by considering highly symmetric spacetimes [1] or by subjecting the RSET to additional approximations [2-4] that the backreaction problem becomes tractable. In black hole spacetimes, one usually proceeds by a perturbative expansion of the metric about a background classical solution. This leads to the reduced order semiclassical equations which describe, to first order in ℏ , the quantum backreaction induced by a quantum field on the classical geometry [5-7]. In these equations, the source term is the renormalized expectation value for a quantum field in a given quantum state on a fixed classical background spacetime.  \nIt turns out that in most situations of interest, even the calculation of the renormalized expectation value of the stress-energy tensor (RSET) for a particular quantum state on a fixed background is a technically challenging endeavour. For black hole spacetimes, there are three main approaches to this calculation, the CandelasHoward approach [8] and its extensions (see for example [9-11]), the so-called 'pragmatic mode-sum prescription' [12-14] and the method developed by authors of this paper called the 'extended coordinate method' [15-18]. While all these methods have different approaches and advantages, one drawback of both the original CandelasHoward approach and the extended coordinate approach is that they were developed specifically for the field in the Hartle-Hawking state [19]. Working with this state  \noffers a computational advantage since one can employ Euclidean techniques that enforce a discretization of the frequency spectrum. Anderson, Hiscock and Samuel [10] generalized the Candelas-Howard method to include a direct method for computing the RSET for a scalar field in the Boulware state [20]. Although the RSET in the Unruh, Boulware and Hartle-Hawking state were computed in Ref. [18], only the Hartle-Hawking state was directly renormalized, the RSET for the other two states were computed by leveraging the fact that the difference between two quantum states is regular. In this paper, we generalize the extended coordinate method to directly compute the RSET for a scalar field in the Boulware state.  \nA natural question then is: if the RSET for a field in the Boulware state can be computed by first renormalizing in the Hartle-Hawking state and employing a state subtraction scheme, then why do we need a direct method for renormalizing in the Boulware state? There are several justifications for a direct renormalization prescription for Boulware. First, when employing the indirect method one needs to compute both the Euclidean modes (to do the renormalization in the Hartle-Hawking state) and the Lorentzian modes (to do the state subtraction scheme). Employing a direct method for Boulware requires only one set of modes and hence is computationally less expensive. Second, there are black hole spacetimes that do not admit a Hartle-Hawking state, for example, the extreme Reissner-Nordstrom spacetime, and hence one requires a direct method for computing the RSET in such situations. Finally, moving away from black hole spacetimes, if we consider semiclassical effects in the spacetime of a star, then the natural vacuum state is the Boulware state (there is, in fact, no other reference state). Hence this paper sets the foundation for renormalization in those contexts.  \nThe approach we take is similar to that taken in the extended coordinate method for the Hartle-Hawking state [15-17]. We take an arbitrarily high-order expansion of the Hadamard parametrix which encodes the shortdistance behaviour of the two-point function, and we rewrite the terms in this expansion as mode-sums of the same form as the Boulware propagator. While the spirit is the same, the mathematical development is significantly trickier because the frequency spectrum is continuous, rather than the discrete frequency spectrum that one obtains in the Euclidean approach to defining the Hartle-Hawking state. As well as presenting a computational challenge, the continuous frequency spectrum also means that the Fourier transform of the terms in the Hadamard parametrix, which is needed to express these as mode-sums, yields functions of the frequency that are not ordinary functions but rather generalized functions. When we treat these mode-sums rigorously as generalized functions, it becomes completely transparent how an infra-red cutoff appears in these ω integrals and indeed it is completely explicit that the RSET is independent of the choice of this cutoff. A similar cutoff was introduced  \nin Ref. [10] but it was introduced by necessity to eliminate an infra-red divergence in the WKB approximation. While the approaches are ultimately equivalent, it is clear here that the source of the cutoff is a consequence of treating rigorously the generalized functions that arise in the Fourier transform of the Hadamard parametrix.  \nAs a consistency check, we contrasted our numerical results with previously known results for the sub-extremal and extremal Reissner-Nordstrom spacetime [18, 21] and found excellent agreement. We pay particular attention to the regularity of the RSET at the horizon of the extreme Reissner-Nordstrom black hole. Numerical evidence suggesting its regularity was found in [21] for massless fields (see [22-24] for discussions about backreaction effects). A clear advantage of the extended coordinate method over some previous approaches is that it does not invoke the Wentzel-Kramers-Brillouin (WKB) approximation in any step of the calculation, allowing to generate very accurate results for the RSET near the horizon. Our results reinforce the existing numerical evidence of the regularity of the RSET in extremal black holes, and extends those results to massive fields as well.  \nThis paper is organized as follows. In Sec. II, we show how we express the expansion of the Hadamard parametrix as a mode-sum of the same type as the Boulware propagator for a scalar field in a static, spherically symmetric spacetime. In Sec. III, we give the explicit expressions for the renormalized expectation values for the vacuum polarization and the stress-energy tensor for a scalar field in the Boulware state. In Sec. IV, we apply the new method to compute the RSET for a scalar field on the Reissner-Nordstrom spacetime, including the extremal case where a direct calculation is necessary. In Sec. V, we employ our accurate RSET results to solve the reduced-order semiclassical equations to compute the backreaction on the extremal Reissner-Nordstrom black hole spacetime. Finally, in Sec. VI, we conclude with some discussion and future prospects.", 'II. RENORMALIZATION PRESCRIPTION IN THE BOULWARE STATE': "In this section, we will briefly outline the extended coordinate approach to calculating the RSET for a quantum scalar field in the Boulware quantum state [20], propagating on a static, spherically symmetric black hole spacetime. The Boulware state is vacuum according to a stationary observer far from the black hole. It is useful to treat the Boulware state as a thermal state at zero temperature and to employ Euclidean methods. We consider Euclidean line elements of the form:  \nds 2 = f ( r ) dτ 2 + dr 2 /f ( r ) + r 2 ( dθ 2 +sin 2 θdϕ 2 ) . (2)  \nIf this line element described a black hole spacetime with event horizon r = r + and surface gravity κ = 1 2 f ' ( r + ), then for the Hartle-Hawking state we would impose the periodicity condition τ = τ +2 π/κ + on the Euclidean  \ntime. Imposing this periodicity discretizes the frequency spectrum of the field modes which is a major computational advantage. However, in the Boulware case, the field temperature is zero and so no such identification is made in the Euclidean time. As a result, the Euclidean geometry for a non-extremal black hole contains a conical singularity at the horizon r = r + and indeed the quantum expectation values will diverge there as a result. Notwithstanding that we do not have a discrete frequency, the Euclidean approach is still useful here since the Green function does not require an ' i ϵ ' prescription to avoid singularities on the lightcone; there are no divergences on the lightcone since there is no lightcone in a Euclidean space. In particular, the wave equation now satisfies an elliptic equation  \n( □ -µ 2 -ξ R ) ϕ = 0 , (3)  \nwhere □ is the d'Alembertian operator with respect to the Euclidean metric, µ is the field mass with dimensions of inverse length, R is the Ricci curvature scalar of the background spacetime and ξ is the coupling strength between the field and the background geometry. The corresponding Euclidean Green function has the following mode-sum representation (with r = r ' for simplicity) on this black hole spacetime,  \nG ( x, x ' ) = 1 8 π 2 ∞ ∑ l =0 (2 l +1) P l (cos γ ) ∫ ∞ -∞ dω e iω ∆ τ g ωl ( r ) , (4)  \nwhere ∆ x ≡ x ' -x ∼ O ( ϵ ) is the coordinate separation, γ is the geodesic distance on the 2-sphere and P l ( z ) is the Legendre polynomial of the first kind. We have denoted by g ωl ( r ) := p ωl ( r ) q ωl ( r ) /N ωl the one-dimensional radial Green function evaluated at the same radius. The radial modes p ωl ( r ), q ωl ( r ) are solutions of the homogeneous equation:  \n[ d dr ( r 2 f ( r ) d dr ) -r 2 ( ω 2 f ( r ) +( µ 2 + ξ R ) ) -l ( l +1) ] χ ωl ( r ) = 0 , (5)  \nwith p ωl ( r ) the solution regular at the horizon and q ωl ( r ) regular at the outer boundary (usually spatial infinity or the cosmological horizon). The normalization constant is given by  \nN ωl = -r 2 f ( r ) W{ p ωl ( r ) , q ωl ( r ) } , (6)  \nwhere W{ p, q } denotes the Wronskian of the two solutions.  \nIn the coincidence limit ∆ x → 0 (i.e. γ → 0 and ∆ τ → 0), the mode sum (4) diverges. To renormalize this mode sum, we will adapt the so-called extended coordinate implementation of Hadamard renormalization developed in Refs. [15-17]. For a quantum state satisfying the Hadamard condition, the local singularity structure  \nof the two-point function when x and x ' are sufficiently close is encapsulated in the Hadamard parametrix  \nK ( x, x ' ) = 1 8 π 2 ( ∆ 1 / 2 ( x, x ' ) σ ( x, x ' ) + V ( x, x ' ) log ( 2 σ ( x, x ' ) ℓ 2 ) ) (7)  \nwhere σ ( x, x ' ) is Synge's world function which measures half the square of the geodesic distance between x and x ' , ∆ 1 / 2 ( x, x ' ) and V ( x, x ' ) are symmetric biscalars that are regular in the coincidence limit x ' → x . Moreover, these biscalars depend only on the local geometry and the field parameters, not on any global considerations such as the quantum state. The lengthscale ℓ here is an arbitrary lengthscale associated with the renormalization ambiguity. It arises since V ( x, x ' ) is a solution to the homogeneous wave equation and so one is free to add arbitrary factors of V ( x, x ' ) to K ( x, x ' ) and it will remain a parametrix for the wave operator. Since σ ( x, x ' ) → 0 as x ' → x , it is clear that K ( x, x ' ) is singular at coincidence.  \nThe fact that the mode-sum (4) does not converge at coincidence is because the two-point function possesses singularities precisely of the type in Eq. (7) in the shortdistance limit. The goal is to subtract K ( x, x ' ) from G ( x, x ' ) in such a way that the coincidence limit can be meaningfully taken. The only way to proceed in this approach is to express K ( x, x ' ) as a mode-sum of the same type as the mode-sum representation of Eq. (4) and then to subtract mode-by-mode. This renormalized mode-sum will converge in the coincidence limit.  \nExpressing the Hadamard parametrix as a mode-sum is non-trivial. The approach of Refs. [15, 16] is to expand the parametrix in a judiciously chosen set of extended coordinates. In that context, it was assumed the field was in the Hartle-Hawking state and the extended coordinates were adapted to that state by ensuring the time dependence had the appropriate thermal periodicity. This must be amended here for the Boulware state. We define  \ns 2 = f ( r ) ∆ τ 2 +2 r 2 (1 -cos γ ) , (8)  \nand then, rather than expand the parametrix in coordinate separation ∆ x , we expand in s and ∆ τ (we take radial points together ∆ r = 0). The result of expanding the parametrix up to O ( ϵ 2 m log ϵ ) in these new extended  \ncoordinates is  \nK ( x, x ' ) = 1 8 π 2 ( m ∑ a =0 a ∑ b =0 D (r) ab ( r ) ∆ τ 2 a +2 b s 2 b +2 + m ∑ a =1 a ∑ b =1 D (p) ab ( r )∆ τ 2 a -2 b s 2 b -2 + m -1 ∑ a =1 a -1 ∑ b =0 T (r) ab ( r ) ∆ τ 2 a +2 b +2 s 2 b +2 + m -1 ∑ a =0 a ∑ b =0 T (l) ab ( r ) s 2 a -2 b ∆ τ 2 b log ( s 2 ℓ 2 ) + m -1 ∑ a =1 a ∑ b =0 T (p) ab ( r ) s 2 a -2 b ∆ τ 2 b ) + O ( ϵ 2 m log ϵ ) . (9)  \nThe coefficients D (r) ab ( r ) and D (p) ab ( r ) come from the expansion of the direct part of the parametrix U/σ and correspond to the coefficients of terms in the expansion that are rational in ∆ τ 2 and s 2 and those that are polynomial in ∆ τ 2 and s 2 , respectively. Similarly, T (r) ab ( r ) and T (p) ab ( r ) are the coefficients of terms that are rational and polynomial coming from the tail V log(2 σ/ℓ 2 ) of the parametrix, while T (l) ab ( r ) is the coefficient of terms in the tail that contain a logarithm. Expressions for the first few of these coefficients are given in Appendix C.  \nFormally speaking, this expansion need only be computed up to m = 1 in order to remove the divergences from the Green function. Nevertheless, by expanding up to higher order terms, including terms that formally vanish in the coincidence limit, and expressing these terms as a mode-sum permits us to have precise control over the rate of convergence of the renormalized mode sum. On the other hand, if we only subtract the singular terms, the mode sum would indeed converge, however, it would converge slowly. Furthermore, the distinction between the different types of terms that occur in the expansion is an important one. In particular, the terms that are polynomial in the variables ∆ τ 2 and s 2 do not require a mode-sum representation. For example, a multipole expansion of terms that are polynomial in s 2 would involve finite multipole moments and hence cannot affect the large l behaviour. Hence the terms above that are polynomial in the expansion parameters are left as-is. The rational and logarithmic terms will be expressed as mode sums. This is the focus of the next section.", 'A. Regularization Parameters for the Direct Part': 'Starting with the direct part of the Hadamard parametrix, which involves terms of the form ∆ τ 2 a +2 b /s 2+2 b for a ≥ b ≥ 0. We take as an ansatz for  \nthe mode-sum representation of this term the following  \n∆ τ 2 a +2 b s 2+2 b = ∞ ∑ l =0 (2 l +1) P l (cos γ ) ∫ ∞ -∞ dωe iω ∆ τ Ψ ωl ( a, b | r ) , (10)  \nwhere the generalized functions Ψ ωl ( a, b | r ) we call the regularization parameters for the direct part. The goal is to invert the above ansatz and hence obtain an explicit closed-form representation for the regularization parameters. A double integral expression for the regularization parameters follows immediately from the completeness of the Legendre polynomials and the Fourier transform of the delta distribution, yielding  \nΨ ωl ( a, b | r ) = 1 4 π ∫ 1 -1 dx ∫ ∞ -∞ dt e -iωt t 2 a +2 b P l ( x ) ( f t 2 +2 r 2 (1 -x )) 1+ b , (11)  \nwhere we have relabeled the integration variables by x = cos γ and t = ∆ τ . If we now define  \nρ 2 = f 2 r 2 , (12)  \nthen we can rewrite our integral expression above as  \nΨ ωl ( a, b | r ) = 1 4 π 1 (2 r 2 ) b +1 ∫ 1 -1 dxP l ( x ) ( -1) b b !2 b × ( 1 ρ ∂ ∂ρ ) b ∫ ∞ -∞ dt e -iωt t 2 a ( ρ 2 t 2 +1 -x ) . (13)  \nNow it is this t -integral that distinguishes the regularization parameters in the Hartle-Hawking and the Boulware states. In the former, the t -integral above would be over a finite domain on account of the time periodicity of a thermal state. But in this case, the integral is over the entire real line which means that for a ≥ 1, the integral does not converge in the sense of ordinary functions and must be considered in the sense of generalized functions. The most straightforward way to proceed is to start with the a = 0 integral which is an ordinary function and is given by  \n∫ ∞ -∞ dt e -iωt ( ρ 2 t 2 +1 -x ) = π e -| ω | √ 1 -x/ρ ρ √ 1 -x . (14)  \nFor a ≥ 1, we can take (2 a )-derivatives with respect to ω of the result above. Since the exponent on the right-hand side involves | ω | and not ω , we pick up delta distributions and their derivatives. The result is  \n∫ ∞ -∞ dt e -iωt t 2 a ( ρ 2 t 2 +1 -x ) = π ( -1) a e -| ω | √ 1 -x/ρ (1 -x ) a -1 / 2 ρ 2 a +1 -2 π ( -1) a a ∑ p =1 (1 -x ) a -p ρ 2 a -2 p +2 δ (2 p -2) ( ω ) , (15)  \nwhich is clearly not an ordinary function. Substituting this back into (13) gives  \nΨ ωl ( a, b | r ) = 1 (2 r 2 ) b +1 ( -1) a + b b !2 b +2 ( 1 ρ ∂ ∂ρ ) b { 1 ρ 2 a +1 A ωl ( a | ρ ) -2 a ∑ p =1 δ (2 p -2) ( ω ) ρ 2 a -2 p +2 B l ( a, p ) } , (16)  \nwhere  \nA ωl ( a | ρ ) = ∫ 1 -1 P l ( x ) e -| ω | √ 1 -x/ρ (1 -x ) a -1 / 2 dx, (17) B l ( a, p ) = ∫ 1 -1 P l ( x )(1 -x ) a -p dx. (18)  \nThe second integral here is a standard one. The result is that integral vanishes whenever l ≥ a -p + 1 and  \notherwise gives  \nB l ( a, p ) = ( -1) l 2 a -p +1 Γ( a -p +1) 2 Γ( a -p + l +2)Γ( a -p +1 -l ) . (19)  \nThe A ωl ( a | ρ ) integral can also be easily performed and there are several different equivalent representations. The tidiest form involves derivatives of the modified Bessel functions  \nA ωl ( a | ρ ) = 2 a +3 / 2 ∂ 2 a ∂z 2 a I l +1 / 2 ( z/ 2) K l +1 / 2 ( z/ 2) , (20)  \nwhere  \nz = | ω | √ 2 ρ = | ω | r √ f . (21)  \nHowever, this is not the best representation for numerically computing this object. For numerical purposes, except when z >> 1, we find the following representation much more efficient  \nA ωl ( a | ρ ) = 2 a √ 2 π ( -1) l { Γ( a + 1 2 ) 2 2 ˆ F 3 ( { a + 1 2 , a + 1 2 } { a + l + 3 2 , a -l + 1 2 , 1 2 } ω 2 2 ρ 2 ) -Γ( a +1) 2 | ω | √ 2 ρ 2 ˆ F 3 ( { a +1 , a +1 } { 3 2 , a + l +2 , a -l +1 } ω 2 2 ρ 2 ) } , (22)  \nwhere 2 ˆ F 3 is the regularized general hypergeometric function [25]. Now using the fact that,  \n( 1 ρ ∂ ∂ρ ) b ρ -α = ( -1) b 2 b Γ( α/ 2 + b ) Γ( α/ 2) ρ -α -2 b , (23)  \nwe can easily express the derivatives in terms of hypergeometric functions. Putting all of the above together, and after some simplification, we obtain the following expression for the regularization parameters  \nΨ ωl ( a, b | r ) = I ωl ( a, b | r ) -(2 r 2 ) a f a + b +1 ( -1) a + l b ! a ∑ p =1 2 a -p ( a -p )!( a -p + b )! ( a -p + l +1)!( a -p -l )! ( f 2 r 2 ) p δ (2 p -2) ( ω ) , (24)  \nwhere, for later convenience, we have defined  \nI ωl ( a, b | r ) = (2 r 2 ) a -1 / 2 f a + b +1 / 2 2 a -2 √ 2 π ( -1) a + l Γ( b +1) { Γ( a + 1 2 )Γ( a + b + 1 2 ) 2 ˆ F 3 ( { a + 1 2 , a + b + 1 2 } { a + l + 3 2 , a -l + 1 2 , 1 2 } ω 2 r 2 f ) -Γ( a +1)Γ( a + b +1) | ω | r √ f 2 ˆ F 3 ( { a +1 , a + b +1 } { 3 2 , a + l +2 , a -l +1 } ω 2 r 2 f ) } . (25)  \nSubstituting this back into (10) and performing the ω integration in the terms involving delta distributions gives  \n∆ τ 2 a +2 b s 2+2 b = ∞ ∑ l =0 (2 l +1) P l (cos γ ) ∫ ∞ -∞ dω e iω ∆ τ I ωl ( a, b | r ) + 1 f b +1 a ∑ p =1 ( a -p + b )! b !( a -p )! ( 2 r 2 f ) a -p ∆ τ 2 p -2 (cos γ -1) a -p . (26)  \nIn arriving at the last line of this expression, we made use of the identity  \nn ∑ l =0 ( -1) l (2 l +1) P l (cos γ ) ( n + l +1)!( n -l )! = (1 -cos γ ) n 2 n ( n !) 2 , (27)  \nwhich can be proved by induction. In the coincidence  \nlimit, the last line of (26) is nonzero only when p = a = 1.', 'B. Regularization Parameters for the Logarithmic Terms in the Tail': "In this subsection, we derive a mode-sum representation for the terms in the tail part of the Hadamard parametrix which in the partial coincidence limit ∆ r = 0 involves terms of the form ∆ τ 2 b s 2 a -2 b log( s 2 /ℓ 2 ). For later convenience, we will write the logarithmic terms as  \n∆ τ 2 b s 2 a -2 b log( s 2 /ℓ 2 ) = ∆ τ 2 b s 2 a -2 b log( ρ 2 ∆ τ 2 +1 -x ) +∆ τ 2 b s 2 a -2 b log(2 r 2 /ℓ 2 ) . (28)  \nThe second term here is polynomial in ∆ τ 2 and s 2 and requires no mode-sum representation. This term contains the arbitrary lengthscale ℓ which encodes the renormalization ambiguity. This term does not contribute to the vacuum polarization since a ≥ 1 in these terms and hence they all vanish in the coincidence limit. For the first term above, we assume the mode-sum ansatz  \n∆ τ 2 b s 2 a -2 b log( ρ 2 ∆ τ 2 +1 -cos γ ) = ∞ ∑ l =0 (2 l +1) P l (cos γ ) ∫ ∞ -∞ dω e iω ∆ τ χ ωl ( a, b | r ) (29)  \nand as before, we can invert to obtain a formal double integral expression for the regularization parameters  \nχ ωl ( a, b | r ) = (2 r 2 ) a -b 4 π ∫ 1 -1 dxP l ( x ) × lim α → a ∂ ∂α ∫ ∞ -∞ e -iωt t 2 b ( ρ 2 t 2 +1 -x ) α -b dt, (30)  \nwhere we have expressed the log terms as the limit of a derivative with respect to an exponent α , and where t and x are defined as before. The time integral above is a complicated generalized function but we can express it in terms of the derivative of a simpler integral as  \nχ ωl ( a, b | r ) = (2 r 2 ) a -b 4 π 1 2 b ( 1 ρ ∂ ∂ρ ) b ∫ 1 -1 dxP l ( x ) × lim α → a ∂ ∂α Γ( α -b +1) Γ( α +1) ∫ ∞ -∞ e -iωt ( ρ 2 t 2 +1 -x ) α dt. (31)  \nThe advantage of this form is that the time integral can now be performed explicitly yielding,  \nχ ωl ( a, b | r ) = (2 r 2 ) a -b √ 2 π 1 2 b ( 1 ρ ∂ ∂ρ ) b × ∫ 1 -1 dxP l ( x ) F ωl ( a, b | x, r ) (32)  \nwhere  \nF ωl ( a, b | x, r ) = lim α → a ∂ ∂α { 2 α ρ α -1 / 2 Γ( α -b +1) | ω | α +1 / 2 Γ( -α )Γ( α +1) (1 -x ) α/ 2+1 / 4 K α +1 / 2 ( | ω | √ 1 -x ρ )} , (33)  \nand K α +1 / 2 ( z ) is the modified Bessel function of the second kind. From the series representation of the Bessel function, we obtain  \nF ωl ( a, b | x, r ) = π √ 2 lim α → a ∂ ∂α { Γ( α -b +1) Γ( α +1) sec πα (2 ρ ) 2 α Γ( -α ) | ω | 2 α +1 ∞ ∑ k =0 1 k !Γ( -α + k + 1 2 ) ω 2 k (1 -x ) k (2 ρ ) 2 k -Γ( α -b +1) Γ( α +1) sec πα (1 -x ) α +1 / 2 Γ( -α )2 ρ ∞ ∑ k =0 1 k !Γ( α + k + 3 2 ) ω 2 k (1 -x ) k (2 ρ ) 2 k } . (34)  \nThe second term here is straightforward to deal with in the limit as α → a . In particular, it is evidently an entire function of α . In fact, if we consider the Taylor series for 1 / Γ( -α ) about a positive integer a , we have  \n1 Γ( -α ) =( -1) a +1 a !( α -a ) + ( -1) a +1 a ! ψ ( a +1)( α -a ) 2 + O ( α -a ) 3 , (35)  \nwhere ψ ( z ) = Γ ' ( z ) / Γ( z ) is the polygamma function, and so it is clear that the only term in the last line of (34)  \nthat survives the limit α → a is when the α -derivative acts on 1 / Γ( -α ).  \nThe first term in (34) is more subtle owing to the presence of the | ω | -2 α -1 term. Clearly this is not an ordinary function since it is singular at ω = 0. It is a generalized function whose definition involves a regularization at ω = 0 (not to be confused with the renormalization prescription required to regularize the coincidence limit of the point-split two-point function). One way to treat this is to simply introduce an infra-red cutoff in the ω -  \nral in (28), but this is not the rigorous approach. We will see that, if we treat | ω | -2 α -1 rigorously as a generalized function, a cutoff will emerge naturally and indeed it will be clear that the overall result will be independent of this choice. The salient point is that, as a function of the exponent α , the definition (which involves its regularization as a distribution) of | ω | -2 α -1 has a simple pole at integer values of α . However, this pole is canceled by the zero in 1 / Γ( -α ) at integer values (see Eq. (35)) so that the combination 1 / (Γ( -α ) | ω | 2 α +1 ) is actually an entire function. We require an expansion of this combination about the integer α = a up to order ( α -a ) 2 if we are to compute the limit appearing in (34). We already have a series for 1 / Γ( -α ). Next, we note the Laurent series for | ω | -2 α -1 which can be found in Ref. [26],  \n| ω | -2 α -1 = -δ (2 a ) ( ω ) (2 a )!( α -a ) + | ω | -2 a -1 + δ (2 a ) ( ω ) (2 a )! 2 log λ + O ( α -a ) , (36)  \nwhere the λ appearing here is an arbitrary inverse lengthscale required because ω has dimensions of inverse length. This inverse lengthscale will be shown later to be related to an infra-red cutoff in certain integrals over the frequency. The corresponding result for Eq. (36) in Ref. [26] is actually for a dimensionless generalized function and  \nso doesn't contain the additional log λ term, but applying the result from [26] to ˆ ω = ω/λ gives Eq. (36). It is essential to understand that | ω | -2 a -1 appearing on the right-hand side is to be understood in terms of its regularization as a generalized function. For example, for φ ( ω ) a test function of compact support, we have  \n( | ω | -2 a -1 , φ ( ω )) = ∫ ∞ 0 dω ω -2 a -1 [ φ ( ω ) + φ ( -ω ) -2 a -1 ∑ m =0 φ (2 m ) (0) (2 m )! ω 2 m -2 φ (2 a ) (0) (2 a )! ω 2 a θ ( λ -ω ) ] , (37)  \nwhere θ ( z ) is the step function. The last term in the regularization requires this step function since otherwise the integral would not converge at infinity.  \nNow combining (35) and (36) gives us the desired expansion  \n1 Γ( -α ) | ω | 2 α +1 = ( -1) a a ! (2 a )! δ (2 a ) ( ω ) +( α -a )( -1) a a ! [ ( ψ ( a +1) -2 log λ ) (2 a )! δ (2 a ) ( ω ) -1 | ω | 2 a +1 ] + O ( α -a ) 2 . (38)  \nEmploying these limits in Eq. (34) gives, after some tedious simplifications,  \nF ωl ( a, b | x, r ) = π √ 2 ( a -b )! { (( ψ ( a -b +1) + ψ ( 1 2 -a + k ) + 2 ln(2 ρ/λ )) a ∑ k =0 δ (2 a -2 k ) ( ω )(1 -x ) k (2 ρ ) 2 a -2 k k !(2 a -2 k )!Γ( k -a + 1 2 ) -∞ ∑ k =0 | ω | 2 k -2 a -1 (1 -x ) k (2 ρ ) 2 a -2 k k !Γ( k -a + 1 2 ) + ∞ ∑ k =0 ω 2 k (1 -x ) k + a +1 / 2 k !Γ( k + a + 3 2 )(2 ρ ) 2 k +1 } , (39)  \nwhere we also made use of the distributional identity  \nω 2 k δ (2 a ) ( ω ) = { (2 a )! (2 a -2 k )! δ (2 a -2 k ) ( ω ) k ≤ a 0 k > a (40)  \nto truncate the infinite sums involving delta distributions. Substituting this expression back into (32) and performing the integration yields  \nχ ωl ( a, b | r ) = √ π (2 r 2 ) a -b ( a -b )! ( -1) l 2 b ( 1 ρ ∂ ∂ρ ) b { ∞ ∑ k =0 Γ( k + a + 3 2 )2 k + a +3 / 2 ω 2 k k ! Γ( k + a + l + 5 2 ) Γ( k + a -l + 3 2 )(2 ρ ) 2 k +1 -∞ ∑ k =0 k !2 k (2 ρ ) 2 a -2 k | ω | 2 k -2 a -1 ( k + l +1)!( k -l )!Γ( k -a + 1 2 ) +[ ψ ( a -b +1) + ψ ( k -a + 1 2 ) + 2 log(2 ρ/λ )] a ∑ k =0 k !2 k (2 ρ ) 2 a -2 k δ (2 a -2 k ) ( ω ) (2 a -2 k )!( k + l +1)!( k -l )!Γ( k -a + 1 2 ) } . (41)  \nThe second term is a generalized function of ω whenever k ≤ a and an ordinary function otherwise. In the former case, it is understood in the sense of (37). It would be better if we could could express the definition in (37) in a way that is independent of the test function. It can be shown that the regularization in Eq. (37) is equivalent to defining  \nthe generalized function as  \n2 ∑ k =0 ( a -k )! k ! ( a -k -b )!Γ( k -a + 1 2 )( k + l +1)!( k -l )!(2 ρ 2 ) k ∑ p =0 p = a -k δ (2 p ) ( ω ) λ 2 p -2 a +2 k (2 p )!(2 p -2 a +2 k ) ,  \n̸  \n1 | ω | 2 j +1 = θ ( ω 2 -λ 2 ) | ω | 2 j +1 +2 ∞ ∑ k = j λ 2 k -2 j δ (2 k ) ( ω ) (2 k )!(2 k -2 j ) , j > 0 ∈ Z (42)  \n̸  \nwhere λ > 0 is the arbitrary cut-off. This can be straightforwardly proven by splitting the integral in (37) into intervals (0 , λ ) and ( λ, ∞ ), and then employing a Taylor series about ω = 0 for the test function on the interval (0 , λ ). Employing definition (42) in the second term of (41) for k ≤ a gives  \n̸  \nχ ωl ( a, b | r ) = √ π (2 r 2 ) a -b ( a -b )! ( -1) l 2 b ( 1 ρ ∂ ∂ρ ) b { ∞ ∑ k =0 Γ( k + a + 3 2 )2 k + a +3 / 2 ω 2 k k ! Γ( k + a + l + 5 2 ) Γ( k + a -l + 3 2 )(2 ρ ) 2 k +1 -θ ( ω 2 -λ 2 ) a ∑ k =0 k !2 k (2 ρ ) 2 a -2 k | ω | 2 k -2 a -1 ( k + l +1)!( k -l )!Γ( k -a + 1 2 ) -∞ ∑ k =0 ( k + a +1)!2 k + a +1 | ω | 2 k +1 ( k + a + l +2)!( k + a +1 -l )!Γ( k + 3 2 )(2 ρ ) 2 k +2 -2 a ∑ k =0 k !2 k (2 ρ ) 2 a -2 k ( k + l +1)!( k -l )!Γ( k -a + 1 2 ) ∞ ∑ p = a -k λ 2 p +2 k -2 a δ (2 p ) ( ω ) (2 p )!(2 p +2 k -2 a ) +[ ψ ( a -b +1) + ψ ( k -a + 1 2 ) + 2 log(2 ρ/λ )] a ∑ k =0 k !2 k (2 ρ ) 2 a -2 k δ (2 a -2 k ) ( ω ) (2 a -2 k )!( k + l +1)!( k -l )!Γ( k -a + 1 2 ) } .  \nFinally computing the derivatives and simplifying gives  \nχ  \nωl  \n(  \na, b  \n|  \nr  \nwhere  \nJ ωl ( a, b | r ) = √ π (2 r 2 ρ 2 ) a -b ( a -b )!2 2 a ( -1) l { -a -b ∑ k =0 θ ( ω 2 -λ 2 ) ( a -k )! k ! | ω | 2 k -2 a -1 ( a -k -b )!Γ( k -a + 1 2 )( k + l +1)!( k -l )!(2 ρ 2 ) k -( -1) b (2 ρ 2 ) a +1 | ω | ( a +1)! b ! 2 ˆ F 3 ( { a +2 , b +1 } { 3 2 , a -l +2 , a + l +3 } ω 2 2 ρ 2 ) + ( -1) b 2 a +1 / 2 ρ 2 a +1 Γ( a + 3 2 )Γ( b + 1 2 ) 2 ˆ F 3 ( { a + 3 2 , b + 1 2 } { 1 2 , a -l + 3 2 , a + l + 5 2 } ω 2 2 ρ 2 ) } . (45)  \nPutting this together gives  \n∞ ∞  \n∆ τ 2 b s 2 a -2 b log( ρ 2 ∆ τ +1 -cos γ ) = ∫ -∞ dω e iω ∆ τ ∑ l =0 (2 l +1) P l (cos γ ) χ ωl ( r ) = ∫ ∞ -∞ dω e iω ∆ τ ∞ ∑ l =0 (2 l +1) P l (cos γ ) J ωl ( a, b | r ) + √ π (2 r 2 ρ 2 ) a -b ( a -b )!2 2 a { a -b ∑ k =0 ( -1) a -k ( a -k )!∆ τ 2 a -2 k (1 -cos γ ) k 2 k k !( a -k -b )!(2 a -2 k )!Γ( k -a + 1 2 )(2 ρ 2 ) k [ [ ψ ( a -b +1) -ψ ( a -b -k +1) + ψ ( a -k +1) + ψ ( k -a + 1 2 ) + 2 log(2 ρ/λ ) ] -b ∑ k =1 ( -1) b ( b -k )!( k -1)!∆ τ 2 b -2 k (1 -cos γ ) k + a -b 2 k + a -b ( k + a -b )!(2 b -2 k )!Γ( k -b + 1 2 )(2 ρ 2 ) k + a -b -2 a -b ∑ k =0 ( a -k )!(1 -cos γ ) k 2 k k !( a -k -b )!Γ( k -a + 1 2 )(2 ρ 2 ) k m -k ∑ p =0 p = a -k ( -1) p ∆ τ 2 p λ 2 p -2 a +2 k (2 p )!(2 p -2 a +2 k ) } . (46)  \n) = J ωl ( a, b | r ) + √ π (2 r 2 ρ 2 ) a -b ( a -b )! 2 2 a ( -1) l × { a ∑ k =0 k !( a -k )! δ (2 a -2 k ) ( ω ) [ ψ ( a -b +1) -ψ ( a -b -k +1) + ψ ( a -k +1) + ψ ( k -a + 1 2 ) + 2 log( 2 ρ λ ) ] ( a -k -b )!(2 a -2 k )!Γ( k -a + 1 2 )( k + l +1)!( k -l )!(2 ρ 2 ) k -a ∞ } (44)  \n̸  \n(43)  \nEven though the sum over p in equation (46) is technically up to ∞ , we actually only need terms up to order ϵ 2 m in the Hadamard parametrix. In that expression, the order of the expansion is given by (1 -cos γ ) k ∆ τ 2 p and is ϵ 2 k +2 p . Hence, we have truncated the sum at p = m -k . In arriving at the expression above, we have also made  \nuse of the identity (27).", 'C. Mode-sum expression for the parametrix': "We now have mode-sum representations of all the singular terms in (9) yielding the following mode-sum representation of the Hadamard parametrix  \n(47)  \nK ( x, x ' ) = 1 8 π 2 ∫ ∞ -∞ dω e iω ∆ τ ∞ ∑ l =0 (2 l +1) P l (cos γ ) k ( m ) ωl ( r ) + 1 8 π 2 { m ∑ a =1 a ∑ b =1 D (p) ab ( r )∆ τ 2 a -2 b s 2 b -2 + m -1 ∑ a =1 a ∑ b =0 T (p) ab ( r )∆ τ 2 b s 2 a -2 b + m -1 ∑ a =0 a ∑ b =0 T (l) ab ( r )∆ τ 2 b s 2 a -2 b log(4 f/ ( ℓ 2 λ 2 )) + m ∑ a =1 a ∑ b =0 D (r) ab ( r ) (2 r 2 ) -b -1 b ! ( ρ 2 ) a + b +1 a ∑ p =1 ( -1) a + p ( a -p + b )! ρ 2 p ( a -p )! ∆ τ 2 p -2 (1 -cos γ ) a -p + m -1 ∑ a =1 a -1 ∑ b =0 T (r) ab ( r ) (2 r 2 ) -b -1 b ! ( ρ 2 ) a + b +2 a +1 ∑ p =1 ( -1) a +1+ p ( a +1 -p + b )! ρ 2 p ( a +1 -p )! ∆ τ 2 p -2 (1 -cos γ ) a +1 -p + m -1 ∑ a =0 a ∑ b =0 T (l) ab ( r ) √ π (2 r 2 ρ 2 ) a -b ( a -b )!2 2 a [ a -b ∑ k =0 ( -1) a -k ( a -k )!∆ τ 2 a -2 k (1 -cos γ ) k 2 k k !( a -k -b )!(2 a -2 k )!Γ( k -a + 1 2 )(2 ρ 2 ) k [ ψ ( a -b +1) -ψ ( a -b -k +1) + ψ ( a -k +1) + ψ ( k -a + 1 2 ) ] -b ∑ k =1 ( -1) b ( b -k )!( k -1)!∆ τ 2 b -2 k (1 -cos γ ) k + a -b 2 k + a -b ( k + a -b )!(2 b -2 k )!Γ( k -b + 1 2 )(2 ρ 2 ) k + a -b -2 a -b ∑ k =0 ( a -k )!(1 -cos γ ) k 2 k k !( a -k -b )!Γ( k -a + 1 2 )(2 ρ 2 ) k m -k ∑ p =0 p = a -k ( -1) p ∆ τ 2 p λ 2 p -2 a +2 k (2 p )!(2 p -2 a +2 k ) ]} + O ( ϵ 2 m log ϵ ) ,  \nwhere we have defined for compactness  \nk ( m ) ωl ( r ) = m ∑ a =0 a ∑ b =0 D (r) ab ( r ) I ωl ( a, b | r ) + m -1 ∑ a =1 a -1 ∑ b =0 T (r) ab ( r ) I ωl ( a +1 , b | r ) + m -1 ∑ a =0 a ∑ b =0 T (l) ab ( r ) J ωl ( a, b | r ) . (48)  \nWhile the expression (47) looks very complicated, we note that many of the terms vanish in the coincidence limit.", 'III. VACUUM POLARIZATION AND COMPONENTS OF THE RSET': "When computing vacuum polarization, the coincidence limit kills almost all of the finite terms in the Hadamard parametrix above, we obtain  \n⟨ ˆ ϕ 2 ⟩ ren = 1 8 π 2 { ∫ ∞ -∞ ∞ ∑ l =0 (2 l +1) ( g ωl ( r ) -k ( m ) ωl ( r ) ) dω -√ π m -1 ∑ a =1 a ∑ b =0 T (l) ab ( r ) ( a -1)!2 2 a f a -b Γ( -a + 1 2 ) λ 2 a } -1 8 π 2 [ D (p) 11 ( r ) + D (r) 10 ( r ) f + D (r) 11 ( r ) f 2 + T (l) 00 ( r )(log( f/ ( ℓ 2 λ 2 )) -2 γ E ) ] . (49)  \n̸  \nThe mode-sum terms above can be made to converge as fast as we like by choosing m appropriately large. It can be shown that the renormalized vacuum polarization is completely independent on the choice of the cutoff (we report the calculations in appendix A for clarity). This property can also be proven for the regularization terms for the stress energy tensor.  \nThe components of the RSET may be written in the form [17]:  \n⟨ ˆ T µ ν ⟩ ren = -w µ ν -( ξ -1 2 ) w ; µ ; ν +( ξ -1 4 ) □ w δ µ ν + ξR µ ν w -1 8 π 2 v 1 δ µ ν . (50)  \nwhere:  \nw ( r ) = ⟨ ˆ ϕ 2 ⟩ ren ( r ) , w µν ( x ) ≡ lim x → x ' [ W ( x, x ' ) ; µν ] , W ( x, x ' ) = G ( x, x ' ) -K ( x, x ' ) , (51)  \nand  \nv 1 = 1 720 R µνρσ R µνρσ -1 720 R µν R µν -1 24 ( ξ -1 5 ) □ R + 1 8 ( ξ -1 6 ) 2 R 2 + 1 4 µ 2 ( ξ -1 6 ) R + 1 8 µ 4 , (52)  \nwhich must be included in (50) to ensure the RSET is conserved [27].  \nDue to the symmetry of the background spacetime w is a function of r only, therefore once it has been calculated numerically to high accuracy on a suitably dense grid, the derivatives of w appearing in Eq. (50) are now easily and accurately obtained by differentiating an interpolation function for w . Considering the components of w a b , we have, by virtue of the wave equation satisified by W ( x, x ' ), that [28]:  \nw r r = -w τ τ -w θ θ -w ϕ ϕ -ξR w -µ 2 w -3 4 π 2 v 1 . (53)  \nFor the remaining non-zero components, w τ τ and w θ θ = w ϕ ϕ , it is advantageous to express these in terms of mixed derivatives at x and x ' of W ( x, x ' ) and derivatives of w ( r ) = ⟨ ˆ ϕ 2 ⟩ ren , using Synge's Rule [17]:  \nlim x → x ' [ W ( x, x ' ) ; µ ' ν ] = 1 2 w ; µν ( x ) -w µν ( x ) . (54)  \nThe required mixed time derivatives and mixed angular derivatives may be obtained by taking appropriate derivatives of Eq (9). Therefore, the RSET is determined by the calculation of three mode sums Eq.(49) along with  \nthe following two mode sums:  \n̸  \n[ g ττ ' W ,ττ ' ] = 1 8 π 2 f { ∫ ∞ -∞ ∞ ∑ l =0 (2 l +1) ω 2 ( g ωl -k ( m ) ωl ) dω -√ π m -1 ∑ a =0 a =1 a ∑ b =0 T (l) ab ( r ) a !2 2 a f a -b Γ( 1 2 -a )( a -1) λ 2 a -2 } + 1 4 π 2 [ T (p) 10 + D (p) 22 + 1 f ( T (p) 11 + D (p) 21 ) + D (r) 20 f 2 + D (r) 21 f 3 + D (r) 22 f 4 + T (r) 10 f 2 +( f T (l) 10 + T (l) 11 )(log( f/ ( ℓ 2 λ 2 )) -2 γ E +3) ] , (55)  \n[ g ϕϕ ' W ,ϕϕ ' ] = 1 16 π 2 r 2 { ∫ ∞ -∞ ∞ ∑ l =0 (2 l +1) l ( l +1) ( g ωl -k ( m ) ωl ) dω + √ π m -1 ∑ a =2 a ∑ b =0 T (l) ab r 2 ( a -b )( a -2)!2 2 a f a -b -1 Γ( 3 2 -a ) λ 2 a -2 } + 1 4 π 2 [ D (p) 22 + T (p) 10 -D (r) 20 f 2 -2 D (r) 21 f 3 -3 D (r) 22 f 4 -T (r) 10 f 2 + T (l) 11 f + r 2 T (l) 10 (log( f/ ( ℓ 2 λ 2 )) -2 γ E +1) ] . (56)  \nInserting the expressions in Eqs. (49), (55), (56) into Eq. (50) results in a natural splitting of the RSET components into a numeric and an analytic part:  \n⟨ ˆ T µ ν ⟩ ren = ⟨ ˆ T µ ν ⟩ numeric + ⟨ ˆ T µ ν ⟩ analytic , (57)  \nwhere ⟨ ˆ T µ ν ⟩ numeric contains the terms in Eqs. (49), (55), (56) that are in curly brackets, while ⟨ ˆ T µ ν ⟩ analytic is comprised of those terms in square brackets. We give explicitly the expression for ⟨ ˆ T µ ν ⟩ analytic in Appendix D. With these definitions, both ⟨ ˆ T µ ν ⟩ numeric and ⟨ ˆ T µ ν ⟩ analytic are independently conserved. In fact ⟨ ˆ T µ ν ⟩ analytic is equal, up to the addition of a conserved tensor, to the equivalent analytic component obtained by the AHS method [10], which in turn reproduces Page's approximation for a massless field. The fact that our direct method reproduces the AHS results, whose derivation depended on the use of the WKB approximation to the radial Green function, is surprising, particularly when contrasted with the results of the extended coordinates approach for a thermal field [17]. In that case the equivalent analytic component does not reproduce the AHS equivalent and neither the analytic nor numerical components are independently conserved, rather only their combination is conserved.", 'IV. APPLICATION TO THE REISSNER-NORDSTR OM SPACE-TIME': 'In this section, we apply the renormalization scheme presented developed herein to the case of scalar fields in the Boulware state in the Reissner-Nordstrom black hole space-time:  \nds 2 = -( r -r + ) ( r -r -) r 2 dt 2 + r 2 ( r -r + ) ( r -r -) dr 2 + r 2 ( dθ 2 +sin θ 2 dϕ 2 ) , (58)  \nwhere r ± = M ± √ M 2 -Q 2  \nWe will first verify this new method by reproducing results for the non extremal case, which were obtained via the state subtraction method in [18]. We then apply the method to the case of an extremal black hole, where we investigate the regularity of the Boulware state on the event horizon for both a massless and massive scalar field. While the regularity of the RSET has previously been examined by Anderson et al. [21], we are not aware of any previous exact results for a massive field.  \nBefore presenting results for both the non-extremal and extremal cases, we will discuss the numerical implementation of the method for both cases.', '1. Non-Extremal Modes': "As seen from Eqs. (49-56), the mode sums required to compute the RSET can be expressed in terms of derivatives of ⟨ ˆ ϕ 2 ⟩ ren . In particular, the only mode sum we need to calculate is that for the one-dimensional Green function g ωl (4), defined as the normalized product of solutions to the radial equation (5).  \nThus, to obtain the RSET with sufficient accuracy (as, for example, to check its regularity at the horizon) we need to generate highly-precise results for the vacuum polarization and, in turn, for the Euclidean modes.  \nFor the non-extremal case the computation of the p ωl ( r ) modes is simplified by recasting the radial equation, (5), into a confluent Heun form. This is achieved by expressing the dependent variable χ ωl ( r ) as:  \nχ ωl ( r ) = e -( ω -¯ ω ) r e ωr ∗ Y ( r )  \nand introducing a new independent variable:  \nz = r + -r r + -r -.  \nThe radial equation then takes the form of the confluent Heun equation:  \nz ( z -1) Y '' ( z ) + ( δ 1 ( z -1) + δ 2 z + z ( z -1) δ 3 ) Y ' ( z ) +( q 2 z -q 1 ) Y ( z ) = 0 (59)  \nwhere  \n(60)  \nq 1 = l ( l +1) + r 2 + ( ω -¯ ω ) 2 -( r + + r -)( ω -¯ ω ) -2 r + ¯ ω q 2 = ( r 2 + -r 2 -)( ω -¯ ω ) 2 -2( r + -r -)¯ ω δ 1 = 1 + 2 ωr 2 + r + -r -δ 2 = 1 -2 ωr 2 -r + -r -δ 3 = -2¯ ω ( r + -r -) z = r + -r r + -r -  \nand with  \n¯ ω = √ ω 2 + µ 2 . (61)  \nIf we let H ( q 1 , q 2 , δ 1 , δ 2 , δ 3 ; z ) denote the solution to Eq. (59) that is analytic in the vicinity of z = 0 and normalized to unity there, then we may express the p ωl solution as:  \np ωl ( r ) = e -( ω -¯ ω ) r e ωr ∗ H ( q 1 , q 2 , δ 1 , δ 2 , δ 3 ; z ) . (62)  \nComputing the q ωl ( r ) modes is computationally harder. While these modes can still be written in confluent Heun form, the Heun functions with the appropriate boundary conditions for q ωl ( r ) are not built into Mathematica or Maple. There are several options one can consider for computing q ωl ( r ) but we found it most efficient to simply numerically integrate the radial equation inwards from a large r value using an asymptotic expansion for the initial conditions (see for example [10] for details of the asymptotic expansion). The initial conditions were optimized so that the asymptotic expansion solved the wave equation to our working precision with the least number of terms in the asymptotic expansion and for the smallest reasonable r value at which this precision could be achieved. Using this approach, the mode solutions and their derivatives that were generated were accurate to at least 30 significant digits. We tested this accuracy by checking the constancy of the Wronskian over the radial grid for the solution pairs { p ωl ( r ) , q ωl ( r ) } .", '2. Extremal Modes': "The extremal black hole metric in its Euclideanized form is  \nds 2 = ( 1 -M r ) 2 dτ 2 + ( 1 -M r ) -2 dr 2 + r 2 ( dθ 2 +sin 2 θdϕ 2 ) , (63)  \nfor which the mode equation (5) amounts to  \nΦ '' ωl + 2Φ ' ωl r -M -[ r 4 ω 2 ( r -M ) 4 + µ 2 r 2 + l ( l +1) ( r -M ) 2 ] Φ ωl = 0 . (64)  \nFor convenience, we follow [29] and adopt the dimensionless coordinate η = r/M -1, for which the wave equation is  \nη 4 Φ '' ωl +2 η 3 Φ ' ωl -[ ω 2 M 2 ( η +1) 4 + l ( l +1) η 2 + µ 2 M 2 η 2 ( η +1) 2 ] Φ ωl = 0 . (65)  \nThe radial solutions regular at the horizon and at radial infinity obey, respectively, the asymptotic expansions  \np ωl = η 2 ωM e -ωM/η ∞ ∑ k =0 a k η k , η → 0 (66)  \nand  \nq ωl = η -2 ωM e ωM/η e -√ ω 2 + µ 2 Mη ∞ ∑ k =0 b k η -k -β 0 , η →∞ . (67)  \nReplacing these relations in (65) we find the following recursion relations for the a k , b k coefficients,  \na 1 = [ ˜ µ 2 + l ( l +1) + 2˜ ω (˜ ω -1) ] a 0 2˜ ω , a 2 = [ ˜ µ 2 +2˜ ω 2 ] a 0 2˜ ω + [ ˜ µ 2 + l ( l +1) + 2˜ ω (˜ ω -3) -2 ] a 1 4˜ ω , a k = ( ˜ µ 2 + ˜ ω 2 ) a k -3 2˜ ωu + ( ˜ µ 2 +2˜ ω 2 ) a k -2 ˜ ωu + [ k (1 -k ) + ˜ µ 2 + l ( l +1) -4˜ ωk +2˜ ω (1 + ˜ ω ) ] a k -1 2˜ ωk , (68)  \nfor the p ωl mode, where ˜ ω = ωM , ˜ µ = µM , with a 0 = 1. For the q ωl mode we have, instead,  \nβ 0 =1 -2˜ ω + √ ˜ ω 2 + ˜ µ 2 + ˜ ω 2 √ ˜ ω 2 + ˜ µ 2 , b 1 = [ ˜ µ 2 + l ( l +1) + β 0 (1 -β 0 -4˜ ω ) +2˜ ω ( 1 + ˜ ω -√ ˜ ω 2 + ˜ µ 2 )] b 0 × [ 2˜ µ 2 -2 β 0 √ ˜ ω 2 + ˜ µ 2 +4˜ ω ( ˜ ω -√ ˜ ω 2 + ˜ µ 2 )] -1 , b k = ( β 0 + v -2) ˜ ωb k -2 ( β 0 + k -1) √ ˜ ω 2 + ˜ µ 2 -2˜ ω ( ˜ ω -√ ˜ ω 2 + ˜ µ 2 ) -˜ µ 2 [ -k 2 + k (3 -2 β 0 -4˜ ω ) + l ( l +1) + ˜ µ 2 -β 2 0 + β 0 (3 -4˜ ω ) + 2˜ ω ( 3 + ˜ ω -√ ˜ ω 2 + ˜ µ 2 ) -2 ] b k -1 × [ 2 ( β 0 + k -1) √ ˜ ω 2 + ˜ µ 2 -4˜ ω ( ˜ ω -√ ˜ ω 2 + ˜ µ 2 ) -2˜ µ 2 ] -1 , (69)  \nwith b 0 = 1.", '3. Accuracy Issues': 'The numerical implementation of the direct method is more challenging than the extended coordinate method developed for thermal states, which uses a discretized frequency spectrum.  \nMost of the issues are related to the sampling of the continuous frequency spectrum. To numerically perform the integral in ω , one needs to interpolate the integrand over a predefined ω -grid. This grid must be fine enough so that the interpolation function works well, ensuring the result is independent of the grid choice. This becomes challenging when the integrand is very large and rapidly decaying. In this context, we found it necessary to push the cutoff λ to values no smaller than 1 / 10 and to increase the density of the ω -grid immediately after the cutoff. This is done to ensure the interpolation function can properly approximate the rapidly decaying behavior of the cutoff dependent terms in the integrand, which are much larger at smaller ω .  \nAnother issue is related to the cancellation of the cutoffdependent terms. We have shown analytically that the renormalized vacuum polarization does not depend on the choice of the cutoff, as the cutoff-dependent terms ultimately cancel each other. However, achieving this outcome numerically is more challenging and can significantly impact the accuracy of the results. For the case of a thermal state, one can increase the accuracy of the results, by increasing the expansion order of the singular field ( m in Eq. (9)). Moreover, this is extremely efficient numerically as only a very modest number of modes are required. However for the direct approach, as the magnitude of the cutoff-dependent terms increases with increasing order, so does the potential error introduced by the numerical cancellation of these terms. To counteract this, as the expansion order increases, one has to also increase the density in the grid in ω immediately after the cutoff. In effect if one increases the expansion order, one has to also significantly increase the number of numerical modes required, making the process increasingly inefficient. Hence, there is a balance to be struck between the accuracy of the results and the efficiency of the method. For the results contained in this paper, we have found that a third order expansion ( m = 3) and summing/integrating up to l = 40/ ω = 5 gives sufficiently accurate results, with the conservation equation satisfied to approximately 9-10 decimal places.  \nThere is an important exception to this choice and that is in the vicinity of the event horizon. All the cutoffdependent terms vanish on the event horizon and therefore as we approach the horizon the error introduced by these terms also decreases. This is particularly true in the extremal case, where the metric function f ( r ) has a double root. Therefore, for points close to the horizon (up at a distance of approximately M/ 10) we employ a sixth order expansion to obtain our results and find that the conservation equation is satisfied to approximately 13-15 decimal places for a massive field and 10-11 deci-  \nmal places for the massive case. In addition, particularly near the horizon, the integral is dominated by the low ω behavior. Hence, it is crucial to select the smallest point in the ω grid accordingly and extrapolate the integrand to ω = 0 to avoid losing important contributions.  \nBesides the problems related to the continuous ω spectrum, we have found different representations for I ωl ( a, b | r ) and J ωl ( a, b | r ). The most amenable to numerical implementation are the ones reported in the derivations above, which involve differences of two hypergeometric functions. However, these representations fail whenever z = | ω | r √ f ≫ 1, due to catastrophic cancellation. For the region of the domain where z ≫ 1, the representations in terms of derivatives of modified Bessel functions and generalized incomplete gamma functions are more efficient and not subject to numerical errors. We include them in appendix B for convenience.', 'B. Non-Extremal Results': 'In Fig. (1) we compare the results of the direct method outlined in this paper to the corresponding results obtained via a state-subtraction approach in [18], for the calculation of the RSET for a massless, minimally coupled scalar field in the Boulware state. The background spacetime on which this calculation is performed is the non-extremal Reissner-Nordstrom spacetime with Q = 2 M/ 10. As is evident from the plot in Fig. (1), there is good agreement between the two approaches, verifying the new method presented in this paper. However, the state-subtraction method appears to be more accurate, as measured by the conservation equation, in the vicinity of the event horizon. This suggests that for space-times where there exists a Hartle-Hawking state to use as a reference state, taking the state-subtraction approach outlined in [18] is the appropriate choice. However, for other space-times, such as the extremal Reissner-Nordstrom spacetime considered in the next section, where we do not have a Hartle-Hawking state to exploit, the direct approach of this paper is required for the RSET calculation for scalar field in the Boulware state.  \nReturning to the splitting into independently conserved tensors in (57), we observe that the leading-order divergence of the RSET at the horizon is determined entirely by ⟨ ˆ T µ ν ⟩ analytic ∝ f -2 , with ⟨ ˆ T µ ν ⟩ numeric ∝ f -1 amounting to a sub-leading correction. This indicates that ⟨ ˆ T µ ν ⟩ analytic indeed captures exactly the singular behaviour of the Boulware state at event horizons, and further justifies its use in backreaction analyses as an analytic approximation to the exact RSET [30].', 'C. Extremal Results': "In this section we present the results of the extended coordinates approach to the calculation of the RSET for a  \nFIG. 1: RSET components in the Boulware state for a massless, minimally coupled scalar field calculated via the direct method (line) and the state subtraction method (dots). The background spacetime is the non-extremal Reissner-Nordstrom spacetime with Q = 2 M/ 10.  \n<!-- image -->  \nscalar field in the Boulware quantum state on an extremal Reissner-Nordstrom spacetime.  \nIn Figs. (2) and (3) we present the components of the RSET in the exterior region of the black hole for a massless and massive field ( µM = 1 / 10) respectively. For the massless case, on the horizon ⟨ ˆ T a b ⟩ numeric vanishes and the RSET is given precisely by ⟨ ˆ T a b ⟩ analytic with  \n⟨ ˆ T µ ν ⟩| r = M = ⟨ ˆ T µ ν ⟩ analytic | r = M = 1 2880 π 2 M 4 δ µ ν , (70)  \nwhich is equal to the corresponding RSET in the Bertotti-Robinson spacetime [21]. The fact that ⟨ ˆ T µ ν ⟩ analytic is exact on the horizon, suggests that it might serve as a good approximation for the full RSET, at least in the massless case. However this is not the case, with ⟨ ˆ T µ ν ⟩ analytic failing to capture the main features of the full results, except in the immediate vicinity of the horizon. Finally, as r increases we see that all RSET components tend to 0, as would be expected for a zero-temperature state.  \nFor the massive case, the situation is quite different. Here both ⟨ ˆ T µ ν ⟩ analytic and ⟨ ˆ T µ ν ⟩ numeric diverge logarithmically on the horizon, however these divergences cancel to render each of the components of ⟨ ˆ T µ ν ⟩ finite there. Specifically, the analytic terms become, in the near-horizon limit,  \n⟨ ˆ T µ ν ⟩ analytic | r = M = µ 2 8 π 2 [ ( ξ -1 6 ) M 2 diag (1 , 1 , -1 , -1) + µ 2 4 δ µ ν ] log ( r M -1 ) . (71)  \nIn contrast to the massless case, only the ⟨ ˆ T r r ⟩ and ⟨ ˆ T t t ⟩ components agree on the horizon with ⟨ ˆ T θ θ ⟩ taking on a distinct value. For massive fields, the r →∞ behaviour of the RSET components depends on the renormalisation lengthscale ℓ . Hence, there is a natural choice of the  \nFIG. 2: Plot of the components ⟨ ˆ T µ ν ⟩ ren for a conformally coupled massless field in the Boulware state, on an extremal Reissner-Nordstrom spacetime.  \n<!-- image -->  \nlengthscale for massive fields for which the RSET tends to 0 as r →∞ . This choice is given by:  \nℓ = 1 2 µ exp { 3 / 4 + 5 γ E } ,  \nand is the value chosen for the results presented in Figs. (3) and (5). For the massless case, there is no such natural choice, so we simply set ℓ = M .  \nIn Figs. (4) and (5) we focus on the near horizon region and investigate the regularity of the Boulware state in an extremal Reissner-Nordstrom spacetime for a massless and massive field ( µM = 1 / 10) respectively. The quantity ⟨ ˆ T µ ν ⟩ analytic is not regular on the event horizon resulting in an infinite energy density there as perceived by a freely falling observer:  \n⟨ ˆ E⟩ analytic = ⟨ ˆ T r r ⟩ analytic -⟨ ˆ T t t ⟩ analytic f = -µ 2 [ 3 M 2 µ 2 +12 ( ξ -1 6 )] 96 π 2 M 2 ( r M -1 ) -2 + 15 M 4 µ 4 -2 240 π 2 M 4 ( r M -1 ) -1 -1 60 π 2 M 4 log ( r M -1 ) + O ( r M -1 ) 0 . (72)  \nThis quantity possesses a ( r -M ) -1 divergence for a massless field and a ( r -M ) -2 for a massive field. In Figs. (4) and (5), we show that this divergence appears to be precisely canceled by a corresponding divergence from ⟨ ˆ T µ ν ⟩ numeric , leaving a total ⟨ ˆ E⟩ that is finite on the event horizon. Of course, to prove this definitively would require a closed form expression for the quantity ⟨ ˆ E⟩ on the horizon (or equivalently a near horizon expression for ⟨ ˆ T µ ν ⟩ to O ( r -M ) 2 ). However, we feel the numerical results presented here provide compelling evidence for the regularity of the Boulware state.  \nWe have therefore used the novel method presented in the paper to confirm the results of [21] for a massless field  \nFIG. 3: Plot of the components ⟨ ˆ T µ ν ⟩ ren for a minimally coupled field with mass µM = 1 / 10 in the Boulware state, on an extremal Reissner-Nordstrom spacetime.  \n<!-- image -->  \nFIG. 4: Plot of the cancellation of the divergence in the numerical and analytical contributions to the quantity ⟨ ˆ E⟩ for a minimally coupled massless field, leaving a total result that appears to be finite on the event horizon of the extremal Reissner-Nordstrom spacetime  \n<!-- image -->  \n-  \nand provide, to the best of the authors' knowledge, the first exact results for the RSET of a massive scalar field in the extremal Reissner-Nordstrom spacetime.", 'V. BACKREACTION IN EXTREMAL BLACK HOLES': "Computing the RSET in the background of an extremal black hole naturally begs the question of whether extremal horizons are stable under quantum corrections. Extremal black holes, despite being in a static zerotemperature vacuum state and thus not suffering from any evaporation process, generate a RSET which is (generically) non-vanishing at the horizon and will backreact on the background metric. In this section we assume the RSET is perfectly regular and take this effect into consideration by expanding the Einstein equations  \nFIG. 5: Plot of the cancellation of the divergence in the numerical and analytical contributions to the quantity ⟨ ˆ E⟩ for a minimally coupled field with mass µ = M/ 10, leaving a total result that appears to be finite on the event horizon of the extremal Reissner-Nordstrom spacetime  \n<!-- image -->  \n-  \nperturbatively in ℏ ,  \nG µ (0) ν + ℏ G µ (1) ν + O ( ℏ ) 2 = 8 π ( T µ (em) ν + ℏ ⟨ ˆ T µ ν ⟩ ) (73)  \nwhere G µ (0) ν is the Einstein tensor evaluated in the background spacetime and T µ (em) ν is the stress-energy tensor of the electromagnetic field. Since we are considering a real quantum scalar field, the electromagnetic stressenergy tensor is not affected by semiclassical contributions [21]. Following [5] we express the metric components as  \ng tt = -e 2 ψ ( r ) (1 -2 m ( r ) /r ) , g rr = (1 -2 m ( r ) /r ) -1 , (74)  \nwhere m ( r ) denotes the Misner-Sharp mass, and expand the metric functions as  \nψ ( r ) ≃ log [1 + ϵ R ( r )] , m ≃ M ( 1 -M 2 r ) [1 + ϵ M ( r )] , (75)  \nwith ϵ = ℏ /M 2 . To linear order in ϵ , we write the semiclassical Einstein equations in the form  \n2 ( z -1) 2 R ' ( z ) M 2 z 3 = 8 πM 2 ( ⟨ ˆ T r r ⟩ - ⟨ ˆ T t t ⟩ ) , -M ( z ) M 2 z 4 -(2 z -1) M ' ( z ) M 2 z 3 = 8 πM 2 ⟨ ˆ T t t ⟩ , (76)  \nwhere we have used the dimensionless coordinate z = r/M . The first equation is clearly the difference between the radial and temporal equations, and the angular equations are such that the Bianchi identities are satisfied. Once the RSET is specified (here we replace it by the exact, numerical values for massless and massive fields shown in the previous section), we can integrate Eqs. (76) for M and R . The full solutions for the perturbations are numeric and we can only calculate them  \nwithin the region for which we have numerical RSET values at our disposal.  \nFor massless fields in the near-horizon limit z → z e , we can make use of the fact that the energy density and pressures of the exact RSET are given, to leading order in z -z e , by (70) to obtain the following analytic solutions for the perturbations  \nM = -C 0 1080 π + O ( z -z e ) , R = K 0 1080 π + O ( z -z e ) , (77)  \nwhere { C 0 , K 0 } are arbitrary integration constants. For massive fields, however, we cannot follow a similar argument since we ignore the exact values the RSET must reach at the extremal horizon, and extrapolating our numerical results introduces spurious errors. Written in ingoing Eddington-Finkelstein coordinates, the metric components become  \ng vv = g tt = -ϵC 0 1080 π + O ( z -z e ) , g vr = e ψ = 1 + ϵK 0 1080 π + O ( z -z e ) . (78)  \nIn light of these expressions, we observe that the metric receives a correction so that g tt is no longer vanishing at z = z e . For negative C 0 we have g tt > 0, hence the surface z = z e appears to lie inside a trapped region. In this case perturbative corrections have shifted the position of the horizon from its classical value z e = 1 to a new position z H > z e which must be found numerically. The surface gravity associated to this quantum-corrected horizon is equal to  \nκ H = e ψ 2 M ∂ ∂z ( 1 -2 Mm z ) + e ψ M ( 1 -2 Mm z ) ∂ψ ∂z ∣ ∣ ∣ ∣ z H . (79)  \nExtremal horizons are characterized by a vanishing surface gravity, i.e., a quadratic root in g vv . In the Reissner-Nordstrom metric, this comes from the fact that Q = M . The RSET makes the Misner-Sharp mass acquire a radial dependence, but does not affect the charge, allowing for the possibility of breaking the delicate balance of scales required by extremal horizons. In practice what we observe is that the backreacted spacetime, locally around z = z e , becomes like the sub-extremal or super-extremal Reissner-Nordstrom spacetimes depending which sign we choose for the mass perturbations.  \nIn particular, together with shifting the position of the horizon outwards, taking ϵC 0 < 0 makes it acquire a positive surface gravity. The relative position and surface gravity of the quantum-corrected horizon are shown in Fig. 6. For massless fields, we see this non-extremal horizon smoothly tends to a horizon of the extremal kind in the C 0 → 0 limit. In the massive case, since we lack RSET values arbitrarily close to z e , we have not been able to generate reliable results in the C 0 → 0 limit. Nonetheless, for large and negative C 0 values, the curves for the massless and massive cases are indistinguishable.  \nFIG. 6: Position of the perturbed horizon ∆ z = z H -z e and its surface gravity κ H for a semiclassically corrected extremal black hole. Corrections from massive and massless fields shift the Misner-Sharp mass by a constant value, which either breaks the extremal horizon into a pair of inner and outer horizons (if C 0 < 0), or eliminates any horizon entirely (if C 0 > 0). Our results do not change qualitatively as long as ϵK 0 ≪ 1, so we set K 0 = 0 here.  \n<!-- image -->  \nSince we do not have exact results for z < z e , we cannot obtain the quantum-corrected metric in the region z ≪ z H . However, it is clear (since the corrections are perturbative) that we are observing a splitting of the extremal horizon into a pair of outer and inner horizons. By solving (73) to the next order in ℏ , we would find: i) that the vacuum state in which this RSET is evaluated is no longer a zero-temperature state, and ii) that this state is singular at the newly formed inner horizon [31]. For positive C 0 the function g vv has zero roots, hence the extremal horizon disappears entirely and the perturbed metric describes a horizonless (and singular) object, whose associated vacuum state has zero temperature. Finally, between both regimes there is a separatrix solution with C 0 = 0 that preserves both the position and the surface gravity of the extremal horizon. This solution could potentially exist for all orders in ℏ only if we fine-tune the perturbations so that their contribution to the Misner-Sharp mass vanishes at every order in ϵ . Of course, the treatment adopted in this section only deals with static solutions from the start, hence we remain agnostic in regard to the stability of extremal horizons under time-dependent perturbations.", 'VI. DISCUSSION AND CONCLUSIONS': "In this paper, we have provided a generalisation of the extended coordinate prescription to directly compute the renormalized stress-energy tensor (RSET) of scalar fields in the Boulware vacuum. While the extended coordinates approach was first developed for the calculation of expectation values for a scalar field in the Hartle-Hawking state, it can be used to obtain results for other quantum states using a state subtraction approach [18]. However,  \nfor spacetimes where one does not have a Hartle-Hawking state to leverage as a reference state, such as an extremal black hole or the spacetime of a star then a direct approach to compute the RSET is required.  \nIn the extended coordinate approach we take an expansion of the Hadamard parametrix in extended coordinates, which encodes the short-distance behaviour of the two-point function, and we rewrite the terms in this expansion as mode-sums of the same form as those in the unrenormalized Green function, which may then be renormalized in a mode-by-mode manner. The main technical challenge in applying this approach to a scalar field at zero temperature arises from having a continuous frequency spectrum and hence the mode sums involve integrals over ω (as opposed to the Hartle-Hawking case where we sums over discrete values of ω ). As result, when we express terms in the Hadamard parametrix as mode-sums, we find that the integrands are not ordinary functions of ω and we must treat these mode-sums rigorously as generalized functions. When we do this, we see that an infra-red cutoff naturally appears in the integrals over ω in such a way so that the final results are independent of this cut-off.  \nHaving developed this approach, we first applied it to the non-extremal Reissner-Nordstrom spactime. In this case we reproduced the results obtained via the state subtraction approach in [18], thereby verifying the new direct approach. We then applied the direct approach the the calculation of the RSET in an extremal ReissnerNordstrom spacetime, for both a massless and massive scalar field. The massless case had previously been considered in [21], however, to the best of the authors' knowledge, the results in this paper are first results obtained for a massive field.  \nOnce we have calculated the RSET components we used these results to investigate the regularity of the Boulware state by calculating the energy density on the horizon as perceived by a freely falling observer. Here we find strong numerical evidence that the Boulware state is regular for both the massless and massive case. However, we stress that these results are indicative only and to say anything conclusive requires exact results on the horizon. We hope to report on this in the near future.  \nDespite lacking a definite proof of the RSET regularity at extremal horizons, we can nonetheless assume that it is indeed regular to give a hint on the backreaction effects it would entail. By solving the static semiclassical equations perturbatively in ℏ [5] we have seen that, unless perturbations are fine-tuned to every order in the expansion, extremal horizons are unstable. If these perturbations contribute positively to the Misner-Sharp mass, the extremal horizon splits into a pair of outer and inner horizons, while in the opposite scenario the horizon disappears completely, leaving no marginally trapped surfaces in the spacetime.  \nThe generalization of the extended coordinate method to the Boulware state opens the window towards obtaining the RSET of scalar fields in previously unexplored sit-  \nuations, particularly, in stellar spacetimes. Recent works that considered perfect fluid stars show indications that the RSET becomes a dominant contribution in stellar interiors as these stars approach their maximum compactness limits [32-34]. Since previous works were based on different RSET approximations, a calculation of the exact RSET will allow to address their validity. These explorations have great implications for the theoretical plausibility of ultracompact stars sourced by semiclassical effects.", 'ACKNOWLEDGMENTS': 'The authors thank the Institute for Fundamental Physics of the Universe for hosting their visit under the Team Research program Dark Compact Objects and The Quantum Vacuum . JA acknowledges funding from the Italian Ministry of Education and Scientific Research (MIUR) under the grant PRIN MIUR 2017-MB8AEZ. The authors would like to thank Marc Casals for helpful conversations.', 'Appendix A: Cut-Off Independence': 'There is a subtle cancellation happening between the terms dependent on λ in J ωl ( a, b | r ) and the finite terms in (47) which also depend on λ . We can make this cancellation explicit by analyzing the terms in the renormalized vacuum polarization that contain the cutoff, which we refer to as w ( λ ) :  \nw ( λ ) = T (l) 00 ( r ) 8 π 2 [ ∫ ∞ -∞ dω θ ( ω 2 -λ 2 ) | ω | -1 +log( λ 2 ) ] + √ π 8 π 2 m -1 ∑ a =1 a ∑ b =0 T (l) ab ( r ) f a -b 2 2 a Γ( -a + 1 2 ) [ a ! ∫ ∞ -∞ dω θ ( ω 2 -λ 2 ) | ω | -2 a -1 -( a -1)! λ 2 a ] . (A1)  \nWe have separated the T ( l ) 00 ( r ) and a > 0 terms in J ωl ( a, b | r ) and used the identity (27) to truncate the l sum as well as introduced k = 0, only condition for these terms being nonzero at coincidence. We can show that the renormalized vacuum polarization is independent on the choice of the cutoff by verifying that  \n∂ ∂λ ( w ( λ ) ) = T (l) 00 ( r ) 8 π 2 ( -2 λ + 2 λ ) + √ π 8 π 2 m -1 ∑ a =1 a ∑ b =0 T (l) ab ( r ) f a -b 2 2 a Γ( -a + 1 2 ) ( -2 a ! λ -2 a -1 +2 a ! λ -2 a -1 ) = 0 , (A2)  \nwhere we made use of  \n∂ ∂λ ( θ ( ω 2 -λ 2 )) = -( δ ( ω + λ ) + δ ( ω -λ ) ) . (A3)  \nSimilarly, it can be demonstrated that the RSET is also completely independent of the choice of the cutoff.', 'Appendix B: Other representations for direct and tail regularization parameters': "In this appendix, we report different but equivalent representations for the functions I ωl ( a, b | r ) and J ωl ( a, b | r ). Above, we have presented the derivation for these expressions in terms of hypergeometric functions. The alternative representations below involve the modified Bessel functions I ν ( z ) and K ν ( z ) and the generalized incomplete gamma  \nfunction Γ( a, z 0 , z 1 ) = ∫ z 1 z 0 t a -1 e -t dt .  \nI ωl ( a, b | r ) = r 2 b +1 f -b -1 / 2 (2 b b !) (2 r 2 ) b +1 b ∑ p =0 ( -1) a 2 p ( p ) 2 b -2 p Γ( b -p +1) ∂ 2 a + p ∂ω 2 a + p ( ω p I l + 1 2 ( ωr √ f ) K l + 1 2 ( ωr √ f )) = 1 2 ( -1) a 2 2 b b !( f ) b +1 / 2 r b ∑ p =0 2 p ( p ) 2 b -2 p Γ( b -p +1) 2 a + p ∑ k =0 ( 2 a + p k ) p ! ω -2 a + k ( -2 a + k )! k ∑ d =0 ( k d )( r √ f ) k ( -1) d 2 k × 2 k -2 d ∑ c =0 ( k -d c ) I l +1 / 2 -k + d +2 c ( ωr √ f ) 2 d ∑ n =0 ( d n ) K l +1 / 2 -d +2 n ( ωr √ f ) = lim x →-1 2 l ( -1) a + b ( 2 r 2 ) a -1 2 2 b ! f a + b + 1 2 { b ∑ k =1 l ∑ j =0 ( -1) j +1 j ! ( l j )( 1 2 ( l + j -1) l ) 2( -k + j +1)! ∂ b -k ∂x b -k ( (1 -x ) a + b -1 2 e -ω √ 2(1 -x ) r 2 f ) + l ∑ k = b k -b ∑ j =0 k ! ( l k )( 1 2 ( k + l -1) l ) ( -1) -2 b + k -j ( k -b )! ( k -b j ) ( ω √ 2 r 2 f ) -2 a -2 k +2 j -1 Γ ( 2 a +2 k -2 j +1 , 0 , 2 rω √ f ) } , (B1)  \n(B2)  \nJ ωl ( a, b | r ) = -√ π (2 r 2 ρ 2 ) a -b ( a -b )!2 2 a ( -1) l a -b ∑ k =0 θ ( ω 2 -λ 2 ) ( a -k )! k ! | ω | 2 k -2 a -1 ( a -k -b )!Γ( k -a + 1 2 )( k + l +1)!( k -l )!(2 ρ 2 ) k +(2 r 2 ) a -b ( -1) b Γ(1 + a -b ) r f 1 / 2 a -b +1 ∑ k =0 ( -1) k 2 a -b +2 ( a -b +1 k ) (2 l +2( a -b ) -4 k +3) ( l -k + 1 2 ) a -b +2 × ∂ 2 b ∂ω 2 b I l + 3 2 + a -b -2 k ( ωr √ f ) K l + 3 2 + a -b -2 k ( ωr √ f ) = -√ π (2 r 2 ρ 2 ) a -b ( a -b )!2 2 a ( -1) l a -b ∑ k =0 θ ( ω 2 -λ 2 ) ( a -k )! k ! | ω | 2 k -2 a -1 ( a -k -b )!Γ( k -a + 1 2 )( k + l +1)!( k -l )!(2 ρ 2 ) k +(2 r 2 ) i -j ( -1) j ( r √ f ) 2 j +1 Γ(1 + i -j ) i -j +1 ∑ k =0 ( -1) k 2 i -j +2 ( i -j +1 k ) (2 l +2( i -j ) -4 k +3) ( l -k + 1 2 ) i -j +2 × 2 j ∑ p =0 ( 2 j p ) ( -1) p 2 2 j 4 j -2 p ∑ q =0 ( 2 j -p q ) I l +3 / 2+ i -3 j -2 k + p +2 q ( ωr √ f ) 2 p ∑ n =0 ( p n ) K l +3 / 2+ i -j -2 k -p +2 n ( ωr √ f ) = -lim ω ' → ω f a -b 4 π a -b ∑ k =0 2 l +2 ( -1) a -k ( a -b k ) ∂ 2 a -2 k ∂ω ' 2 a -2 k l ∑ n =0 n ∑ j =0 π ω ' ( -1) n -j ( l n )( 1 2 ( l + n -1) l )( n j ) × ( 2 r 2 f ) k    Γ ( 2( k -j + n +1) , 0 , 2 ω ' r √ f ) ( √ 2 ω ' r √ f ) 2( k -j + n +1) -2 k -j + n θ ( λ 2 -ω 2 ) k -j + n +1    .", 'Appendix C: Hadamard Coefficients': "Below we list the coefficients D (p) ij , D (r) ij , T (l) ij , T (r) ij and T (p) ij to 2nd order.  \nD (r) 00 = 2 , D (r) 10 = -f ( r 2 f '' -2 rf ' +2 f -2 ) 12 r 2 ,  \nD (r) 11 = f ( r 2 ( f ' 2 ) -4 f ( rf ' +1) + 4 f 2 ) 2  \n24 r , D (r) 20 = 1 2880 r 4 ( f [ -5 r 2 ( -f ' 2 ) ( r 2 f '' -2 rf ' -2 ) -8 f 2 ( 3 r 3 f (3) -7 r 2 f '' +19 rf ' +10 ) + f ( 9 r 4 f '' 2 -20 r 2 f '' +86 r 2 f ' 2 +4 rf ' ( 3 r 3 f (3) -14 r 2 f '' +20 ) +4 ) +76 f 3 ] ) , D (r) 21 = -1 2880 r 4 ( f [ r 4 ( f ' 2 + f ' 4 ) + r 2 f ( -30 rf ' 3 + f ' 2 ( 11 r 2 f '' -30 )) +4 f 3 ( 11 r 2 f '' -52 rf ' -40 ) -2 f 2 ( 10 r 2 f '' -67 r 2 f ' 2 + f ' ( 22 r 3 f '' -80 r ) -28 ) +104 f 4 ] ) , D (r) 22 = f 2 ( r 2 f ' 2 -4 f ( rf ' +1) + 4 f 2 ) 2 1152 r 4 , D (p) 11 = -f ' 6 r , D (p) 21 = f ( -9 rf ' 2 +6 rf ( rf (3) -2 f '' ) +2 f ' ( 7 r 2 f '' + f +5 )) 720 r 3 , D (p) 22 = 7 r 2 f ' 2 -10 rf ' + rf (9 rf '' +4 f ' ) -3 f 2 +3 720 r 4 , T (l) 00 = 6 µ 2 r 2 -(6 ξ -1) ( r 2 f '' +4 rf ' +2 f -2 ) 12 r 2 , T (l) 10 = 1 480 r 4 ( r 2 f '' ( (10 ξ (3 ξ -1) + 1) r 2 f '' -10(6 ξ -1) ( 2 ξ + µ 2 r 2 )) +4(5 ξ (24 ξ -5) + 1) r 2 f ' 2 +2 rf ' ( (5 ξ -1) r 3 f (3) +2 ( 60 ξ 2 -5 ξ -1 ) r 2 f '' +40(1 -6 ξ ) ξ +10 µ 2 (1 -12 ξ ) r 2 ) +2 rf ( (80 ξ (3 ξ -1) + 6) f ' + r ( (20 ξ (3 ξ +2) -8) f '' + r ( (5 ξ -1) rf (4) ( r ) + (40 ξ -7) f (3) ))) -40 fξ ( 6 ξ +3 µ 2 r 2 -1 ) +4(10 ξ (3 ξ -1) + 1) f 2 +2 ( 15 ( 2 ξ + µ 2 r 2 ) 2 -2 ) ) T (l) 11 = f 480 r 4 ( -2 ( r 2 f '' +2 ) ( (1 -5 ξ ) r 2 f '' +10 ξ +5 µ 2 r 2 -2 ) +2(1 -10 ξ ) r 2 f ' 2 + rf ( (16 -80 ξ ) f ' + r ( (12 -80 ξ ) f '' + r ( rf (4) ( r ) + (6 -20 ξ ) f (3) ) +20 µ 2 )) +8(5 ξ -1) f 2 + rf ' ( (10 ξ -1) r 3 f (3) +4(20 ξ -3) r 2 f '' +20(6 ξ -1) ) ) T (r) 10 = f ( r 2 f ' 2 -4 f ( rf ' +1) + 4 f 2 ) ( (1 -6 ξ ) (2 -r ( rf '' +4 f ' ) -2 f +2) -6 µ 2 r 2 ) 576 r 4 , T (p) 10 = ( f -1) ( (1 -6 ξ ) (2 -r ( rf '' +4 f ' ) -2 f ) -6 µ 2 r 2 ) 144 r 4 , T (p) 11 = f ( rf ' -2 f +2) ( (1 -6 ξ ) (2 -r ( rf '' +4 f ' ) -2 f ) -6 µ 2 r 2 ) 144 r 4 .  \n,", 'Appendix D: Analytical Part of the RSET': "In this appendix, we give explicit expressions for ⟨ ˆ T µ ν ⟩ analytical that arises in the split in Eq. (57) of the RSET into a numerical and analytical part. The analytical part is separately conserved and given by  \n⟨ ˆ T µ ν ⟩ analytical = -¯ w µ ν -( ξ -1 2 )¯ w ; µ ν +( ξ -1 2 ) □ ¯ w + ξ R µ ν ¯ w -1 8 π 2 v 1 , (D1)  \nwhere  \n¯ w = -1 8 π 2 { D (r) 10 f + D (r) 22 f 2 + D (p) 11 + T (l) 00 ( r )(log( f/ ( ℓ 2 λ 2 )) -2 γ E ) } (D2)  \n¯ w τ τ = 1 4 f ' ¯ w ' -1 4 π 2 { T (r) 10 f 2 + T (p) 10 + T (p) 11 f + T (r) 10 f 2 + D (p) 22 + D (p) 21 f + D (r) 20 f 2 + D (r) 21 f 3 + D (r) 22 f 4  \n+(  \nf T (l) 10 + T (l) 11 )(log( f/ ( ℓ 2 λ 2 )) -2 γ E +3) } (D3)  \n¯ w ϕ ϕ = ¯ w θ θ = 1 2 r f ¯ w ' -1 4 π 2 { T (p) 10 + T (l) 11 f -T (r) 10 f 2 + D (p) 22 -D (r) 20 f 2 -2 D (r) 21 f 3 -3 D (r) 22 f 4  \n+ r 2 T (l) 10 (log( f/ ( ℓ 2 λ 2 )) -2 γ E +1) } (D4)  \n¯ w r r = -¯ w τ τ -2¯ w ϕ ϕ -ξ R ¯ w -µ 2 ¯ w -3 4 π 2 v 1 ,  \nwhere explicit expressions for the Hadamard coefficients in terms of f are found in the previous appendix.  \n- [1] P. Anderson, Phys. Rev. D 28 , 271 (1983).\n- [2] L. Parker and J. Z. Simon, Phys. Rev. D 47 , 1339 (1993), arXiv:gr-qc/9211002.\n- [3] R. Parentani and T. Piran, Phys. Rev. Lett. 73 , 2805 (1994), arXiv:hep-th/9405007 [hep-th].\n- [4] L. Gao, P. R. Anderson, and R. S. Link, Phys. Rev. D 108 , 105011 (2023), arXiv:2308.11040 [gr-qc].\n- [5] J. W. York, Phys. Rev. D 31 , 775 (1985).\n- [6] E. E. Flanagan and R. M. Wald, Phys. Rev. D 54 , 6233 (1996).\n- [7] P. Taylor and C. Breen, Phys. Rev. D 103 , 025006 (2021), arXiv:2007.09405 [gr-qc].\n- [8] K. W. Howard and P. Candelas, Phys. Rev. Lett. 53 , 403 (1984).\n- [9] P. R. Anderson, W. A. Hiscock, and D. A. Samuel, Phys. Rev. Lett. 70 , 1739 (1993).\n- [10] P. R. Anderson, W. A. Hiscock, and D. A. Samuel, Phys. Rev. D 51 , 4337 (1995).\n- [11] C. Breen and A. C. Ottewill, Phys. Rev. D 85 , 084029 (2012).\n- [12] A. Levi and A. Ori, Phys. Rev. D 91 , 104028 (2015).\n- [13] A. Levi and A. Ori, Phys. Rev. Lett. 117 , 231101 (2016).\n- [14] A. Levi, E. Eilon, A. Ori, and M. van de Meent, Phys. Rev. Lett. 118 , 141102 (2017).\n- [15] P. Taylor and C. Breen, Phys. Rev. D 94 , 125024 (2016).\n- [16] P. Taylor and C. Breen, Phys. Rev. D 96 , 105020 (2017).\n- [17] P. Taylor, C. Breen, and A. Ottewill, Phys. Rev. D 106 ,  \n065023 (2022).  \n- [18] J. Arrechea, C. Breen, A. Ottewill, and P. Taylor, Phys. Rev. D 108 , 125004 (2023).\n- [19] J. B. Hartle and S. W. Hawking, Phys. Rev. D 13 , 2188 (1976).\n- [20] D. G. Boulware, Phys. Rev. D 11 , 1404 (1975).\n- [21] P. R. Anderson, W. A. Hiscock, and D. J. Loranz, Phys. Rev. Lett. 74 , 4365 (1995), arXiv:gr-qc/9504019.\n- [22] P. R. Anderson, W. A. Hiscock, and B. E. Taylor, Phys. Rev. Lett. 85 , 2438 (2000).\n- [23] D. A. Lowe, Phys. Rev. Lett. 87 , 029001 (2001).\n- [24] P. R. Anderson, W. A. Hiscock, and B. E. Taylor, Phys. Rev. Lett. 87 , 029002 (2001).\n- [25] DLMF, ' NIST Digital Library of Mathematical Functions ,' https://dlmf.nist.gov/ , Release 1.2.1 of 202406-15, f. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, B. V. Saunders, H. S. Cohl, and M. A. McClain, eds.\n- [26] I. M. Gel'fand and G. E. Shilov, ' Generalized Functions, V. 1: Properties and Operations ,' (Academic Press, 1964).\n- [27] R. M. Wald, ' Quantum field theory in curved spacetime and black hole thermodynamics ,' (The University of Chicago Press, 1994).\n- [28] M. R. Brown and A. C. Ottewill, Phys. Rev. D 34 , 1776 (1986).\n- [29] P. R. Anderson, Phys. Rev. D 41 , 1152 (1990).  \n(D5)  \n- [30] J. Arrechea, C. Barcel'o, R. Carballo-Rubio, and L. J. Garay, Phys. Rev. D 107 , 085005 (2023), arXiv:2212.09375 [gr-qc].\n- [31] N. Zilberman, M. Casals, A. Ori, and A. C. Ottewill, Phys. Rev. Lett. 129 , 261102 (2022), arXiv:2203.08502 [gr-qc].  \n- [32] W. A. Hiscock, Phys. Rev. D 37 , 2142 (1988).\n- [33] I. A. Reyes and G. M. Tomaselli, Phys. Rev. D 108 , 065006 (2023), arXiv:2301.00826 [gr-qc].\n- [34] J. Arrechea, C. Barcel'o, R. Carballo-Rubio, and L. J. Garay, Phys. Rev. D 109 , 104056 (2024), arXiv:2310.12668 [gr-qc]."}