Exact equilibrium distributions in statistical quantum field theory with rotation and acceleration: scalar field

We derive a general exact form of the phase space distribution function and the thermal expectation values of local operators for the free quantum scalar field at equilibrium with rotation and acceleration in flat space-time without solving field equations in curvilinear coordinates. After factorizing the density operator with group theoretical methods, we obtain the exact form of the phase space distribution function as a formal series in thermal vorticity through an iterative method and calculate thermal expectation values by means of analytic continuation techniques. We separately discuss the cases of pure rotation and pure acceleration and derive analytic results for the stress-energy tensor of the massless field. The expressions found agree with the exact analytic solutions obtained by solving the field equation in suitable curvilinear coordinates for the two cases at stake and already - or implicitly - known in literature. In order to extract finite values for the pure acceleration case we introduce the concept of analytic distillation of a complex function. For the massless field, the obtained expressions of the currents are polynomials in the acceleration/temperature ratios which vanish at $2\pi$, in full accordance with the Unruh effect.


I. INTRODUCTION
The goal of thermal quantum field theory is to calculate mean values of physical quantities at thermodynamic equilibrium, henceforth denoted as thermal expectation values. This is a well known subject for the familiar thermodynamic equilibrium, described by the grand-canonical ensemble density operator: T 0 being the temperature and µ 0 the chemical potential coupled to a conserved charge Q, and Z the partition function. However, this is not the only possible form of global 1 thermodynamic equilibrium in special relativity. It has become a topic of interest the calculation of thermal expectation values for the general form of global thermodynamic equilibrium in flat space-time, which is described by the density operator: where the P µ are the total four-momentum operators and J µν the total angular-momentum boost operators; b is a timelike four-vector and ̟ an antisymmetric constant tensor which is known as thermal vorticity. Its form shows that the most general global equilibrium density operator in Minkowski spacetime involves all 10 generators of its maximal symmetry group, the proper ortochronous Poincaré group, with 10 constant coefficients. The operator (2) and its derivation from the generally covariant form as well as the relations between thermodynamic and hydrodynamic fields and the parameters b, ̟ have been discussed in detail elsewhere [1][2][3][4][5]. Here we just recollect that (2) implies that the four-temperature field β is given by: whose amplitude β 2 is the inverse of the local comoving temperature T and its direction defines a four velocity u. The thermal vorticity can be decomposed into two space-like fields as: where, in the case of global equilibrium only: A µ being the acceleration field and ω µ the kinematic vorticity field. The calculation of thermal expectation values of quantities such as the stress-energy tensor, conserved charge currents or spin density matrices with the operator (2) has a theoretical interest in the framework of statistical quantum field theory in general spacetimes [6][7][8] but it is also phenomenologically relevant for relativistic fluids undergoing very strong accelerations and rotations, such as the QCD plasma in high energy nuclear collisions [9]. Indeed, with b and ̟ depending on the space-time point x, the operator (2) is the lowest order approximation of the local thermodynamic equilibrium density operator [2,5,[10][11][12][13] which is needed to calculate those mean values in relativistic fluids.
For free boson fields, it can be expected that any thermal expectation values arises from the Bose-Einstein distribution function in its covariant form (Juttner distribution): f (x, p) = 1 e β·p+ζq − 1 (6) with β in (3) including all the effects of acceleration and rotation. In fact, it is known that a non-vanishing thermal vorticity implies the appearance of quantum corrections to the quantities obtained by integrating the (6). They have been calculated perturbatively [3,4,[14][15][16] for small values of ̟ with the operator expansion technique or the functional approach [17,18]. However, not much is known about the full global equilibrium expression except for the two special cases of pure acceleration [19,20] and rotation [21] where the exact solution is obtained by solving field equations in suitable curvilinear coordinates.
In this work, we present a method to calculate the exact expression of the distribution function in statistical quantum field theory in global thermodynamic equilibrium with rotation and acceleration without solving the field equations in curvilinear coordinates, that is the exact form of (6) with all quantum corrections. The method is based on the factorization of (2) and an iterative procedure to determine the thermal expectation values of creation and annihilation operators for imaginary thermal vorticity. The physical solutions are thereafter obtained extracting the analytic part of the formal series found, a mathematical procedure introduced in this work and named analytic distillation, followed by an analytic continuation to real thermal vorticity. We finally compare the obtained results with those calculated by solving field equations and find complete agreement.

Notations
In this paper we adopt the natural units, with = c = K = 1. The Minkowskian metric tensor g is diag(1, −1, −1, −1); for the Levi-Civita symbol we use the convention ǫ 0123 = 1. We will use the relativistic notation with repeated indices assumed to be saturated. However, contractions of a single index, e.g. β µ p µ will be sometimes denoted with a dot, i.e. β · p; similarly, the contraction of two indices, e.g. ̟ µν J µν will be sometimes denoted by a colon, i.e. ̟ : J. Operators in Hilbert space will be denoted by an upper hat, e.g. O.

II. THERMAL EXPECTATION VALUES OF ANNIHILATION AND CREATION OPERATORS
For free fields, the building block to calculate any statistical quantity (mean values, correlations) is the mean value of the combination of one creation and one annihilation operator of four-momentum eigenstates: Tr ρ a † (p) a(p ′ ) .
Its value depends, of course, on the density operator; if ρ = (1/Z) exp[− H/T ] the result is a well known one -the familiar Bose-Einstein and Fermi-Dirac distributions -and the derivation is based on the use of the known commutation relations between ρ and a, a † . However, if the density operator is (2), finding an exact expression is not trivial. Before setting out to do that, it should be first pointed out that any other combination of creation or annihilation operators, such as a † (p) a † (p ′ ) , a(p) a(p ′ ) or combination of creation/annihilation operators of particles and antiparticles in case of a charged scalar field, will vanish. This happens because the density operator (2), involving just the generators of Lorentz transformation, translations and charge, does not change the number of particles and antiparticles. Hence: a † (p) a † (p ′ ) = a(p) a(p ′ ) = b † (p) b † (p ′ ) = b(p) b(p ′ ) = a(p) b(p ′ ) = a(p) b † (p ′ ) = 0.
In this section we will present a general method to calculate (7) with the density operator (2). We will first set ̟ imaginary and factorize the density operator (2) by using a group theoretical method; then we will calculate (7) for imaginary ̟ by means of an iterative method and give the general solution in terms of a uniformly convergent series. Finally, we will discuss the analytic continuation to real ̟.
A. Factorization of the density operator As has been mentioned, we take advantage of Poincaré group relations which make it possible to factorize the density operator (2) into the product of two independent operators.
Let us start with the following very simple observation concerning the composition of translations and Lorentz transformations in Minkowski space-time. Let x be a four-vector and apply the combination T(a) Λ T(a) −1 , T(a) being a translation of some four-vector a and Λ a Lorentz transformation. The effect of the above combination on x reads: x → x − a → Λ(x − a) → Λ(x − a) + a = Λ(x) + (I − Λ)(a) = T((I − Λ)(a))(Λ(x)).
Since x was arbitrary, we have: where φ are the parameters of the Lorentz transformation. By taking φ infinitesimal, we can obtain a known relation about the effect of translations on angular momentum operators: The left hand side of (8) can now be worked out by using the above relation: Hence, combining (9) with (8), we have obtained the factorization: Now, since it holds by setting: and taking into account that: we have: Therefore, the right hand side of the eq. (11) becomes: Finally, the eq. (10) becomes: where we denotedṼ The eq. (12) can be read as the factorization of the exponential of a linear combinations of generators of the Poincaré group. For this reason, it must be derivable also by using the known formulae of the factorization of the exponential of the sum of matrices exp[A + B] in terms of exponentials of commutators of A and B. Indeed, it can be shown, by using the commutation relations of P and J, that one precisely gets the eq. (12) for any vector V and tensor φ, either real or complex. Hence, the formula (12) can be applied to factorize the density operator (2), neglecting Q for the moment, by setting φ = i̟: It is important to point out thatb(̟) is in general a complex vector, however, according to eq. (14), it is real when ̟ is imaginary, that is when we deal with actual Lorentz transformations.

B. Iterative solution with imaginary ̟
We can now take advantage of this factorized form to calculate the mean value (7). The idea is to seek the solution for pure imaginary ̟ first and then to analytically continue it to real ̟, the physical case. The reason is that with imaginary ̟ = −iφ one deals with an actual Lorentz transformation, for which commutation relations with creation and annihilation operators are known. Specifically, a Lorentz transformation operator Λ = exp[−iφ : J/2] acts as 2 : For the translation-like operators, the commutation relations can be calculated regardless of whether the vectorb is real or complex: Let us then write the analytic continuation of the relation (13) to ̟ = −iφ, so: By using the relations (15) and (16), the analytic continuation of (7) can be worked out as follows: where, in the last equalities, we have used the ciclicity of trace. The above derivation can be written by using the shorthand of replacing the trace with the density operator: We can now use the known commutation relation between creation and annihilation operators where ε = p 2 + m 2 , to obtain: where Λp stand for the space part of the four-vector Λp. This is the basic equation we have to solve. First, we observe that if Λ = I, we get back to the familiar thermal field theory and the (17) can be solved algebraically leading to the well known Bose-Einstein distribution function, multiplied by 2ε in view of our adopted covariant commutation relations. In the more general case, we will obtain a consistent solution of (17) by an iterative method. The idea is to approximate the left hand side with the second term of the right hand side with the Dirac delta, and correct the expression obtained by inserting the first term calculated with the previous approximation. Let us see how this is accomplished in formulae. We start with: This approximated solution implies: Now we can plug the above expression again in the (17) and get, after multiplying both sides by eb ·p ′ : From eq. (18) we can now calculate an updated approximated expression of a † (Λ(p)) a(p ′ ) : and plugging it again in eq. (17); we obtain: Iterating, we eventually obtain the solution of (17) in the form of a series: The expression (19) can be further transformed by using the orthogonality of Λ: Note that, the orthogonality relation holds for complex ̟ and complex vectors x, y, so this equality would not be affected by choosing ̟ real or imaginary. We can then rewrite eq. (19) as: It is possible to find alternative forms of the (20) without the appearance of Lorentz transformations in the exponent. Indeed, it can be shown that (see Appendix A), for generally complex ̟ such that Λ = exp[̟ : J/2]: implying: By substituting the (22) into the (20), we obtain its alternative form: Some comments are now in order. First of all, it can be checked that indeed the found series in either form (20) or (23) are actual solution of the equation (17) by direct substitution. Besides, it can be readily realized that they are an extension of the Bose distribution function by taking the limit ̟ → 0, which implies Λ → I. Since, from (14) b(0) = b one has, from (23): An important point is that, because of the Dirac deltas with argument Λp, the forms (20,23) for finite ̟ make sense only if ̟ is imaginary, that is only if Λ is an actual Lorentz transformation. Therefore, it is possible to make analytic continuations only after the integration in the momentum p ′ is done. Once the integration is carried out, the convergence of the analytically continued series is not trivial (see later discussion). We just remark that ifb(̟) is time-like and future oriented, each exponent in the (17) is negative, what is more apparent in the form (19).
To conclude, we observe that, should the density operator include a chemical potential term coupled to a conserved charge like in eq. (2), the iterative derivation can be repeated and one has one more exponential factor: An important question is whether the (20), or one of its equivalent forms (23) (24), is the unique solution of equation (17) or if there are possibly more solutions. In general, it can be readily shown that any solution S(p, p ′ ) of the (20) can be written as: where S 0 (p, p ′ ) is a particular solution of the (17), such as (24), and H(p, p ′ ) is the general solution of the associated homogeneous equation: We will see that the existence of non-trivial solutions of the homogeneous equation (25) depends on the Lorentz transformation Λ; if Λ is a rotation, we will show in section VII that there are no non-trivial solutions; conversely, if Λ is a pure boost, there are non-trivial solutions, as shown in section VI. Nevertheless, it is a general feature that solutions of the homogeneous equation (25) are non-analytic around the identity, that is for φ = i̟ = 0, as shown in section V. Therefore, the particular solution found by iteration (17) may, in principle, contain non-analytic contributions which would obviously hinder a process of analytic continuation from imaginary to real ̟. If we could single out an analytic part around Λ = I of the particular solution (17), the function would be unique and its analytic continuation possible by construction. This method of analytic distillation is discussed in section V.

III. THE COVARIANT WIGNER FUNCTION AND THE PHASE SPACE DISTRIBUTION FUNCTION
The expectation value of the quadratic combination (7) is all we need to calculate every field-related quantity, like the stress-energy tensor or any other conserved current. However, it is convenient to make use of a very useful concept, the covariant Wigner function [22]: where ": :" stands for the normal ordering of creation and annihilation operators and a density operator is understood in the mean value. Note that the covariant Wigner function is real, but it is not positive definite. The covariant Wigner function (26) allows to express the mean particle current as a four-dimensional integral: Plugging the free scalar field expansion: e −ip·x a(p) + e ip·x b † (p) (28) in the eq. (26) we get: In the above equalities, we have taken advantage of the symmetric integration in the variables p, p ′ . The expression (29) makes it apparent that the variable k of the Wigner function W (x, k) is not on-shell, i.e. k 2 = m 2 even in the free case. This makes the definition of a particle distribution function f (x, p)à la Boltzmann not straighforward in a quantum relativistic framework. From equation (29), we can infer that W is made up of three terms which can be distinguished for the characteristic of k. For future time-like k = (p + p ′ )/2, only the first term involving particles is retained; for past time-like k = −(p + p ′ )/2, only the second term involving antiparticles; finally, for space-like k = (p − p ′ )/2 the last term with the mean values of two creation/annihilation operators is retained. In symbols: Multiplying the previous expression by θ(k 0 )θ(k 2 ) selects the particle contribution to the covariant Wigner function and by θ(−k 0 )θ(k 2 ) the anti-particle one. It is thus possible to define the particle covariant Wigner function W + (x, k) and an antiparticle counterpart W − (x, k). Without loss of generality, we can confine ourselves to particles and write their contribution to the current in eq. (27) by using the (29): To obtain the last expression, we have taken advantage of the hermiticity of the density operator, implying: which makes it possible to swap the integration variables p and p ′ . The formula (31) identifies the particle distribution function or phase space density f (x, p) as the real part of a complex distribution function f c (x, p): where We will later see that the complex distribution function is a very useful tool. For antiparticles the distribution function f (x, p) is obtained from (32) replacing a † (p) a(p ′ ) with b † (p) b(p ′ ) . Thus, by using (27) and (31): which shows that the current can be written as an integral of four-momenta, which are on-shell when using the phase space densities and not on shell when using the Wigner function.
We can now turn to the calculation of the mean value of the stress-energy tensor. From the Lagrangian: one can obtain the so-called canonical stress-energy tensor operator: and so determine its normal ordered mean value : T µν C : . For this purpose, we again make use of the covariant Wigner function and the known expression [22]: Now, it can be shown by using the equations of motion of the free scalar field, that: Since, from the covariant Wigner function definition (26): : we finally obtain, by using (35) and (36): : We note in passing that, since T µν is a quadratic operator of creation and annihilation operators of momentum eigenstates, just like the Wigner operator, its normally ordered mean value corresponds to the subtraction of the vacuum expectation value from the plain mean value, that is: : In order to derive the relation between the stress-energy tensor and the phase space distribution function, we first write down the derivative of f c in the eq. (33): Now we can plug the (29) in (37) and by the same method used for the current, we can express the stress-energy tensor as a function of the phase space distribution function: : The above expression makes it apparent that the gradient terms are quantum correction to the "classical" expression of the stress-energy tensor in kinetic theory in that they require -after restoring the natural constants -, 2 factors to have the same dimension. Also, the phase space distribution functions may themselves have quantum corrections depending on . It should also be emphasized that the relation between stress-energy tensor and phase space distribution, unlike sometimes believed, depends on the particular stress-energy tensor operator. This relation indeed is non-invariant under pseudo-gauge transformations except at homogeneous thermodynamic equilibrium [23]. The equation (20) or its equivalent forms (23), (24) can now be used to calculate the complex phase space distribution function (33) for the general global equilibrium. It is important to keep in mind that the expectation values were obtained for imaginary ̟ and that an analytic continuation to real ̟ must be eventually carried out. By using the eq. (23), it turns out to be: where we have used the orthogonality of Λ and its linearity. Since: and: we have: where we have defined c ≡ ̟ · x and used the eq. (14) as a definition of tilde-transformed vector. Therefore, by using (41), (43) and (22), which hold for any tilde-transformed vector, we obtain: We can substitute this result into the (40):

Now note that, by definition ofb andc
where β is the four-temperature vector of global thermodynamic equilibrium in eq. (3). Another very useful expression of this vector is:β which is a consequence of (42) and (44). We can finally write the complex phase space distribution function as: For antiparticles, thef c (x, p) is obtained by replacing ζ → −ζ. We are now going to study the above series and highlight some of its main features.

IV. PROPERTIES OF THE DISTRIBUTION FUNCTION
From the function (46) we can calculate most quantities of interest, once analytic continuation to real ̟ is made. We first observe that in the limit ̟ → 0 the series is real and it boils down to the Bose-Einstein distribution function: where the last equality applies if ζ is lower than the mass so that the series converges. This suggests that the series expresses the familiar quantum statistics expansion and its first term is the Boltzmann limit. Otherwise, if ̟ = 0, the (46) is a series of analytic functions of the various ̟ µν taken as complex variables and if the series uniformly converges in some region, it defines an analytic function therein. In fact, in the physical case with ̟ real, the series in eq. (46) is not always convergent, as we will see, whetherb(̟) is time-like or not. While this does not prevent an analytic continuation of the phase space distribution function, it requires a careful analysis to obtain it. Particularly, it turns out that the properties of the series are different according to which components appear in the thermal vorticity tensor ̟ and a separate treatment is needed.
In the physical case, it is worth pointing out that Reβ is itself a Killing vector, just like β. Looking at the last equation in (45) we obtain: Each term of the last series is symmetric in µ ↔ ν if k is even and antisymmetric if k is odd. Thus, because of i k , the real part of the right hand side series selects antisymmetric terms and ∂ ν Reβ µ turns out to be an antisymmetric tensor, which proves our statement. Conversely, the imaginary part of ∂ νβµ is a symmetric tensor. This observation leads to the conclusion that f (x, p) is not a solution of the Boltzmann equation. From the eq. (46), withβ n ≡β(−n̟): If we now multiply by p µ we get: The first term on the right hand side vanishes as Reβ n is a Killing vector, while the second term in general is non vanishing.
It is also important to point out that the form (46) is indeed a resummation of all possible quantum corrections at all orders in of the Bose-Einstein distribution function. Looking at (45) it can be realized that, restoring natural constants, β is a classical vector taking into account the eqs. (3) and (4), whereas the k-th term in the series is proportional to k to make ̟ : J adimensional. Hence, an expansion in ̟ of the (46) with fixed β, such as (setting ζ = 0 for simplicity): is in fact the full semi-classical expansion in , which has been purposely restored for the occasion; each power of ̟ involves a corresponding factor to the same power. The series above can be rearranged as an asymptotic power series in or ̟, which, as we will see in sections VI and VII, is mostly divergent. Notwithstanding, it could be used to work out the quantum corrections of, say, the stress-energy tensor, with the eq. (39), at some fixed order in or ̟. It is worth remarking the similarity of this situation with that of the general, non-equilibrium, gradient expansions in relativistic hydrodynamics [24], whose the above expression is, in a sense, a special case (with vanishing shear tensor) because:

V. MATHEMATICAL INTERLUDE: ANALYTIC DISTILLATION AND SERIES RESUMMATION
The solution of the equation (17) found by iteration, the eq. (24), may not be unique if the associated homogeneous equation (25) has non-trivial solutions, as we discussed in section II. However, it is possible to show that if (25) has non-trivial solutions they are non-analytic about the identity, that is for To prove it, suppose a non-trivial analytic solution around φ = 0 exists, F (p, p ′ , φ) (without loss of generality, we can assume that φ is a single-valued variable). We can then write: and equate both sides of (25): We can now expand the right hand side around φ = 0 and equate each power of φ separately. For the constant term at φ = 0 we get: whose solution is F 0 (p, p ′ ) = 0. At first order we have, taking into account that F 0 (p, p ′ ) = 0: whence F 1 (p, p ′ ) = 0. Iterating, by taking derivatives of higher order, it can be shown that all functions F n (p, p ′ ) vanish, hence the only analytic solution of the equation (25) around φ = 0 is 0. In Appendix B a particular non-analytic solution of the (25) is found. It should be now quite clear that non-analytic terms in ̟ = 0 could be hidden in the particular solution (24) and all quantities derived therefrom. This makes analytic continuations quite problematic unless an unambiguous procedure to single out an analytic part is pinpointed. We call this procedure analytic distillation and we define it as follows: Definition. Let f (z) be a function on a domain D of the complex plane and z 0 ∈D a point where the function may not be analytic. Suppose that asymptotic 3 where n can take integer negative values. If the series formed with the common coefficients in the various subsets restricted to n ≥ 0 has a positive radius of convergence, the analytic function defined by this power series is called analytic distillate of f (z) in z 0 and it is denoted by dist z0 f (z).
It is worth dwelling upon this definition. First, note that asymptotic power series in a point of a function can be constructed [25] with iterative limits: However, these limits may depend upon the argument of the complex number z (e.g. exp(−1/z) for z = 0). Indeed, it is a well known fact that asymptotic expansions are different in different angular sectors centered in z 0 -the well known Stokes phenomenon [25]. The above definition prescribes to retain the part of the power series asymptotic expansion which is common to the various sectors or subsets D i covering the domain where the function is originally defined:ā if there is no common part, the distillate is simply zero. The strong requirement in the definition is the convergence of the asymptotic series, which is not usually the case, but it will apply in our cases of interest, as we will see. Also, from the definition it turns out that, if the function is analytic in z 0 , then dist z0 f (z) = f (z). If, on the other hand, the function has an isolated pole in z 0 , the distillation extracts the positive powers of the Laurent series.
A simple example of analytic distillation is that of a function which is the sum of an analytic function and a non-analytic function in z = 0, for instance: The analytic distillation of f in zero is because it can be readily shown that dist 0 e −1/z = 0, as its power series about zero in the domain Re(z) > 0 is vanishing.
The calculation of an analytic distillate involves that of an asymptotic power series, which is not always a simple task. When the function presents itself as a series of simple functions, however, there is an important result by D. B. Zagier [26,27] which can be cast as follows: Theorem 1. Let f : (0, +∞) → C be a C ∞ complex valued function on the positive real axis with the asymptotic power series in Then, the function defined by the series: 3 We denote asymptotic equality with the symbol ∼.
has the asymptotic expansion for x → 0 + : and ζ is the Riemann Zeta function. This theorem can be extended to a complex function. Let F (z) be a function of the complex variable z with the asymptotic power series about z = 0 in a domain of the complex plane including the real positive axis: Defining: we can apply the previous theorem to the function f (x) defined as: We have, from the above definition, for the function f (x): where Γ ϕ is the argz = ϕ line on the complex plane. Therefore, for fixed ϕ we have, by using the theorem 1 and the above results: A n ζ(−n)z n which concludes the proof. Now, according to the definition of distillate, we can conclude that, if the asymptotic expansion of G(z) is convergent: The integral I Γϕ , in general, depends on the phase of the complex variable z by construction. According to the definition of distillate then, under the same hypothesis of convergence: However, if the function F (z) is analytic in the complex plane, including z = 0, then it does not depend on the integration path because the integral on the arc Γ R vanishes for R → ∞ owing to the hypothesis of fast decay of f in the theorem (see figure 1). In this case, the integral I Γ φ does not depend on the argument of z, hence the asymptotic expansions is independent of the domain or angular sector and we have: where The theorem 1 has an important extension to asymptotic expansions with negative powers [26]: Theorem 2. Let f : (0, +∞) → C be a C ∞ complex valued function on the positive real axis with the asymptotic power series in Then, the function defined by the series: has the asymptotic expansion for x → 0 + : It can be seen that, if a −1 = 0, the asymptotic expansion for x → 0 + has a logarithmic term and it is thus not an asymptotic power series. Nevertheless, if a −1 = 0, this theorem provides, again, an asymptotic power series which is fit for analytic distillation according to the definition above. Furthermore, if a −1 = 0, the theorem 2 can be extended to complex variables with the same argument following theorem 1 with: n=−M A n z n is analytic, the integral can be done on any path and no dependence on the argument of z arises in the asymptotic expansion in the complex plane.

VI. STUDY OF THE SERIES AND COMPARISON WITH KNOWN RESULTS: ACCELERATION
We are now going to study a special case of global equilibrium which is obtained by setting the constant vector b and the antisymmetric tensor ̟ in eq. (2) to with T 0 and a real positive constants. The resulting density operator can be written: with K z ≡ J 30 being the generator of a Lorentz boost along the z axis. The combination H − a K z can be seen as the generator of translations along its flow lines [28]. This form, and its relation with the Unruh effect, has been recently studied in some detail in refs. [14-16, 19, 20]; we only mention here that the constant a is the constant acceleration of the comoving observer with velocity u = β/ β 2 with a hyperbolic world-line through the origin and T 0 the temperature measured by a comoving thermometer with the same world-line. From the equation (48) we get: In this case, there is a point Hence, from eqs. (42) and (14): and, because of the (45):β Note that: so that if a/T 0 = 2π, at the Unruh temperature, the above matrix vanishes and so doesβ, making the distribution function f c (x, p) in (46) independent of x and p, though divergent.
An important result for global equilibrium (2) is that the mean values of local operators depend on x only through the four-temperature vector, that is [20]: where the mean value is calculated with the density operator (2) with the parameters in the subscript and β(x) given by the eq. (3). For instance, for a scalar operator, this implies that it may depend only on the scalars that can be formed with β and its only non-vanishing derivative, that is ̟. Looking at the equation (4),(5), they are: However, in the pure acceleration case, we have w = 0 and we are just left with two independent scalars, β 2 and α 2 ; moreover, α 2 = A 2 /T 2 turns out to be uniform throughout [19]. Thus, once a scalar function is calculated in x = 0, it can be inferred everywhere by replacing β 2 (0) with β 2 (x). Therefore, the calculation of O(0) β(0)=b,̟ is sufficient to determine the whole function. For vector or general tensor fields, this conclusion holds, because it is possible to reduce the calculation of a vector or a tensor field to a set of scalar fields once decomposed onto the vectors β, the tensor ̟ and the metric tensor g. Altogether, in the pure acceleration case, we can determine every local mean value by setting x = 0 in the distribution function and study the series: where, from the (51) The resulting series is: where ε is the energy and p z the contravariant component of the momentum along the z direction. The above series is not convergent for a = 0 because the limit of the sequence is undefined. Note that the above series, as it stands, is not an asymptotic series because the terms are not an asymptotic sequence for a → 0 according to the definitions [25].
Since the simple replacement of the real acceleration in the phase space distribution function does not give rise to a convergent series, one may wonder whether an analytic continuation can be obtained through different methods, like e.g. Borel resummation. However, this would not cure the problem of (56). Moreover, we should consider that the non-convergence of the series in the physical case could arise from non-analytic terms in ̟ = 0, specifically a/T 0 = 0, which are hidden in the form (56), as discussed in some detail in section V. Indeed, since in the acceleration case a non-vanishing vector v = −ix 0 = (0, 0, 0, i/a) exists such that: according to the formula (51), there are non-trivial solutions of the homogeneous equation (25) and an instance can be easily found which is actually non-analytic for a = 0 (see Appendix B). Therefore, the idea is to single out the physical solution through the method of analytic distillation in a = 0. We will not apply this procedure to the distribution function (56) itself, but to some integrals thereof, because a comparison can be done with calculations of thermal expectation values carried out with a completely independent method.

A. Distillation and comparison with known results
Quantum field theory at thermodynamic equilibrium with the density operator (49) can be approached by solving field equations in Rindler coordinates and calculating the expectation values of relevant creation and annihilation operators. Since Wigner function is defined with a normal ordering, all expressions obtained in section III and, consequently, the distribution function (46) involve a normal ordering in Cartesian coordinates. Therefore, we should obtain the same thermal expectation values in Rindler-based calculations with the subtraction of the Minkowski vacuum contribution; particularly, all obtained results should vanish when a/T 0 = 2π, that is at the Unruh temperature [19].
To start with, we will compare the mean value of ψ † ψ for a single real uncharged scalar field in the massless neutral case (hence ζ = 0) which has been calculated analytically in Rindler coordinates [20]. According to the (26), (29) and the definition of f c (x, p) in (33), we have: where we have again setβ n ≡β(−n̟). Note that the momentum integral of f c (x, p) is real, which can be proved from the definition (33) taking advantage of the hermiticity of the density operator. As has been mentioned, it suffices to calculate the above integral at x = 0, and yet the series (56) expressing f c is divergent in the physical case. However, we can write the integrand as a series for imaginary acceleration, that is setting a/T 0 = −iφ and obtain the analytically continued mean value of the field squared : ψ 2 (x) : I : : where, by using (55): where the scalar product is the usual Minkowskian between complex vectors. Note thatb n is indeed real and future time-like for any real φ. With such a vector, the series in (57) is uniformly convergent for real values of φ and can be integrated term by term. The result is, for a massless field: therefore: : which has the correct known limit of T 2 0 /12 for φ = 0. We now have to go back to the physical case, i.e. to continue the above expressions to imaginary φ or real accelerations of the function defined by the series: However, as we have pointed out right before this subsection, the above function may include non-analytic contributions which should be subtracted away through the procedure of analytic distillation in φ = 0 described in section V. Indeed, the function S 2 (φ) is peculiar on the complex plane; it is analytic if Re φ = 0 because it can be readily shown that the series is uniformly convergent, whereas for imaginary φ one has: which is certainly divergent as 1/ sin 2 (n Imφ/2) ≥ 1; thus, the function S 2 (φ) is non-analytic everywhere on the imaginary axis, although it is convergent for φ = 0. Moreover, the series is not an asymptotic series for φ → 0 [25] and cannot even be resummed with the Borel method because zeroes of the denominators are overall dense. Nevertheless, for real φ, the series G 2 (φ) fulfills the requirements of the theorem 2 of section V with a −1 = 0, hence an asymptotic power series around φ = 0 can be found which can be extended in the complex plane according to the discussion following theorem 2 and the equation (47). However, the function has poles on the imaginary axis, so the integral term in the asymptotic expansion I Γ /z must be dependent on the argument of z and does not contribute to the distillate. This is confirmed by a separate application of the theorem 2 to the positive real axis and the negative real axis by simply changing φ → −φ in the function S 2 (φ). One finds, being G 2 (φ) an even function: Note that the asymptotic power series of G 2 (φ) is finite, because of the vanishing of the ζ(−n) for even n in the (47). This makes it possible to make an analytic distillation in φ = 0 of the function S 2 : The physical solution is then obtained by setting φ = ia/T 0 : : and, taking into account that β 2 (0) = 1/T 2 0 and a 2 /T 2 0 = −α 2 , we have, in a generic point x: : This result coincides with the exact solution in Rindler coordinates found in refs. [19,20] with subtraction of Minkowski vacuum contribution. Indeed, as we have pointed out in section III, the normally ordered mean value of an operator which is a quadratic combination of creation and annihilation operators of momentum eigenstates coincides with the mean value with subtraction of the vacuum expectation value. Therefore we expect, according to the general features of the Unruh effect and the conclusions of ref. [19], the normally ordered thermal expectation value (61), as well as any other quadratic operator in the fields to vanish precisely at the Unruh temperature, that is a/T 0 = √ −α 2 = 2π and this actually occurs. Furthermore, since (57) is an integral of the complex distribution function f c (x, p), we argue that this characteristic must extend to any integral of the distribution function after analytic distillation and continuation; this will be demonstrated for a large class of integrals in the Appendix C. As pointed out below eq. (52), for a/T 0 = 2π the equation (52) yieldsβ = 0, hence, for ζ = 0 a complex distribution function (46) which is a divergent constant independent of x and p.
In order to confirm the agreement found, we have compared the expectation values of the canonical stress-energy tensor with an exact calculation in Rindler coordinates, reported in Appendix D. If the distribution function were known, the stress-energy tensor could be calculated using the eq. (39), which select only the real part of a momentum integral. But we only know the distribution function through eq. (46) that is only valid for imaginary value of thermal vorticity. Therefore to obtain the physical expectation values we have to compute the momentum integrals by plugging the distribution function representation (46) into (39) for imaginary thermal vorticity, then continue it to physical values and at last take the real part: : where the imaginary expectation value is the series: : The momenta in the integrand of eq. (63) can be expressed as derivatives with respect toβ n : Now the momentum integral is the same as eq. (59) and we obtain: The derivatives can be worked out by using: and one eventually obtains: : where: As has been pointed out earlier in this section (see eq. (53)), it is sufficient to calculate the stress-energy tensor in x = 0. From the (52) one can write the vector fieldβ n (x) as a function of the Minkowskian coordinates t and z: whence:β and: The derivatives of theβ n needed to work out the (66) in x = 0 are: .
From the equations (65),(66) and using the above derivatives we get: Because of the symmetries of the density operator (49) and the fact that the stress-energy tensor can only depend on β and ̟ and not explicitely on x, according to the discussion about equation (53), it may only have this form [20]: where u µ = β µ / β 2 and the thermal functions ρ, p and A can only depend on the Lorentz scalars β 2 and α 2 . From the (68) and the (62) we then get the physical thermal functions: . The series are similar to the one in eq. (60) and they can be tackled in the same manner. The analytic distillation and the continuation to physical acceleration with a = −iT 0 φ of these series are discussed in Appendix C. The result is: The form of the thermodynamic coefficient in any point x can be inferred from the previous result simply by taking into account that β 2 (0) = 1/T 2 0 and a 2 /T 2 0 = −α 2 , so we have: These coefficients nicely coincide with those calculated by solving field equations in Rindler coordinates in the right Rindler wedge (reported in Appendix D), where the four-temperature is a time-like vector. Moreover, they all vanish at the Unruh temperature T 0 = a/2π that is for |α| = 2π, as expected for a quadratic operator according to the discussion following equation (61). These results precisely coincide with a perturbative expansion in acceleration to quadratic order in refs. [29] and to fourth order presented in refs. [14,15] for the same statistical operator in (49). This result bears out the observation made in refs. [14] that the vanishing of the stress-energy tensor at the Unruh temperature is achieved, in the perturbative expansion in a/T at the fourth-order which led the same authors to conclude that the exact expression must be polynomial [15], namely that the (70) were indeed the complete solutions. These authors also pointed out that these same results can be derived within a non-perturbative approach in a quantum field theory over a space with a conical singularity [16], see for instance [30] and other references in [16].

VII. STUDY OF THE SERIES AND COMPARISON WITH KNOWN RESULTS: ROTATION
We are now going to study a special case of global equilibrium which corresponds to system rigidly rotating along the z axis and which is obtained from eq. (2) by setting: where ω and T 0 are positive constants. More in general, without specifing the axis of rotation, we can choose with ω constant vector. It can be shown that ω has the physical meaning of a constant angular velocity [31] and T 0 is the temperature measured on the axis of rotation. This state is also known as rotational thermodynamic equilibrium. The density operator becomes: Because of equation (14), it turns out thatb(̟) = b. In the physical case with real ̟, that is real ω, the eq. (45) can be rewritten as:β where R(iω/T 0 ) = exp[ω · J/T 0 ] is the rotation around the axis defined by ω by an imaginary angle iω/T 0 . Thus: and, consequently:β Therefore, the analytically continued distribution function is, from eq. (46): which converges only if (x × p) ·ω < 0, that is only if the orbital angular momentum along the ω direction is negative. Otherwise, the series is divergent, which seems to imply that there is an infinite number of particles with (x × p) ·ω > 0.
To understand the source of this problem, and to prove that the distribution function found is the actual exact solution, it is very convenient to diagonalize the density operator (72), following the method of ref. [32]. Settingω =k as the direction of the z axis, this can be done by using single particle states with eigenvalues p T , p z (transverse and longitudinal momentum) and the eigenvalue of the angular momentum M which takes on integer values. The relation between momentum eigenstates and eigenstates |p T , p z , M is well known and leads to the relations between creation operators: along with its inverse: where ϕ is the azimuthal coordinate of p. One can then plug the (75) in the field expansion in cartesian coordinates (28) and obtain the corresponding solution in cylindrical coordinates. The operator on the exponent of (72) can then be written as: where H 0 is a divergent constant. Unlike the Hamiltonian, the spectrum of the operator (77) is not bounded from below because the single particle energy can be as low as the mass m, whereas the angular momentum eigenvalue M can take on any value, so ε − M ω has no lower bound. This causes most statistical mean to be hopeless divergent because the probability of large occupation numbers is exponentially large when ε − M ω < 0 which is clearly the case when M is positive and the energy is low enough; this is the origin of the divergence of the distribution function f (x, p). This problem is due to the absence of boundary conditions for the field; if the scalar field has Dirichlet boundary condition ψ = 0 at a radius R < 1/ω the spectrum is bounded from below because the transverse momentum eigenvalues are discrete and depending on M , and the following inequality holds [32]: Curing the divergence inherent in the rotational equilibrium without field boundary conditions goes beyond the scope of this work. However, we can check that the solution found (23) is the same as the known one for the basis (p T , p z , M ). By using the same method as for a † (p) a(p ′ ) one can calculate a † (p T , p z , M ) a(p ′ T , p ′ z , M ′ ) with the density operator (72); it is a much easier task for J z is diagonal on the eigenstates created by a † (p T , p z , M ) [32]. The result is diagonal in the mean values: Note that the series is convergent only if ε > M ω, see above discussion. For an imaginary ω/T 0 = −iφ, that is for an actual Lorentz transformation, the sum (78) is always convergent and the result is: The relation (78) can be used to calculate a † (p) a(p ′ ) by means of the (75), for ω/T 0 = −iφ: Let us now use take (24) with Λ = exp[−iφJ z ]. Since the rotation around the z axis leaves p T and p z unchanged, the vector Λ n p has cylindrical coordinates p T , ϕ + nφ and p z . Hence: Therefore, beingb = b = (1/T 0 , 0) for rotations, as we have seen, the eq. (23) becomes: which is the same as (80). The proof is concluded. A remarkable feature of the calculation with imaginary angular velocity is that the series involved in (79) is always convergent. Thereby, the problem of the divergence due to the lack of lower bound of the operator (77) is evaded by going to imaginary ω and it is possible that a suitable analytic continuation to real ω yields finite results. It should be stressed, however, that the physical meaning of the expressions found is limited because the divergence should be cured in the physical case with different methods.
This proof has another noteworthy consequence, that it is that there is no non-trivial solution of the homogeneous equation (25) for the density operator (72); in other words the (24) is the only solution of the equation (17). The reason is simply that a † (p) a(p ′ ) is the only possible expression corresponding to the the solution (79), which is in turn the only possible expression of a † (p T , p z , M ) a(p ′ T , p ′ z , M ′ ) which can be obtained algebraically [32] in the basis diagonalizing the density operator. This is consistent with what we described as a sufficient condition to obtain a non-trivial solution of the homogeneous equation (25) as, in the case of rotation, there is no vector v such that b = b = (I − Λ)(v) with b as in eq. (71) (see Appendix B).

A. Comparison with known results
We now move to the calculation of thermal expectation values similar to those of the previous section. First we define a suitable tetrad for the rotational equilibrium withω =k. Keeping in mind the decomposition (4), we define the following four-vectors (see also [3,4]): The tetrad {u, α, w, l} is an orthogonal, non normalized basis of Minkowski spacetime built upon the four-temperature β and thermal vorticity ̟ and can be used to decompose all vectors and tensor fields. Moreover, denoting with r 2 = x 2 + y 2 the distance from the rotation axis we have that the time component of the four-velocity u reads: and: Formulae (81) allow to transform the widely used r 2 , ω, T 0 variables used in the rotational equilibrium problems into the Lorentz invariants discussed in the previous section. Note that, unlike in the acceleration case, in the rotational case the calculation of the thermal expectation value of a scalar operator at x = 0, that is r = 0 is not sufficient to know its value everywhere because at r = 0 we have α 2 = 0 and the dependence on one of the three independent Lorentz scalars is lost. We now turn to the mean value of the field squared in the massless neutral case, see eq. (57). With φ = iω/T 0 , we haveβ and the series (57) becomes: : For real φ the series is absolutely and uniformly convergent, because it is bounded by the series of 1/n 2 ; this is an expected result according to the discussion on the convergence of the series (79). However, for complex φ the series is densely divergent. To show this, let Re φ be a rational multiple of 2π, so that Re φ = (N/L)2π with N and L irreducible integers. The n-th denominator of the series in (83) vanishes when n = KL with K integer and if: making the series in (83) a divergent one. Since K and N are arbitrary integer numbers, the curves (84) are dense in the complex plane. Altogether, the function (83) cannot then be defined except on the real axis. The divergence of the series (83) for non-real φ is evidently a consequence of the divergence of the distribution function (74) for (x × p) · ω > 0 and its physical reason has already been discussed. Even though there is no proper domain of existence of the complex function (83), thus no asymptotic expansion for complex φ, we can, notwithstanding, extract a finite analytic distillate of the function defined by the series (83) for real φ by using the theorems of the section V, confining ourselves to the positive and negative real axis and keeping in mind that the physical meaning of the solution will be limited (this point will be resumed later on).
We start rewriting the right hand side of (83) in a form which is suitable for the application of the theorem 2. By using the auxiliary real parameter t we have: : Exchanging the limit and the series, for real values of φ and t, is possible because as a function of t the (85) is a uniformly convergent series of continuous functions. We can now apply the theorem 2 to the function g(φ) = R(t, φ)/φ 2 which leads to the following result: where the upper sign applies to φ > 0 and the lower to φ < 0 and: Note that, again, the term a −1 = 0, which excludes logarithmic terms in the asymptotic series. The asymptotic power series of R(t, φ) near φ = 0 (86) has, again, a finite number of terms due to the fact that there are no odd powers in the expansion and ζ(−n) = 0 for even n ≥ 2. Finally, by taking the limit t → φ we obtain: : ψ 2 : I = lim t→φ R(t, φ) ∼ 1 12 It can be seen that the right hand side of (87) can be expanded as a power series in φ around φ = 0 for |φ| < 1/rT 0 by using the geometric series. Hence, altogether, the (87) is the resummed asymptotic power series of the function : ψ 2 : , about φ = 0, which is suitable for analytic distillation. The distillation just removes the last term in (87) and, going to physical angular velocity with φ = iω/T 0 , we have: : ψ 2 : = 1 12 Finally, taking advantage of relations (81), the (85) and (87) yield: : which coincides with the expression found in the pure acceleration case (61). This is indeed a remarkable and unexpected result. Starting from two considerably different statistical operators (49) and (72) we have obtained, after a long mathematical derivation, the very same expression of the mean value of the field squared in terms of the Lorentz invariant thermodynamic variables. On one hand, this confirms the covariant structure of the theory and how deep is the connection between angular velocity and acceleration in relativity. On the other hand, it is quite surprising that in the rotating case the (89) does not feature any dependence on w 2 .
Like in the pure acceleration case, we are going to calculate the thermal expectation values of the stress-energy tensor.
To compare with expressions found in literature [21,33] we considered the general form of the stress-energy tensor operator of the scalar field, which is trace-less for m = 0 for ξ = 1/6: where the canonical stress-energy tensor is given by the (34). The mean values can be calculated with the same method of the previous section for pure acceleration: we study the series associated to a complex mean value : T µν ξ : I for ω = −iφ/T 0 , we take the analytic distillate and continue the result to physical angular velocity and lastly we take the real part. We just need to modify the relation (65) to account for the term proportional to ξ in the (90); this is simply done by adding the following terms to the quantities in (66): , ∆t ξ = −ξ 1 n 2 (β n ·β n ) .
From the (73), with φ = iω/T 0 , we obtain: Plugging these expressions into the (65),(66) and (91) we can work out every component of : T µν ξ : I and express it as a series; the explicit forms are given in Appendix E. Each series has the same mathematical features as the one in (83) and can be resummed likewise, that is by using theorem 2 followed by the analytic distillation and continuation to the physical values. The details of the procedure and the results of each components in terms of the physical rotation are given in Appendix E. Here we report the results in terms of the thermal coefficients associated to the stress-energy tensor. Indeed, the stress-energy tensor can be decomposed on the tetrad {u, α, w, l}, which can only depend on the Lorentz scalars (54): with ∆ µν = g µν − u µ u ν . While a symmetric stress-energy tensor has ten degrees of freedom, the thus defined scalar coefficients are eleven. The redundancy of (92) is owing to the fact that the tetrad {u, α, w, l} is no longer a basis for r = 0 because there α = l = 0. For this reason we introduced a term p∆ µν which would be redundant at finite r. By projecting the stress-energy tensor onto the suitable pair of vectors (see Appendix E), we obtain the expressions of each scalar function. The isotropic pressure is identified in the projection of T µν along the direction of l, which has the following continuation to physical values of ω: Since l 2 = −α 2 w 2 , the identification of p and G l from the above equation is unambiguous. Once p has been determined, all the other coefficients are obtained by subtraction and one finally finds: Taking into account the relations (81), these results can be compared with those calculated in refs. [21,33] by solving the Klein-Gordon field equation in rotating coordinates without enforcing boundary conditions at finite radius r and, in spite of the different methods used, precise agreement is found.
Remarkably, by setting w = 0 and l = 0, as well as ξ = 0, the resulting canonical stress-energy tensor is the same function of β, α found in the pure acceleration case and particularly the functions in (70) are the same functions of β 2 and α 2 . This confirms the previous finding for : ψ 2 : and the deep relation between the two examined cases. It appears, because of the covariant structure of (2), that from the study of the pure rotation, one can also deduce the exact solution of thermal expectation values in the pure acceleration case, which is certainly a remarkable and unexpected fact. It should be pointed out, however, that the physical meaning of the (93) is limited in the case of pure rotation; these results have been obtained by removing the inherent divergences pertaining to the density operator (72) with the method of analytic distillation, and not by setting appropriate boundary conditions or other physical meaningful constraints which would make the operator H − ω J z bounded from below. Indeed, the terms in α 2 /β 2 , α 4 /β 4 , w 4 /β 4 and α 2 w 2 /β 4 are actually independent of T 0 and do not vanish in the limit T 0 → 0 what would be expected in the pure rotation case for a well-behaved, bounded from below operator H − ω J z in the formula (72).

VIII. SUMMARY
In summary, we have derived a general exact form of the phase space distribution function and the thermal expectation values of local operators of the free quantum scalar field in the most general case of thermodynamic equilibrium in Minkowski space-time, with the four-temperature being a Killing field, that is including rotation and acceleration. The presented derivation does not make use of the solutions of the Klein-Gordon equations in curvilinear coordinates but it is just based on the plane wave expansion of the field in Minkowski space-time and it is therefore suitable for the general case including both rotation and linear acceleration. The crucial steps of the derivation are a factorization of the general density operator (2) using Poincaré group algebra and an iterative method to obtain the thermal expectation values of quadratic combinations of creation and annihilation operators. The general form of the phase space distribution function has been written as a formal series including all quantum corrections to the classical term. We have studied the series in two major cases of non-trivial equilibrium, the pure acceleration and the pure rotation and compared with exact known results obtained solving Klein-Gordon equation in Rindler and rotating coordinates respectively. In the case of pure acceleration, the iterative method introduces undesired non-analytic terms to be subtracted while in the rotation case the absence of boundary conditions at finite radius < 1/ω implies unavoidable physical divergences. The resummation of the series in these two cases required the introduction of a new operation on complex functions, defined as analytic distillation, in order to extract the analytic part in the pure acceleration case and analytically continue the finite solution for imaginary angular velocity in the pure rotation case. In the former case, the method of analytic distillation leads to expressions in the massless case which are the same as obtained in Rindler coordinates and, remarkably, automatically vanish at the Unruh temperature; mathematically, a new class of complex polynomials is generated by this physics problem, which all vanish for z = 2πi. In the case of pure rotation, on the other hand, it should be pointed out that the analytic distillation and continuation provides a finite result which is in agreement with analytic calculations in literature, but with limited physical meaning.

ACKNOWLEDGMENTS
We are greatly indebted to D. Dorigoni for letting us know about Zagier's theorem and his work on Lambert series. We warmly thank D. Basile and D. Dorigoni for illuminating discussions about asymptotic series and the extraction of an analytic part. We thank F. Colomo and C. Dappiaggi for equally useful discussions.

Appendix A: Recurrence formula
We want to work out Λ −sb (̟) for a non-negative integer s and Λ = exp[̟ : J/2] being ̟ complex in general. By using the eqs. (14) and the expansion of Λ, we have: This expression can be worked out as follows: where we set n = h + k. We can write the last obtained expression with binomial coefficients: whereb is given by the (14) with φ = i̟, then: solves the equation (25) for any function G(p ′ , φ). Note that, sinceb is real for real Λ, v(φ) must be real as well. Since: taking into account that v is real, we have: This functional equation, applying to any p, p ′ is solved by the combination: If Λ is a rotation, thenb = b (see Section VII) and since b is proportional to the time unit vectort, there is no real vector v(φ) solving the (B1). If, on the other hand, Λ is a pure boost, say along the z axis, then we have that (B1) has non-trivial solutions. By using the equation (51) and the definition of x 0 = (0, 0, 0, −1/a) (see section VI) it can be checked that: is a solution of the (B1) with φ = ia/T 0 . In this case, it turns out that, according to the above general proof, H(p, p ′ ) is non-analytic in φ = 0.
We want to study such integrals for the pure acceleration case, with imaginary acceleration a/T 0 = −iφ and demonstrate explicitly that -in the massless neutral case -after analytic distillation and continuation they give rise to expressions vanishing at the Unruh temperature T 0 = a/2π. This might be a consequence of the distribution function (46) being independent of p and x when evaluated at the Unruh temperature, as shown in in section VI. Since that is a general feature of the distribution function, we can expect that the expressions for the massive and charged case will also be vanishing, as dictated by the Unruh effect. Moreover, it is sufficient to calculate the integrals (C1) in x = 0 because the whole dependence on x of any tensor field is fully constrained (see discussion after eq. (53)).
By using eq. (59) and eq. (67), it can be realized that the equation (C1) in x = 0 gives rise to linear combinations of the following series: As shown in the text specifically for S 2 (φ), these are series of analytic functions if Re φ = 0 and they are uniformly convergent in the same domain, so they define analytic functions when Re φ = 0. Conversely, these series are divergent whenever φ is purely imaginary. Nevertheless, we can obtain asymptotic power series about φ = 0. As far as the series S 2m+2 (φ) in (C2) is concerned, it can be shown that its asymptotic expansion about φ = 0 is given by a series of odd powers of φ, hence its contribution to the physical expectation values, after distillation, continuation, and extraction of the real part, vanishes.
We are then left with S 2m+2 (φ) only. With φ > 0 real and positive, we can use the theorem 2 to obtain an asymptotic power series of the function G 2m+2 (φ): By using the generalized version of the Bernoulli polynomials B n (t) defined by [34] x the function f can be written as a power series about φ = 0: which converges for any complex φ with |φ| < 2π [34]. Owing to the parity of the function f , the coefficient B (2m+2) n (m + 1) appearing in (C3) must be vanishing for odd n. The conditions of the theorem 2 for its application to the function G 2m+2 (φ) are fulfilled and we have, in its notation: as well as where in the last step we iteratively used [35, Formula 1 §1.4.5, p. 146]. Therefore, the asymptotic expansion of G 2m+2 (φ) reads: 2m+2+n (m + 1) (2m + 2 + n)! ζ(−n)φ n .
Since a −1 = 0 and the function f in eq. (C3) is analytic in the region Re φ > 0, we can extend the above expansion to complex φ in the same region according to the discussion following theorem 2 (see also figure 1 in section V). Notice that the above sum effectively stops at n = 0 because, for positive n, the Zeta function is non-vanishing only when n is odd, but in this case the generalized Bernoulli polynomial coefficients vanish. Then, by taking advantage of the relation between the Riemann Zeta function and the Bernoulli numbers the asymptotic expansion of the hyperbolic series can be written as: For instance, in the cases m = 0, 1 the asymptotic expansions are: Repeating this procedure in the region where Reφ < 0, we find that the asymptotic expansion of S 2m+2 (φ) is the same as in eq. (C4) with the opposite sign of the term proportional to φ 2m+1 , which is hence removed by the distillate. Now we can use the analytic distillation method to the function S 2m+2 (φ) and continue it to imaginary values of φ, that is to a real acceleration. Like in section VI the result of distillation are the polynomials: which are to be continued to imaginary φ. We now prove that the polynomials (C5) have two zeroes in φ = ±2πi for any m ≥ 0, hence their analytic continuation for real accelerations a/T 0 = −iφ vanish precisely at the Unruh temperature. If φ = 2πi is a zero of B 2m+2 then the following identity must be true: To prove the identity (C6), we start showing that for every integer m > 0 the generalized Bernoulli polynomials have a zero in m: With z ∈ C and α, β ≥ 0 integers, the following ratio of Euler Gamma functions has an exact power series representation [36]: Then, for m > 0 and setting α = β = m, the above equation becomes: If we evaluate the previous expression in z = 0 and we use the reflection formula for the Gamma function, we have: For z = 0 the sum (C8) contains just one term, hence we have: (m) = 0, which proves the equation (C7).
To proceed towards the proof of (C6), let us now consider the function which involves the generators of the Bernoulli polynomials and has a power series expansion: which is convergent for |x| ≤ 2π. We are going to show that, for any integer m ≥ 0: which proves the identity (C6). Consider the even terms of (C10), i.e. d 2k . Splitting the product in (C9) as x 2 x e x/2 e x − 1 2m+2 , it can be realized that the second term is odd and does not contribute to the d 2k while the first term gives rise to the generalized Bernoulli polynomials: From this expansion and eq. (C7) it follows that for every integer m ≥ 0 the coefficient d 2m+2 vanishes: 2m+2 (m + 1) = 0.
To prove that d 2m+2 equates the right hand side of (C11), we need to write the full power series expansion of f in (C10), which can be obtained with the product of the series of the two functions enclosed in the brackets of (C9). The first factor in the right hand side of (C9) can be expanded as: while the second factor is the generating function of the generalized Bernoulli polynomials of order 2m + 2: x e x/2 Hence, by making the Cauchy product of (C12) and (C13) in the (C9) we obtain: Equating this expression with the (C10), we obtain the expression (C11) for the vanishing coefficient d 2m+2 , which finally proves the identity (C6). Taking advantage of the (C6), the polynomial (C5) can also be written as: which manifestly shows that B 2m+2 vanishes at φ = ±2πi.
Appendix D: Stress-energy tensor for a massless scalar field in Rindler coordinates We calculate the thermal expectation value of the stress-energy tensor of the massless scalar field with the density operator (49) by using the solutions of the Klein-Gordon equation in Rindler coordinates, in order to compare with the results (70). The task will be accomplished by taking advantage of calculations presented in refs. [19,20].
The symmetries of (49) constrain the mean value to be of the form (69). On the other hand, from quantum field theory and using the equation of motion for a neutral massless field, the canonical stress-energy tensor may be written as: We kept the covariant derivatives for future convenience. Combining this expression with (69) the following expression are obtained: where u µ is the four-velocity β µ / β 2 and α µ = β 2 A µ where A µ is the acceleration field, see section I. Note that in this section the subscript R specifies that normal ordering applies to Rindler's creation and annihilation operators and not to Minkowski's, like it is understood in the rest of the paper. Since they have non-trivial Bogoliubov transformations, the two normal-ordering are not equivalent. From the four-temperature in (50) it ensues (Cartesian components): where z ′ = 1 + az and k = (z ′2 − t 2 ). The thermal expectation values of the fields can be calculated from the solution of the Klein-Gordon equation in Rindler coordinates. In fact, this calculation has already been partially carried out in [20], where it was found that: with α 2 = α µ α µ . With these results, the (D1) become: and we just need to calculate one more thermal expectation value. Before doing that, it is useful to remind some useful relations in Rindler coordinates. The relation between Cartesian and Rindler coordinates (called here τ and ξ) is given by: t = e aξ a sinh(aτ ), z ′ = e aξ a cosh(aτ ).
Using these expressions, we can show that: Also, in Rindler coordinates, k 2 = a −2 e 2aξ and β 2 = e −2aξ T −2 0 . Now we can work out the thermal expectation value in eq. (D2c): : α µ ∇ µ ψ 2 : R = β 2 : a µ ∇ µ ψ 2 : R = β 2 a 2 e −4ξa : d ψ dξ In the right Rindler wedge, the scalar massless field is given by [37]: K n (z) being the modified Bessel function of the second kind and k T and x T are the transverse momenta and coordinates. Using the mean value of the number operator in the right Rindler wedge [19]: T 0 − 1 and the formula: we obtain the integrals: where the sum runs over all possible combinations of signs ± in the order of the Bessel functions. The integration over k T can be done by means of the known integral [38]: which applies to Re a > 0 and Re λ < 1 − |Re µ| − |Re ν|. After integrating over ω we get: : and, putting it in the (D2), we have: so that, finally: The last step is to calculate the stress-energy tensor with subtraction of the Minkowski vacuum, which was implied in the normal-ordering used throughout the paper. As has been discussed in detail in ref. [19], the subtraction of Minkowski vacuum corresponds to the subtraction of the normally-ordered Rindler expression at the Unruh temperature which is such that T 2 U = −A µ A µ /4π 2 . So, for the Minkowski normally-ordered quantities (the subscript is dropped) we have: which agree with the equations (70).
As in the case of the series (83), the poles of these series cover the whole complex plane. However, we can still obtain a proper asymptotic power series if we restrict φ to take on only real values. Notice that the variable φ never appears alone but always multiplied by either n or T 0 . So, we can replace φ with an auxiliary variable t every time it appears with a coefficient T 0 . By means of this replacement, new functions R µν (nφ, t, x) are defined and the thermal expectation values of the stress-energy tensor can be obtained by taking the limit: as it was done in the main text for : ψ 2 : I . Indeed, these series are of the same kind of that in eq. (85). Also, it can be realized that for real values of φ, t and x they are all uniformly convergent series of continuous functions and, consequently, we can exchange the limit and the sum: : T µν ξ (x) : I = lim t→φ ∞ n=1 R µν (nφ, t, x).