Renormalized stress-energy tensor for scalar fields in Hartle-Hawking, Boulware and Unruh states in the Reissner-Nordstr\"om spacetime

In this paper, we consider a quantum scalar field propagating on the Reissner-Nordstr\"om black hole spacetime. We compute the renormalized stress-energy tensor for the field in the Hartle-Hawking, Boulware and Unruh states. When the field is in the Hartle-Hawking state, we renormalize using the recently developed ``extended coordinate'' prescription. This method, which relies on Euclidean techniques, is very fast and accurate. Once, we have renormalized in the Hartle-Hawking state, we compute the stress-energy tensor in the Boulware and Unruh states by leveraging the fact that the difference between stress-energy tensors in different quantum states is already finite. We consider a range of coupling constants and masses for the field and a range of electric charge values for the black hole, including near-extreme values. Lastly, we compare these results with the analytic approximations available in the literature.


I. INTRODUCTION
The renormalized expectation value of the quantum stress-energy tensor plays a crucial role in the semiclassical theory of gravity.It governs the quantum backreaction on the classical spacetime geometry via the semiclassical field equations where g ab is the metric of spacetime, G ab is the Einstein tensor and Tab R is the renormalized expectation value of the stress-energy tensor of a quantum field in some quantum state.The tensors H ab , H ab are geometrical terms that are quadratic in the curvature and arise through the point-splitting regularization process that yields Tab R .This regularization process corresponds to an infinite renormalization of the constants Λ, α and β.
The calculation of renormalized expectation values of the stress-energy tensor (RSET) for a particular quantum state is a technically challenging endeavour.For black hole spacetimes, there are three main approaches to this calculation, the most established being the Candelas-Howard approach [1] and its extensions (see for example [2][3][4]).Recently two new, more efficient methods for the calculation of the RSET have been developed, the first being the "pragmatic mode-sum prescription" [5,6].The method has proven indeed to be pragmatic, both in its efficiency and its broader applicability.The second recent development in methods to compute the RSET in black hole spacetimes is known as the "extended coordinate method" [7][8][9], which is extremely efficient and applicable to arbitrary field parameters and arbitrary spacetime dimensions.
Of the three standard quantum states considered in spherically symmetric black hole space-times, the Hartle-Hawking state [10,11] has received the most attention (see for example [1][2][3][4][12][13][14][15]).This is mainly due to the fact that this state can be defined on the Euclidean metric and there are some convenient simplifications that occur in this Euclidean setting, for example, the frequency spectrum is discrete.When it comes to the other two quantum states of interest, the Boulware and Unruh states [16,17], the literature is more sparse.Anderson, Hiscock and Samuel developed a method for the calculation of the RSET of a scalar field with arbitrary mass and coupling in the Boulware state in a general spherically symmetric (Euclidean) metric [2].Jensen, Mc Laughlin and Ottewill calculated the RSET for a massless, conformally coupled field in the Unruh and Boulware states in the Schwarzschild spacetime [18,19].More recently, the pragmatic mode-sum prescription has been applied arXiv:2307.10307v1[gr-qc] 18 Jul 2023 to the calculation of the RSET for the Boulware and Unruh states for a massless minimally coupled field in the Schwarzschild, Reissner-Nordström and Kerr spacetimes [5,6,20].Anderson, Siahmazgi, Clark and Fabbri also applied the pragmatic mode-sum prescription to the calculation of the RSET for a massless minimally coupled field in a spacetime where a black hole forms from the collapse of a null shell [21].The relative lack of RSET results for the Unruh state and the fact that there does not appear to be any results whatsoever in the literature for a massive field in this state is surprising, given that the Unruh state is considered to be the one of most physical interest, that is, the state that models the late-time evolution of the collapse of a spherical body to a black hole [17].
The extended coordinate renormalization method [7][8][9] mentioned above was developed using Euclidean techniques and is applicable only to a quantum field in the Hartle-Hawking state.While it is possible, by considering a field at zero temperature in the Euclidean metric, to extend the method to the Boulware vacuum, it is more efficient to adopt a state subtraction approach using the Hartle-Hawking as our reference state.In the state subtraction approach one leverages the fact that the difference between stress tensors in different quantum states does not require renormalization, therefore by calculating the RSET in the Hartle-Hawking state we may then obtain the RSET in the Unruh and Boulware states without recourse to renormalization.We note that for the case of a spacetime where we can not define a Hartle-Hawking vacuum to use as our reference state (such as in Kerr black holes [22]), we must then consider the direct calculation for the Boulware state through the zero temperature extended coordinate method mentioned above.We hope to present results on this in the near future.
In this paper we employ the extended coordinate method to compute the RSET for a scalar field in the Hartle Hawking state propagating in a general spherically symmetric black hole spacetime.We assume the field has arbitrary mass and coupling to the background curvature.We then show how to compute the RSET in the Boulware and Unruh states by a state subtraction scheme, using the RSET in the Hartle-Hawking state as a reference state.We then apply these results to the particular case of the Reissner-Nordström spacetime and calculate the RSET in all three quantum states for a range of quantum field masses and coupling constants.We also consider a range of electric charge values for the black hole and probe the near extreme case.Finally we compare the exact RSET results with various analytic approximations and judge their reliability.

II. RENORMALIZATION PRESCRIPTION IN THE HARTLE-HAWKING STATE
In this section, we will briefly outline the extended coordinate approach to calculating the RSET for a quan-tum scalar field in the Hartle-Hawking quantum state, propagating on a static, spherically symmetric black hole spacetime.It is convenient to construct the Hartle-Hawking state on the Euclideanized line element (2) where in order to avoid a conical singularity on the event horizon r = r + , it is necessary to impose the periodicity τ = τ + 2π/κ + on the Euclidean time, where is the horizon surface gravity.Imposing this periodicity discretizes the frequency spectrum of the field modes which now satisfy an elliptic equation where E is the d'Alembertian operator with respect to the Euclidean metric, µ is the field mass, 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, where ∆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 nl (r) := κ + p nl (r) q nl (r)/N nl the one-dimensional radial Green function evaluated at the same spacetime point r.The radial modes p nl (r), q nl (r) are solutions of the homogeneous radial equation: where p nl (r) and q nl (r) are regular on the horizon and the outer boundary (usually spatial infinity), respectively.The normalization constant is given by where W{p, q} denotes the Wronskian of the two solutions.
In the coincidence limit ∆x → 0 (i.e.γ → 0 and ∆τ → 0), the mode sum (5) diverges.To renormalize this mode sum, we must find a way to express the locallyconstructed Hadamard parametrix K(x, x ) as a mode sum and subtract mode-by-mode.In [7,8] a mode sum expression for the Hadamard parametrix was derived by first introducing the so-called extended coordinates: For simplicity, the separation in the radial direction, ∆r, is set to zero but it is important to the development that the separation in the other directions is maintained.Expressing the Hadamard parametrix in terms of these extended coordinates permits its decomposition in terms of Fourier frequency modes and multipole moments where, remarkably, the coefficients in this decomposition are expressible in closed form for any static sphericallysymmetric spacetime in arbitrary dimensions.In four dimensions, the result is where Here m denotes the order of the expansion, the coefficients ij (r) and T (r) ij (r) arise in the expansion of the Hadamard parametrix K(x, x ) in the extended coordinates s and while the terms Ψ (+) nl (i, j|r) and χ nl (i, j|r) are the so-called regularization parameters that arise in expressing K(x, x ) as a mode sum.The well-known renormalization ambiguity is expressed as an arbitrary lengthscale in the regularization parameters χ nl (i, j|r).We find that all the regularization parameters are obtainable in closed form in terms of complicated combinations of special functions.Explicit expressions for each of these are given in [9].
To apply the extended coordinate method to the calculation of the RSET, we first introduce some notation: were we have adopted square brackets [•] to denote the coincidence limit x → x.Then the RSET may be writ-ten in the form [9]: where: We note that since φ2 R is a function of r only, once it has been calculated numerically to high accuracy on a suitaby dense grid, derivatives of φ2 R are easily and accurately obtained by differentiating an interpolation function for φ2 R .Considering next the components of w a b , we have, by virtue of the wave equation satsified by W (x, x ), that [23]: For 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 R , using Synge's Rule [9]: The required mixed time derivatives and mixed angular derivatives may then in turn be expressed as mode sums, given by [9]: where g nl (r) ≡ g nl (r)−k nl (r).Therefore, the application of the extended coordinates method reduces the calculation of the RSET to that of three mode sums, those just presented along with The modes g nl (r) converge like O(l −2m−3 ) for large l, fixed n and O(n −2m−3 ) for large n, fixed l, where we remind the reader that m is the order of the expansion of the singular field.We therefore have precise control over the convergence of the mode sums by choosing the order of the expansion.By choosing a sufficiently high order Hadamard parametrix, very high accuracy in the RSET calculation can be achieved by truncating the sums at a modest number of l and n modes.

III. RSET IN THE BOULWARE AND UNRUH STATES
In this section, we will outline our strategy for computing the RSET in the Boulware and the Unruh states.There are two approaches one might consider.The first is to adapt the extended coordinate approach to decompose the Hadamard parametrix for the Lorentzian spacetime and then apply the Lorentzian equivalent of the mode sum renormalization prescription described in the previous section.While taking this direct approach has some advantages, there are subtle difficulties in dealing with the Hadamard distribution on the Lorentzian sector.The other approach, which is the one we adopt here, is to use the fact that the differences between states does not require renormalization.This is because the Hadamard singularity structure is agnostic to the quantum state.So even though the Euclideanization trick employed herein is relevant only to the Hartle-Hawking state, we can still use this state as a reference to compute the RSET in other quantum states of interest.
To be more explicit, let us start by writing out the Wightman Green function in the Boulware, Hartle-Hawking and Unruh states in terms of their Lorentzian modes.Mode solutions to the Klein-Gordon equation for a massive scalar field have the form where Y lm (θ, φ) are the spherical harmonics and Φ ωl (r) is a solution of the radial equation It is helpful to recast this equation in Schrodinger form by writing Φ ωl (r) = ψ ωl (r)/r, where ψ ωl (r) satisfies where dr * /dr = 1/f is the radial tortoise coordinate and the potential is From here on, we will work in asymptotically flat spacetimes, so spacetimes with a non-zero cosmological constant are excluded from the analysis.Importantly, this potential asymptotes to different values at the horizon and infinity for non-zero field mass.We have Hence the solutions to Eq. ( 21) have the asymptotic forms ψ ωl ∼ e ±iωr * as r → r + and ψ ωl ∼ e ±iωr * as r → ∞.On the exterior, we take as a linearly independent basis the solutions with the following boundary conditions.We label the solution ψ in ωl (r) as the one with boundary condition e −iωr * on past null-infinity and which vanishes on the past event horizon.Note that for |ω| < µ, this would be an exponentially growing mode and the solutions would not be square integrable.Hence for the "in"-mode, we must have the restriction ω > µ.This wave which originates at past null infinity partly reflects back to future null infinity and partly transmits to the future event horizon.In terms of the boundary conditions only on the radial function, this amounts to where B in ωl and A in ωl are respectively the dimensionsless transmission and reflection coefficients for the "in" mode.For the other independent solution, we label ψ up ωl (r) as the one with asymptotic form e iωr * on the past event horizon and which vanishes on past null infinity.This solution represents a wave propagating out of the event horizon and being partly scattered back to the future event horizon and partly transmitted to future null infinity.Then we have the asymptotic forms where B up ωl and A up ωl are again transmission and reflection coefficients, respectively.Note that for the "up" modes with |ω| < µ, the solution is oscillatory near the horizon but exponentially damped at infinity.Hence these modes are square integrable and must be included in the twopoint function.The "up" modes with |ω| < µ are the bound-state modes [16].
Taking as the basis modes to the wave operator then the normalization conditions imply where the inner product is defined by the following integral on an arbitrary Cauchy surface Σ with unit futuredirected normal n a , In deriving the normalization conditions, we move the integral over Σ to an integral over the past event horizon H − plus an integral over past null infinity I − , making use of the asymptotic forms ( 23)- (24).We note the different normalizations for the "in" modes and "up" modes in (25).Moreover, the constancy of the Wronskian of linearly independent solutions of the radial equation implies The details of how to numerically compute the normal-ized modes {Φ in ωl (r), Φ up ωl (r)} are briefly discussed in the next section.
Putting these details together, and performing the trivial m-sum in the spherical harmonics, we get the following representation of the two-point function in terms of the normalized modes for the field in the Boulware state [16], For the field in the Hartle-Hawking state [10] there is a thermal factor coth(ω π/κ + ) in each of the "in" modes and "up" modes.We obtain, Finally, for the field in the Unruh state [17], we pick up a thermal factor on the "up" modes but not the "in" modes: Now if we consider a quantum scalar field in any reference Hadamard state |R , then the RSET in any other Hadamard state |Q is related to the RSET in our reference state by where is the difference between our two-point function in our state |Q and our reference state.The explicit dependence on the field mass, i.e., the µ 2 term appearing in the expression for the stress-energy tensor (12), canceled on application of the wave equation to δG Q ; though obviously the above expression depends on µ implicitly through T a b R and through δG Q itself.The salient point is that δG Q is a smooth homogeneous solution of the wave equation, at least on any region of the spacetime where both states satisfy the Hadamard condition.
In the current context, the reference state is the Hartle-Hawking state for which we can leverage Euclidean techniques to compute T a b HH , then in the Boulware and Unruh states, we have As we can see from the exponential thermal factor in each of these expressions, they are exponentially convergent in ω and reasonably straightforward to compute.
For convenience, we give the explicit mode-sum expres-sions for the RSET in the Boulware and Unruh states in terms of the RSET in the Hartle-Hawking state.For the field in the Boulware state, we have with and where, in order to arrive at these expressions, we have made use of Synge's Rule and the radial wave equation to avoid computing second derivatives of our radial modes at the same spacetime point.
To obtain the analogous expressions for the quantum field in the Unruh state, we simply omit the contributions from the "up" modes in each component in (37).The RSET in the Unruh state also has an off-diagonal component Tt which gives the outgoing flux of Hawking radiation.
It is straightforward to verify that in each state |Q , one obtains for the trace of the conformally-coupled stress-energy tensor which yields the well-known trace anomaly for massless fields.

IV. NUMERICAL IMPLEMENTATION
The main numerical task in the implementation of the prescriptions outlined in Sections II and III is to obtain the radial modes by numerically integrating the radial equation for the quantum state under consideration.The methods for treating the radial differential equation are quite different for the Euclidean modes and the Lorentzian modes, so we discuss each case separately.

A. Euclidean Modes
The two-point function on the Euclidean slice is obtained by the usual separation of variables procedure, the problem essentially reducing to a one-dimensional radial Green function g nl (r, r ), where the mode numbers n and l are discrete.Solving this radial Green function amounts to computing a normalized product of solutions to the homogeneous equation (6).The Hartle-Hawking state can be uniquely defined on this Euclidean slice by the condition of regularity on the event horizon.We denote the radial solution regular on the event horizon to be p nl (r) while the solution q nl (r) is the solution regular at infinity but divergent at the event horizon.To impose regularity of the Green function at the horizon, we must ensure that p nl (r) is evaluated at the smaller of the two radial points, i.e., we take g nl (r, r ) = κ + p nl (r < )q nl (r > )/N nl , where r < ≡ min{r, r }, r > = max{r, r }, and the factor of κ + is incorporated for convenience into our definition of g nl (r, r ) so that it satisfies an equation with a κ + δ(r − r ) source term.The normalization constant comes from the Wronskian N nl = −r 2 f (r)W{p nl , q nl } All that remains is to compute the solutions p nl (r) and q nl (r).The computation of the former is simplified by recasting the radial equation into a confluent Heun form [9]. Since the confuent Heun functions that are regular at a regular singular point are built into many software suites such as Mathematica or Maple, computing these presents no difficulty.In particular, if we let H(q, α, γ, δ, ; z) be the confluent Heun function that solves that is analytic in the vicinity of z = 0 and normalized to unity at z = 0, then it is straightforward to show that: where and where in the expressions above we have used the notation Computing the q nl (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 nl (r) are not built into Mathematica or Maple.There are several options one can consider for computing q nl (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.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 nl (r), q nl (r)}.
We wanted to compute the RSET to a high degree of accuracy which ostensibly requires a large set of {n, l} modes.However, when employing the extended coordinates prescription as outlined in Section II, the number of modes required can be significantly reduced by taking a suitably high order expansion of the singular field.Here we choose to take a 6th order expansion (setting m = 6) in Eq. ( 10) and generate 40 l modes and 15 n modes, which yields the RSET accurate to approximately 10-15 decimal places for the parameter sets considered in this paper.

B. Lorentzian Modes
In this section, we briefly describe the computation of the normalized Lorentzian modes {Φ in ωl (r), Φ up ωl (r)}.The boundary conditions on these modes are expressed in terms of transmission and reflection coefficients via ( 23)-( 24) with ( 27)- (29).All that remains to compute is the relationship between the reflection/transmission coefficients and a pair of convenient numerically computed radial modes { Φin ωl (r), Φup ωl (r)}.Starting first with the "in" modes.As in the Euclidean discussion above, the homogeneous solutions of the radial equation can be expressed in terms of confluent Heun functions, we take as Φin ωl (r) the solution where here and the parameters in the confluent Heun function are obtained from (44) by applying the transformations ω → −iω, ω → −iω.Comparing with the asymptotic forms, we see that As in the Euclidean case, we solve for Φup ωl (r) numerically by integrating the radial equation inwards from a suitably large radius with our initial conditions determined by an asymptotic series.The leading order term in this asymptotic series is and the numerical solutions are related to the normalized modes by With these definitions, the Wronskian condition implies where W ωl is the constant associated with the Wronskian of the numerical radial modes: The frequency integrals in the the RSET components have a rapid convergence so we compute modes on a frequency grid out to ω = 10.The error induced by truncating at this upper bound is very small.We choose a grid that is finely meshed near the lower limit of integration.For the l sums appearing in (37), the convergence with l is very rapid also, except for sums over the "up" modes very close to the horizon.Again because we desired to have very accurate results, we computed 150 l modes for each frequency.The error induced by truncating the sums at l = 150 is tiny for all but the points closest to the horizon.For example, for a massive field in the Reissner Nordström spacetime with Q/M = 0.2 and µ M = 0.1 (with Q and M denoting the electromagnetic charge and the ADM mass, respectively), for a fixed frequency, the error induced by truncating the l sums at l = 150 is less than 38 significant figures for points (r − r + )/M 0.05.
The last numerical issue we wish to briefly mention is that it is necessary to compute Lorentzian bound state modes with ω < µ.These are more difficult to compute than the modes with ω > µ.The main problem is that the asymptotic expansion does not converge to a sufficiently accurate value unless r is chosen to be quite large whence the value of the initial conditions is usually very small.This forces one to increase the working precision of the numerical integrator and hence slows down the computation.The problem is accentuated for larger field masses.Fortunately, large field masses are not the physically relevant cases.

V. RESULTS
We have used the method previously outlined to obtain the RSET of massless and massive scalar fields in the Reissner-Nordström spacetime where The RSET is an ambiguous quantity due to its dependence on the arbitrary lengthscale present in the 4-dimensional Hadamard parametrix.Similar to the Hadamard parametrix, such renormalization ambiguities are local and independent of the vacuum state under consideration.Consequently, all dependence on is already contained in T a b HH .In spherical symmetry, Anderson, Hiscock and Samuel [2] (AHS) showed that -dependent terms amount to a covariantly conserved analytic stress-tensor.Particularizing for the Reissner-Nordström spacetime (53), this contribution takes the form with where clearly, for µ = 0, the terms (55) vanish at the event horizons for Q = 0 and everywhere for Q = 0.The inherent ambiguities in the definition of the RSET can make comparisons between the results of different approaches to its calculation difficult.However, we are able to write down a relationship between results obtained via the AHS approach, which is based on Christensen's DeWitt-Schwinger expansion [24], and the extended coordinate method, which is based on the Hadarmard parametrix.To do this we note that in addition to the ambiguity in the choice of lengthscale mentioned above, Mc Laughlin proved that, for a scalar field, an RSET obtained via Christensen's or by the Hadamard approach will differ by a geometric term, given by [25]: This is clearly a conserved quantity that can be absorbed into the semiclassical field equations via a renormalization of the constants G and Λ. Therefore we have that the RSETs obtained via the two approaches are related by where γ E is Euler's constant and ˆ = EC / AHS is the ratio of the choice of lengthscale made in both approaches.
As the AHS approach is based on the DeWitt-Schwinger expansion, AHS is conventionally taken to be equal to 1/µ for a massive field and is arbitrary otherwise.Eq. (57) then enables meaningful comparisons between the two approaches, in particular in section VI, it will allow us to compare our exact results with the AHS and De-Witt Schwinger approximations to the RSET.
As the space of parameters to explore is large, in the next subsections we consider first how the different vacuum states affect the RSET in the massless and massive cases, to later analyse the impact of varying the coupling ξ and the charge Q.We will compare these exact results with the values predicted by analytic RSET approximations in section VI.

A. Vacuum states
Figure 1 contains plots of the RSET components multiplied by f 2 for {ξ = 0, Q/M = 0.2} and {ξ = 1/6, Q/M = 0.2}, both with µM = 0 (continuous lines) and µM = 0.1 (dashed lines), in the first and second rows, respectively.The third row depicts the nearextremal case {ξ = 1/6, Q/M = 0.99} with µM = 0.For every coupling, charge, and mass, all RSET components in the Hartle-Hawking states are finite at the event horizon (strictly, they are finite at the lowest point in our radial grid, which can be taken as close to r = r + as desired, however the regularity of all components on the horizon was proven in [26]).As expected, every RSET component in the Boulware state diverges at the horizon in a way ∝ f −2 for all the couplings considered, while for the Unruh state T r r and T t t diverge like ∝ f −1 .For large r, the diagonal RSET components in the Unruh and Boulware states approach the same value, in accordance with the asymptotic behaviours found in [27,28].We also find that in a freely falling frame, as expected, T a b is regular on the future event horizon for the Unruh state, regular on the future and past event horizons for the Hartle-Hawking state, and diverges on both the past and future horizons in the Boulware state.
Increasing the mass of the field affects the value and sign of the RSET at large r.In the Boulware and Unruh states, these do not decay to zero asymptotically, as with µM = 0, but approach a constant value instead.These asymptotic values can be identified with the choice of arbitrary length scale.In fact, from the numerical results we find that the RSET for the Boulware and Unruh states approach the following value From inspection of Eq (57), we see that this corresponds to the large r expansion of the difference between the extended coordinates and AHS RSETs (with AHS = 1/µ).
For massive fields, we therefore have a natural choice for the lengthscale EC , for which the Boulware and Unruh RSET decays to 0 as r → ∞, given by We make this choice for the massive field results presented in this Section.For the massless field, the large r asymptotic behaviour is independent of the choice of lengthscale, so for simplicity, in this case we set EC = M , where M is the black hole mass.The aforementioned properties hold for all Q < M .As Q increases, the finite terms relating the various states in Eq. ( 37) decrease in magnitude due to their dependence on the black hole temperature, and the Boulware, Unruh and Hartle-Hawking states converge.The Unruh and Hartle-Hawking states converge faster than the Boulware state, since contributions from the "up" modes to Eq. ( 37) decrease slower in the Q → M limit than contributions from the "in" modes.Results for the extremal case Q = M will appear elsewhere.

B. Varying the coupling
Next we consider the effect of varying the field coupling ξ. Figure 2 shows the diagonal RSET components in the Hartle-Hawking, Boulware and Unruh states (first, second and third row, respectively) for Q/M = 0.6.Components in the Boulware and Unruh states have been multiplied by powers of f for visualization purposes.
Increasing the value of the coupling does not affect notably the magnitude of the RSET, but it does produce a change in the sign of every RSET component with the exception of T t t U and T r r U .As was shown in [9] for the Hartle-Hawking state, RSET values for any coupling can be generated from results in minimal coupling (this follows from definition (12) in spacetimes with vanishing Ricci scalar).With this we find that, for Q/M = 0.6, the T t t HH and T r r HH components change sign at the lowest grid point at ξ ≈ 0.24 and the T θ θ HH component component does at ξ ≈ 0.19.For the first two components, said sign changes occurs at smaller ξ as Q increases, whereas the converse happens to the latter component.For the Boulware state, the analytic RSET approximation from [2] suggests that the sign of the RSET at the horizon is independent from Q (see Sec. VI below).
For diagonal RSETs, the point-wise null energy condition is satisfied as long as the inequalities hold everywhere outside the event horizon [29,30].The Hartle-Hawking state satisfies the null energy condition everywhere outside the event horizon for ξ = {0, 1/8, 1/6}, while for ξ = 1/2 it is violated from r ≈ 2.8r + to r = r + .The Boulware state violates the null energy condition everywhere outside and at the event horizon for ξ = {0, 1/8, 1/6} and satisfies it everywhere for ξ = 1/2.Note that these behaviours depend on the particular choice of renormalization lengthscale (54) and hence should not be considered to be physically meaningful statements.

C. Varying the charge
Varying the black hole charge Q affects both the sign and the magnitude of the RSET components, which become more sensitive to slight variations of Q as the extremal limit Q = M is approached.Figure 3 shows the RSET components for ξ = 0 and the charge values Q/M = {0.4,0.8, 0.9, 0.95, 0.99}.We study the minimally coupled case in detail to allow for a faithful comparison with the analytic approximations presented in Section VI.
For the RSET in the Hartle-Hawking state, increasing Q decreases the asymptotic values in every RSET component at large r.Asymptotically, this state reproduces a thermal bath at equilibrium with the horizon temperature, which decreases with increasing Q.At the horizon r = r + , the components T t t HH and T r r HH increase with Q except for Q/M = 0.99, for which these components are smaller than in the Q/M = 0.95 case.For T φ φ HH , this inversion happens at lower Q and is already noticed in the Q/M = 0.95 case, even becoming negative at the event horizon for Q/M = 0.99.
In the Boulware state, the sign of the RSET components does not change at both extremes of our radial grid or anywhere in between, except for a narrow region in the Q/M = 0.99 case.As the charge is increased, the posi-tive coefficient that controls the divergence of the RSET at r = r + decreases, but this does not change the fact that every component diverges ∝ f −2 , which only become finite at the horizon in the extremal case.
For the Unruh state, again the sign of the RSET components is independent of Q, but their overall magnitude diminishes with the charge, in consistency with this state describing a flux of black-body radiation for large r.At the event horizon, the T t t U and T r r U components decrease and increase with Q, respectively, whereas the T φ φ U component increases for Q/M = {0.4,0.8, 0.9}, decreases with Q/M = 0.95, and becomes negative for Q/M = 0.99, as in the Hartle-Hawking state.
In the extremal case Q = M , it was shown in [31] that the RSET of massless fields at the event horizon takes the same value as the RSET of a conformally invariant field in the Bertotti-Robinson spacetime [32].The magnitude of the RSET at r = r + varies more abruptly for Q/M > 0.9, hence we do not observe they approach the values obtained in the Bertotti-Robinson spacetime for 0.99.To observe this tendency, values even closer to M would need to be explored.

VI. COMPARISON WITH APPROXIMATIONS
Obtaining accurate results for the RSET, even in scenarios of high symmetry like the Reissner-Nordström black hole, proves to be a computationally expensive task: a large amount of highly precise Euclidean and Lorentzian field modes need to be calculated for each collection of parameters {M, Q, ξ, µ}.Once results are available for the RSET, we can find its backreaction on the background metric at O ( ) through Eqs.(1).If one insists in progressing along this line to find the backreacted metric at higher-orders in , it is necessary to compute, every time, new sets of modes propagating over the backreacted metric.Such a scheme proves to be unreasonably time consuming.
Analytic RSET approximations alleviate the difficulties behind computing the RSET and its backreaction since they bypass numerical mode calculations, expressing the RSET (almost) exclusively in terms of the metric components and their derivatives.This simplifies dramatically the complexity of semiclassical analyses at the cost of reducing the physical content encoded within the RSET itself.Furthermore, there is no unique or preferred approximate scheme available in the literature.Instead, we have multiple RSET approximations based on different physical principles.Despite this non-uniqueness, analytic RSETs should nonetheless reproduce the physics of their exact counterparts at least qualitatively, replicating the defining properties of the different vacuum states and yielding correct results at the asymptotic regions of the spacetime.In this section, we will review the various RSET approximations available in spherical symmetry and compare them with the exact results presented in Section V.
FIG. 3: RSET components in the Hartle-Hawking, Boulware, and Unruh states (first, second, and third rows, respectively) with ξ = 0.The blue, green, yellow, red and purple curves correspond to Q/M = {0.4,0.8, 0.9, 0.95, 0.99}, respectively.Left panel is T t t , middle panel is T r r and right panel is T φ φ .
A. Analytic RSETs in spherical symmetry

The Polyakov RSET
In spherical symmetry, we can obtain various analytic RSETs by fixing the field parameters of the theory, namely, the field mass and the coupling.Perhaps the most well-known example is that of conformally invariant fields {µ = 0, ξ = 1/6} in conformally flat backgrounds, where the RSET is fully determined by the local trace anomaly [33].In more generic spherically symmetric spacetimes, the essential features of the propagation of a massless minimally coupled scalar {m = 0, ξ = 0} in four spacetime dimensions can be captured by twodimensional models, described by the line element This connection between 4D and 2D physics manifests upon taking the near-horizon limit r → r + in Eq. ( 21) for the s-wave (l = 0) component of the field.The gravitational potential vanishes in this limit, and the (t, r) sector reduces to the two-dimensional free wave equation, which is conformally invariant (to see this explicitly, we avoid expanding in t in (19), see [34] for details).The emergence of this symmetry allows to express the 2D RSET in closed analytic form [35].These 2D expressions are then identified with a 4D RSET through where {a, b} run over the 4D spacetime indices while {A, B} run over the 2D spacetime indices, and P stands for the Polyakov RSET.The multiplicative factor 1/4πr 2 has been introduced to ensure covariant conservation in 4D and so that this approximation reproduces the adequate Unruh fluxes at infinity [see Eq. (66) below].In black hole spacetimes, the components are where f = (r − r + )(r − r − )/r 2 in the Reissner-Nordström spacetime and angular pressures vanish.The terms T a b P κ+ are temperature-dependent and relate RSETs in different vacuum states in a manner analogous to the integrals in Eqs.(37).Here, they amount to the Schwarzian derivative between null coordinates in two different conformal coordinate systems [34,36].The Boulware state is defined as the state for which such temperature-dependent terms vanish.For the Hartle-Hawking state they equal where For the Unruh state we find At spatial infinity, this component describes the usual Hawking flux emitted by an evaporating black hole.It can be easily checked that (40) reduces to the above expression in the massless case by ignoring backscattering (|B in ωl | = 1) and neglecting l > 0 multipoles.As it properly describes black hole evaporation, the Polyakov approximation has been extensively used in the literature [37][38][39][40].Note that the Polyakov RSET is illdefined at r = 0.This pathology motivates the search for regularized Polyakov RSETs that display a regular behaviour at r = 0 [40] while reproducing (66) at large distances.

The s-wave RSET
It is possible to incorporate into the Polyakov RSET the backscattering effects of the gravitational potential (which we had neglected via the near-horizon approximation) by considering a two-dimensional scalar coupled to a dilaton field (we refer the reader to [34] for details on this approach).This method is hybrid, in the sense that the resulting 2D RSET is non-conserved, such nonconservation being identified as the angular components of the following 4D RSET, where 4D time and radial components have been obtained through relation ( 62).The temperaturedependent terms now acquire the more involved form for the Hartle-Hawking state.For the Unruh state, these terms acquire a dependence on the time coordinate t and the resulting RSET is not covariantly conserved.

The Anderson-Hiscock-Samuel RSET
Dimensional reduction is not the only procedure to derive analytic RSET approximations.In 4D, there is the analytic RSET approximation derived by Anderson, Hiscock and Samuel [2] that incorporates the effects of field mass and curvature coupling.This approximation naturally arises from point-splitting regularization, which results in a separation of the exact RSET into two independently conserved analytic and numeric parts.The analytic portion -or AHS-RSET hereafter-gives the correct trace anomaly in the {m = 0, ξ = 1/6} case.This approximation yields a well-behaved RSET at r = 0, but in exchange it exhibits third and fourth order derivatives of the metric functions and depends on the arbitrary lengthscale .Expressions for this RSET, which we avoid showing here since they are lengthy and opaque, can be found in [2].The AHS-RSET can describe the Boulware and Hartle-Hawking states, but it cannot describe the Unruh state.
As was hinted in [2], the AHS-RSET is not an appropriate approximation for massive fields: it does not reproduce standard results in flat spacetime.Upon evaluating the AHS-RSET in Minkowski spacetime we obtain T r r AHS with T θ θ AHS , and λ is a positive parameter related to an infrared cutoff in some integrals from [2].In the Boulware state, λ can be absorbed in if the field is massless (being an arbitrary parameter otherwise), whereas in the Hartle-Hawking state, λ = κ + exp (−γ E ).In view of (69), we cannot make the AHS-RSET amount to a 4D thermal bath with temperature κ + /2π in flat spacetime by an appropriate choice of .Neither can we identify these components with a renormalization of the cosmological constant, as we did in section V for the exact RSET in the Boulware state.Therefore, we regard the AHS-RSET as an inadequate approximation for massive fields, and consider µ = 0 hereafter when referring to this approximation.
As the reader may have noticed from (69), the AHS-RSET contains an explicit dependence on the temperature.In the Reissner-Nordström spacetime, these terms take the form Notice how these terms differ from those given by the Polyakov (64) and s-wave approximations (68).These differences will have a major impact on the behaviour of approximate RSETs at the event horizon on the different vacuum states.

The DeWitt-Schwinger RSET
The work [2] also derived expressions for the RSET in the DeWitt-Schwinger approximation for fields of large mass.The corresponding expressions [obtained to O m −2 ] are local, have no information about the state of the vacuum, and can be used as an approximation to the complete RSET in situations where the field mass is comparable to the ADM mass M .
The components of the DS-RSET are where the coefficients in the sums are given in Eq. ( 83) in the Appendix.
In the next subsection we explore the accuracy of the aforementioned RSET approximations in the Reissner-Nordström black hole spacetime.

B. Comparing exact and approximate RSETs
Establishing comparisons between RSET approximations is complicated by the presence of renormalization ambiguities.For the sake of brevity, we will focus on the features of the RSET that are independent of these ambiguities.These are: the behaviour at the event horizon in the different states where, for massless fields, the RSET is independent of [see Eq. ( 54)]; the form of temperature-dependent terms, and the trace anomaly.
with vanishing angular pressures, whereas for the s-wave RSET The components of the Polyakov RSET are of the same order of magnitude as the exact RSET and have positive sign at r = r + until Q/M ≈ 0.94.Despite both approximations being regular in the Hartle-Hawking state, the temporal and radial components of the s-wave RSET have opposite sign.This discrepancy is alarming, especially considering that the s-wave RSET only incorporates backscattering effects to the Polyakov RSET, and that both approximations should agree at the event horizon, where the gravitational potential vanishes.It has been argued that such a discrepancy is due to the dimensional reduction anomaly [41][42][43] that accounts for the non-commutativity of quantization and dimensional reduction.
Among the various analytic approximations presented, only the Polyakov RSET can describe the Unruh state, for which diagonal RSET components are simply given by (65).Hence, its components are finite at the horizon, contrary to the divergences appearing in the exact RSET.
Turning now to the AHS-RSET, we have The magnitude of the AHS-RSET components shows excellent agreement with the exact results.Similarly, we observe a sign inversion for certain coupling values.In the extremal limit Q → M , the AHS-RSET becomes independent from ξ, in agreement with the exact results obtained in [31].Clearly, the AHS-RSET is perfectly regular at r = r + when the field is massless.For massive fields, it has an unphysical divergence at the horizon [2] caused by the failure of the WKB approximation there [44].
In fact, for a minimally coupled massless field, the expression above for T t t AHS HH = T r r AHS HH at r = r + is exact.To see this we first note that, using the results in [26] and Eq. ( 12), we have: Moreover only the n = 0 modes contribute to [g φφ W ,φφ ] r=r+ and for a massless field the radial mode functions p 0l , q 0l are given by Legendre polynomials.Hence it is straightforward to show using standard identities involving Legendre functions (see for example [2]) that the n = 0 mode sum contribution to [g φφ W ,φφ ] r=r+ vanishes, yielding: which for the Reissner-Nordström spacetime, evaluates to T t t AHS HH = T r r AHS HH with ξ = 0.The T φ φ HH component on the event horizon contains a contribution from the n = 1 radial modes and is therefore not amenable to a similar calculation, however we note in passing that a quasi-closed form event horizon expressions for all RSET components for a scalar field in the Hartle-Hawking state with arbitrary mass and coupling, are presented in [26].
In summary, the Polyakov and AHS approximations are in good qualitative and quantitative agreement with exact results for massless fields in the Hartle-Hawking state.The s-wave approximation, on the contrary, predicts wrong signs for the energy density and radial pressure.

The Boulware state
To obtain the Polyakov, s-wave and AHS approximations in the Boulware vacuum, we need to subtract the corresponding temperature-dependent terms in each case, resulting in an RSET that is singular at the event horizon.For example, the s-wave RSET returns the following divergent contributions at r = r + , For the analytic AHS-RSET we find In the ξ = 0 case, all RSET components in every analytic approximation diverge with the same signs at the event horizon, coinciding with the sign of the exact RSET for all Q values.For all the coupling values analyzed, the sign of the AHS-RSET also agrees with that of the exact RSET.
The leading-order divergence in the T  64) and ( 68)] scale differently from those obtained in 4D (70).In the extremal limit, all three approximations are singular at the extremal horizon [31,45].

Anomalous trace and massive fields
Another quantity that is worth comparing is the trace of the RSET.The Polyakov RSET predicts the correct trace anomaly in 2D.Due to the conformal symmetry of the dimensionally-reduced theory, the trace of the Polyakov RSET, is state-independent in 4D, thus finite and positive everywhere except when Q is in the range M ≥ |Q| > 8/9M .The s-wave RSET has a temperature-dependent trace instead obtained from the components (67).In the Boulware state, this trace is negatively divergent at r = r + for any Q, in stark contradiction with the trace of the exact RSET and the AHS approximation (see Eq. (82) below).This divergence in the trace is caused by the angular pressures (78).
For the AHS-RSET the trace in the Hartle-Hawking state is (82) In the conformally coupled case, Eq. (81) yields the correct trace anomaly.At the event horizon, (82) is positively divergent for ξ ≤ 1/6 and any charge, whereas it is negatively divergent for ξ > 1/6.Finally, Fig. 4 shows the trace of the exact RSET in the Hartle-Hawking state for fields of various masses, together with the trace of the DS-RSET.Clearly, the largemass approximation becomes better as the field mass increases.However, this approximation does not include temperature-dependent terms and, as such, will not approximate the exact RSET evaluated in the Boulware and Unruh states near the event horizon.

VII. DISCUSSION AND CONCLUSIONS
In spherically symmetric black hole spacetimes, we can employ Euclidean methods to compute the RSET in the Hartle-Hawking state.The extended coordinate method outlined in this article provides an accurate and efficient way of calculating the RSET in such situations.On the other hand, the RSET for the Boulware and Unruh states cannot be defined on the Euclideanized metric.While progress is underway in applying the extended coordinate method to metrics with Lorentzian signature, we adopt an alternative approach towards obtaining the RSET in the Boulware and Unruh states by taking advantage of the fact that the difference between RSETs in two different Hadamard states is finite.Hence we need only ever renormalize in one reference state.In this paper, we used the Hartle-Hawking state as the reference state where renormalization is performed using the extended coordinate method on the Euclidean slice; the RSET in other states involves the components in the reference state plus finite, rapidly converging integrals over products of Lorentzian modes.This way, we were able to generate results for the RSET in the three states, for various couplings and field masses in the Reissner-Nordström spacetime.To the best of our knowledge, our results for massive fields in the Unruh state are the first in the literature.
The scope of this method extends to spacetimes in which there is no preferred thermal state, i.e., situations where the Euclideanized line element (2) has no conical singularities.These encompass stellar spacetimes, where an auxiliary finite temperature state could be used as a conduit towards obtaining the RSET in the Boulware state.In stationary spacetimes that admit no Euclidean line element, as long as the RSET for a single state is known, it would still be possible to generate results for other states in those regions where both states are Hadamard.
Regarding analytic RSET approximations, we have seen that the Polyakov approximation [37] reproduces reasonably well the Hartle-Hawking and Unruh states, whereas it predicts a milder Boulware divegence at the horizon.The s-wave approximation [34] is clearly off in the Hartle-Hawking state and is unable to yield a covariantly conserved RSET in the Unruh state.Note that both approximations only describe massless, minimally coupled fields.The AHS-RSET [2] admits any coupling but does not give standard results in Minkowski spacetime when the field is massive.For massless fields, it is qualitatively correct in the Hartle-Hawking and Boulware vacuums.Finally, massive fields in a state regular at the event horizon are well approximated by the DS-RSET.
in the Hartle-Hawking, Boulware and Unruh vacuum states for different couplings and charge values.

4 FIG. 2 :
FIG.2: RSET in terms of r/M for various couplings and vacuum states for the Reissner-Nordström black hole with Q/M = 0.6.The blue, green, yellow and orange curves correspond to the couplings ξ = {0, 1/8, 1/6, 1/2}, respectively.Left, middle and right columns each represent T t t , T r r and T φ φ .Top, middle and bottom rows correspond to the Hartle-Hawking, Boulware and Unruh states.

1 .r 5 ++
The Hartle-Hawking state and Unruh states In the Hartle-Hawking state we find the following leading-order contributions for the Polyakov RSET at the event horizon, T t t P HH = T r r P HH = r + − 2r − 96π 2 O (r − r + ) ,