Hawking Radiation from Universal Horizons

The persistence of a suitable notion of black hole thermodynamics in Lorentz breaking theories of gravity is not only a non-trivial consistency test for such theories, it is also an interesting investigation {\em per se}, as it might help us identifying the crucial features at the root of these surprising laws governing such purely gravitational objects. In past investigations, controversial findings were presented in this sense. With the aim of settling this issue, we present here two complementary derivations of Hawking radiation in geometries endowed with universal horizons: a novel feature of back holes in Lorentz breaking theories of gravity which reproduces several properties normally characterizing Killing horizons. We find that both the derivations agree on the fact that the Hawking temperature associated to these geometries is set by the generalized universal horizon peeling surface gravity, as required for consistency with extant derivations of the first law of thermodynamics for these black holes. We shall also comment on the compatibility of our results with previous alternative derivations and on their significance for the survival of the generalized second law of black hole thermodynamics in Lorentz breaking theories of gravity.


I. INTRODUCTION
Local Lorentz invariance of space-time (LLI) is at the basis of our present understanding of Nature. It is a fundamental symmetry of the standard model of particle physics as well as a funding pillar of General Relativity. Nonetheless, the search for a definitive theory of quantum gravity has led in the last decades to a renewed questioning about this symmetry and its role in Nature. First of all, being the Lorentz group non-compact we cannot say we have tested it completely, and we also know that LLI is strictly related to ultraviolet divergences in quantum field theory, being the latter a direct consequence of the assumption that the spectrum of field degrees of freedom is boost invariant. Moreover, hints of possible departures from Lorentz invariance came by several quantum gravity scenarios and, albeit non-conclusive, further stimulated the study of the phenomenological consequences of these departures from standard physics (see e.g. [1,2] for extensive reviews).
Missing a definitive theory of quantum gravity, it has been a common approach to take a bottom up attitude and adopt an effective field theory approach to Lorentz violations by adding systematically Lorentz breaking operators to different sectors of the Standard Model and then use current observations to constrain them [1,2]. Lorentz breaking can be limited to the boost sector, in a Riemannian geometry, at the cost of introducing a preferred vector field (normally a unit, timelike one) that selects a preferred frame 1 . This in turn, together with the desire to preserve background independence, led to the development of extensions of General Relativity which include, also in vacuum, the dynamics of an "aether" field and not just of the metric. The most general quadratic action for a metric and the aether is the one of the so called Einstein-Aether theory of gravity [5]. It is possible however to further require that such a vector selects as well a preferred foliation: this can be easily enforced by requiring it to be the gradient of a scalar and hence "hyper-surface orthogonal" and vorticity free. The action, in this case, takes the form of the low energy action of Hořava Gravity (which in principle entails not only mass dimension four operators but also several others, up to mass dimension eight, to guarantee power counting renormalizability).
It is worth stressing that while we do have fairly stringent constraints on deviations from Lorentz invariance in particle physics (see e.g. [2]), we do not have comparable constraints on the gravitational sector. Indeed, currently detected gravitational waves are too low energy to constrain in any way terms of mass dimensions larger than four, and even for the latter, they severely constrain only two of the three basic parameters of the theory [6]. Fortunately, when dealing with a physical theory it is not only possible to test it by using observations and experiments but it is also possible to probe its self-consistency.
Black holes represent the purest gravitational objects in nature and their thermodynamic aspects are still nowadays considered one of the few phenomena allowing us to have a glimpse into the possible merging of the gravitational and quantum realms. It is then not surprising that they can also play a crucial role in testing the self-consistency of Lorentz breaking theories of gravity. Furthermore, proving or disproving the persistence of black holes' thermodynamics in theories with Lorentz invariance would provide invaluable insight to what lies at the root of this tantalizing feature of gravity, and hence it would be relevant beyond the actual interest in such alternative theories with a preferred frame. Eddington 2 famously said that an unavoidable consistency test for a theory is to check whether the second law of thermodynamics holds true. This test seemed however miserably failed by the aforementioned Lorentz breaking theories of gravity when it was realized that violations of LLI seem to directly lead to violation of the Generalized Second Law (GSL) of black hole thermodynamics. The basic mechanism, in this case, is based on the fact that if different fields on a black hole space-time are endowed with different limit speeds, then in general they will experience different horizons and ergoregions. This mismatch can then be used to generate a violation of the second law, for example, by constructing a classical and quantum perpetuum mobile, or by devising a setup whose only outcome is to lower the total entropy of the universe [7][8][9][10].
While this could have spelled doom for this whole class of theories, in a surprising twist it was soon recognized that black hole solutions in these frameworks have an internal structure that is quite different from the standard black holes of General Relativity. First of all, it was shown that classical arguments for GSL violation can be circumvented once the dynamics of the non-relativistic theory (for gravity and matter) is properly taken into account [11]. Also, for what concerns Hawking radiation based arguments, it appeared soon clear that due to the possibility of non-relativistic dispersion relations for matter fields and gravitons, in such theories the Killing horizon does not capture the notion of the causal boundary as in standard black-hole space-times, as it can be realized by including the aforementioned higher order terms in the matter action (so to have not only different limit speeds but also superluminal signalling). Moreover, it was recognized that, static, spherically-symmetric black-hole solutions, in both Einstein-Aether and Hořava theories, contain a special hypersurface that acts as a genuine causal boundary because it traps all possible signals, even those traveling at arbitrarily high speeds 3 . This special hypersurface has been called the "universal horizon". Indeed it was also proven that the typical peeling structure shown by null rays at the Killing horizon in General Relativity is reproduced now at the universal horizon by the generalized causal rays moving under the influence of the metric and the khronon field with the peeling being characterized (albeit not completely fixed) by the surface gravity of the universal horizon, κ uh [14].
It is evident that such finding might play a crucial role in the restoration of the GSL in non-relativistic frameworks as the presence of a universal temperature associated to the UV completion of the framework considered in [7][8][9][10] would solve at the root any possibility of construction of a perpetuum mobile, at least based on Hawking radiation. Indeed, further support for a restoration of black hole thermodynamics via universal horizons was lent by the finding that a generalized first law of black-hole mechanics can be associated to universal horizons [15][16][17], instead of the usual Killing horizon, and by requiring that the temperature associated to the universal horizon has to be fixed by their surface gravity κ uh . Noticeably, such association was derived independently in [15] using the tunneling formalism of Ref. [18], and later on generalized to arbitrary truncated power law dispersion relations, ω ∼ p N for large p, in [17,[19][20][21]. 4 Despite this growing evidence, still many points have remained unclear. In particular, if the conjectured temperature associated with the universal horizon would be relevant for observers outside the black hole, or if instead the emitted radiation would be somehow reprocessed at the Killing horizon in an energy and species dependent way. In this context, the results reported in [24] seemed to surprisingly contradict the accumulated evidence for a universal horizon thermodynamics, as they implied that the temperature associated to Hořava black holes was, as usual, the one set by the peeling surface gravity of the Killing horizon (see e.g. [25] for a detailed discussion of different notions of surface gravity).
An extra motivation to understand the thermal properties of black holes in Einstein-Aether and Hořava Gravity comes from holography in Lifshitz systems through the AdS/CFT duality. These are systems similar to CFTs but which exhibit an anisotropic scaling between time and space. Initially, they were proposed to be described by bulk theories combining a Proca field and a dilaton field. The former condenses along a direction -in a similar fashion to the Aether here -thus breaking boost invariance and leaving a non-relativistic theory. However, it was later shown by [26] that some properties related with conformal anomalies that were not easily captured by these duals were emerging naturally in Hořava Gravity, opening thus the window to a Lifshitz/Lifshitz holography [27]. In this context, a complete understanding of thermality within the gravitational theory in the bulk could lead to important developments and applications for Lifshitz theories at finite temperature on the boundary.
In order to settle the issues described above, in this paper we shall rederive the Hawking temperature associated to spherically symmetric black holes in Hořava gravity using two different methods. We shall show that this is indeed associated to the generalized κ peeling as found in [19][20][21]. Also, we shall show that the different outcome of the investigation reported in [24] is not due to any technical issue but rather to a different choice of the vacuum state.
The plan of the paper is the following. In chapter II we shall revise black holes in Lorentz violating theories. In chapter III we shall review the usual derivation in a relativistic setting of Hawking radiation in a collapsing shell geometry. In section IV we shall apply the same method, by a suitable generalization of physical rays, to the case of a geometry endowed with a universal horizon. We shall show that in this case, the relevant temperature for Hawking radiation is the one set by the universal horizon. To investigate further the compatibility of this result with the investigation reported in [24] we revisit their calculation in chapter V, showing that their and our results can be reproduced with a suitable different choice of the vacuum state. Finally, we provide in chapter VI an extension of our results to more general dispersion relations, and we discuss their implications in the conclusions.

II. LORENTZ VIOLATING BLACK HOLES
We will consider the following action describing Einstein-Aether gravity where Here G is the Newton's constant and R the scalar curvature of the manifold. The couplings c i are dimensionless and parametrize the departures from General Relativity in terms of the dynamics of the aether vector U µ , which is forced to have unit norm U µ U µ = 1 by the Lagrange multiplier λ, thus breaking boost invariance. The coupling c 2 is constrained to be positive to ensure the absence of ghosts [28]. Here we are interested in static black hole solutions, which are only known analytically for aether fields that are hypersurface orthogonal. In that case, c 1 can be always absorbed into redefinitions of the rest of the couplings, so we will fix c 1 = 0 from now on. Incidentally, in that case, this action corresponds to the low energy limit of Hořava Gravity, known as khronometric theory 5 [29], in which the aether is described in terms of a single scalar degree of freedom, the khronon T Note that U µ U µ = 1 is automatic from this and therefore solutions enjoy a preferred time direction given by the integral lines of U µ . As a consequence, the theory enjoys a reduced symmetry group with respect to General Relativity. Instead of full Diffeomorphism invariance, the dynamic is invariant only under those diffeomorphisms which preserve the global time direction, dubbed FDiff [29].
In (1) we have set a possible cosmological constant to vanish. We do this for simplicity since we are interested in understanding Hawking radiation. Although formally this forbids the application of our results to contexts like AdS/CFT, we expect them to be easily translated to those cases. The role of a cosmological constant for the black hole solutions that we are going to study in the following would be only to regulate the large distance behavior of the metric and therefore, the dynamics of internal structures, such as horizons, whose characteristic length is typically much smaller than the curvature radius of a possible AdS embedding, should not be affected.
Regarding the viability of (1) as a realistic theory of Gravity, it must be noted that although Lorentz violations are highly constrained in the matter sector of the Standard Model [2], current bounds on gravitational dynamics are still not as competitive. The allowed parameter space has been quickly reduced in recent years, but there is still room for departures to be accommodated. In particular, the couplings in (1) are constrained to satisfy |c 3 | 10 −15 from the observation of multi-messenger signals from GW170817 [6,30,31], while Solar System tests [32][33][34][35] demand |c 4 | 10 −7 . An additional theoretical constrain, demanding that black holes moving slowly with respect to the aether remain regular except for their central singularity, requires c 3 = c 4 = 0 exactly [36]. On the other hand, measurements of the abundance of elements during Big Bang Nucleosynthesis [37][38][39] bound |c 2 | ≤ 0.1.
Here, however, we will be looking to a different region of the parameter space, not in the physical region but advantageous for our purposes here, given by c 4 = 2c 3 = −2c 2 . For this particular value of the coupling constants, the action (1) admits an analytic spherical static black hole solution given by the Schwarzschild metric with dΩ 2 = dθ 2 + sin 2 θdφ 2 and which will thus allow us to perform computations in an analytical manner, without having to resort to numerical methods, necessary for other values of the couplings. However, the main feature of the solution that leads to the phenomenon of interest for our work here is the presence of the universal horizon, which is present for generic values of all c i . Thus, we expect our analytical findings to also hold in the whole parameter space, at least schematically. Finally, this is the same solution considered in [24], and thus it allows us to lay a direct comparison of our results and methods against those of them.
The metric (3) is complemented by an aether configuration given by and its orthogonal vector As previously advertised, this aether configuration is hyper-surface orthogonal and it is normed to unity, while its orthogonal vector satisfies S µ S µ = −1. Note that there is a sign arbitrariness in choosing S µ which we fix by demanding it to approach S µ = δ µ r for large radii. The metric has a Killing horizon at r kh = 2µ for particles moving at light speed (c 2 = 1) and a time-like Killing vector given by as well as the usual curvature singularity when r = 0.
However, this is the end of the similarities with the relativistic case. Due to the presence of the aether, the solution also contains a universal horizon, sitting at At this radius, the integral lines of U µ wrap over themselves, becoming a compact surface of simultaneity. Moving over the r = µ ≡ r uh surface requires then to move at infinite speed and therefore, no signal can escape from it, thus becoming an event horizon for any motion. Signals moving at speeds c > 1 can penetrate the Killing horizon and eventually escape, but once they cross the universal horizon, they are trapped behind it. The schematic Penrose diagram of the geometry can be seen in figure 1. In the following, we will work with a two-dimensional toy model by throwing away the angular coordinates and focusing on the sub-manifold defined by t and r.
Motion within this space-time is a subtle issue. For free-falling observers moving at the speed of light, things are identical to the usual relativistic setting. They follow null trajectories of the metric (3). However, due to its non-relativistic character, the theory also admits observers which move at different speed c and therefore, in the general case, position dependent speeds c(r) larger than the speed of light.
As previously discussed, the black-hole solutions that we consider are hyper-surface orthogonal and therefore the integral curves of U µ define a global time coordinate which is shared by all observers. Thus, when writing dynamical equations, time derivatives must be identified with derivatives along U µ while space derivatives are projected onto the orthogonal hyper-surfaces where γ µν = g µν − U µ U ν is the projector onto the hyper-surfaces orthogonal to U µ . FDiff invariance then allows for terms which contain only time derivatives or only space derivatives, thus admitting fields with nonrelativistic dispersion relations from explicit violations of LLI. One must be careful, though. Including more than two time derivatives increases the number of initial conditions required to solve the general Cauchy problem for the system and will inevitably lead to instabilities in the form of ghost degrees of freedom. Therefore, we choose to keep our action second order in time derivatives and include only higher orderirrelevant in the IR -operators constructed with spatial derivatives.
In the present work we will thus consider a model of a scalar field coupled to the background geometry and the aether in the following way where we have included the first possible higher derivative operator which breaks Lorentz invariance, weighted by a scale Λ. The corresponding equations of motion are where the factor of 2 is there for future convenience.
In coordinates adapted to the global time defined by the aether, where 1 − 2µ r 1 2 dt AE = U µ dx µ , the equations are explicitly second order in time derivatives and thus one can quantize the scalar field in the standard way, introducing modes with momentum k µ which can be decomposed as and satisfy the dispersion relation, obtained in the eikonal approximation The scale Λ thus signals the magnitude of the momentum at which Lorentz violations become apparent in the motion of φ. Note that the existence of the aether, allowing for the terms with only spatial derivatives, is crucial in order to have a modified dispersion relation.
Before going further, we should make an important remark about the quantization of these modes. In a relativistic theory, it is well known that one can construct different equivalent quantum field theories for the same action by considering different coordinate charts and thus defining the Hamiltonian by using different time coordinates. However, the theories considered here enjoy a preferred time direction, given by the integral lines along the aether field U µ as previously commented. This is important because only when the time direction is chosen to coincide with the preferred direction, the action remains second order in derivatives [40]. This means that only when this time coordinate is used to define the Hamiltonian, the theory remains unitary in flat space-time. Otherwise, the action will contain up to four derivatives along any other time direction, which will lead to the presence of ghosts. Thus hereinafter we will always define the generalized momentum of the modes by using the preferred time coordinate as which thus always ensures the local absence of ghosts in the spectrum.
Through this work, and in order to be able to carry out an analytical analysis, we will rely on perturbation theory on Λ −1 , retaining the next-to-leading order contributions and in some situations, when required, also the next-to-next-to-leading order terms. Under this approximation, the dispersion relation thus becomes Being polynomial, using this dispersion relation will simplify our computation hugely. Later we will discuss the implications of retaining the full dispersion relation (13).
Another important point on understanding the motion of observers in this space-time is the construction of physical free-falling trajectories, i.e the rays of massless fields coupled not only to the metric, but also to the aether and therefore endowed with a modified dispersion relation. These serve not only as the natural clocks on the geometry but also as probe observers moving under the sole influence of g µν and U µ . In order to understand their motion, we consider a ray of the scalar field φ and note that its velocity vector can always be decomposed as [14] V where the sign indicates whether the ray propagates outwards or inwards in the geometry and c(r) = dω dq is the group velocity of the ray. Since we are dealing with a non-relativistic theory, the precise form of c(r) is frame-dependent and must be computed case by case. Note that although we have not displayed it explicitly, c(r) is not only space dependent but it will also change with the energy of the ray as we will see in a moment. Specifying to t and r, we can combine the components of V µ to get which allows us to define the point-wise effective causal-cone coordinates of the ray as The reader can check that these reduce to the usual light-cone coordinates u, v for the case of a relativistic dispersion relation. Following the integral lines ofū andv is equivalent to following the trajectory of outgoing and incoming rays in this geometry. At every point, the pair (ū,v) defines the causal cone of an observer communicating with signals that travel at speed c(r).
At this point, it is worth noting that although the rays satisfy the modified dispersion relation (13) locally at every space-time point, none of the components of k µ are constant along the trajectories. Both ω and q are functions of the radial coordinate r. Nevertheless, they can be related to conserved quantities easily.
Since neither the metric nor the aether depend explicitly on the time coordinate t, it is easy to check that the Killing frequency Ω = χ · k is conserved. Contracting (12) with the Killing vector we thus find This allows us to study the possible solutions to the two branches of this equation -the two branches in ωin terms of r. Their qualitative behavior is shown in figure 2. At distances larger than the Killing horizon, there exist only two solutions for positive Ω, which correspond to the expected positive energy modes as given by a unitary quantum field theory. However, for radius in the range µ < r < 2µ the parabolas are deformed and cross the Ω = 0 axis. We thus find four positive energy solutions. In particular, for r close to the universal horizon at r uh = µ, these solutions clearly separate in two sets, whose behavior is shown in figure 3. Two of the solutions remain soft and cross the universal horizon with finite momentum, while the other two become of order Λ. The extra two modes that now become physical correspond to the orange parabola in figure 2 and therefore they enjoy negative aether energy and, as a consequence, negative norm in the quantum theory. However, they should not be interpreted as ghosts in the usual manner, but instead as the negative energy partners to positive modes that end up trapped behind the universal horizon.
The behavior of the momenta can be computed perturbatively both in Λ and (r − µ) by choosing an appropriate ansatz and using (19). For the soft modes, we take where the subscript distinguishes the branch of the dispersion relation.  Retaining only the first sub-leading correction both for large Λ and r − µ we find for the momenta and group velocities c 0 Note that both modes are regular on the universal horizon, collapsing onto the line Ω = q, which allows them to cross the horizon. We thus identify these modes as incoming modes. Following the discussion in [24] it can be seen that the mode in the would-be positive branch -the branch corresponding to positive energy outside the Killing horizon -corresponds to a physical incoming mode while the would-be negative one represents a mode which tries to escape the gravitational well but eventually falls back.
The other modes have momentum of order Λ. Therefore we change our ansatz to  [24].
which gives, when solved perturbatively where c Λ ± is the group velocity of the modes.
These modes actually diverge on the universal horizon, which prevents them from crossing it. They move infinitely close but still at a finite distance from the universal horizon. Following [24] and [14] we identify them with modes that linger for a long time around the Universal Horizon, thus highly blue-shifted. The would-be positive mode escapes the geometry while the would-be negative mode can be understood as incoming from the inner region and attaching to the horizon. Note that one has to be careful when analytic continuing them through r = µ = r uh , owing to the presence of the absolute value. Finally, it is worth to point that although the mode itself is of order Λ, the solution is still under perturbative control, since it is given as an analytic expansion on Λ −1 with finite and non-vanishing radius of convergence.

III. HAWKING'S RADIATION FROM A COLLAPSING SHELL: THE RELATIVISTIC CASE
Let us for the moment forget about the existence of the aether and the modified dispersion relation and review in this section a standard approach to derive Hawking's radiation from a Schwarzschild black hole. Thus, for the moment we keep the Schwarzschild metric (3) but consider free-falling observers to be light-rays moving along null trajectories, with dispersion relation obtained in the eikonal limit In order to derive the thermal properties of the horizon, we will follow the classic treatment by [41,42]. We will consider a matter configuration in form of an infinitely thin shell moving towards the center of the geometry following a time-like trajectory r = R(τ ), starting from r = R 0 , and where τ is the proper time of the shell. At early times and far enough from the shell, the space-time can be taken to be flat everywhere, and therefore one can construct the usual Minkowsky vacuum |0 M there. For later times, and at any radius r = R(τ ), space-time is divided in two regions. The outer region of the shell will be described by the metric (3), while we take the interior, without any loss of generality, to be flat Eventually, the shell will cross the Schwarzschild radius 2µ and will form a black hole.
We will perform an experiment of sending rays along null trajectories towards the interior of the geometry from I − . These rays will cross the shell, hitting its center and being re-emitted towards the exterior geometry, eventually arriving to I + . Incoming rays will experience a blue-shift, due to falling into the gravitational well. If the shell was static, then this blue-shift would be exactly compensated by the corresponding red-shift of the reflected ray when exiting the geometry. However, in a collapsing situation, the difference between the radius of the shell when the ray enters and exits will generate a non-negligible contribution. Our approach to compute the form of the rays (which is tantamount to the state of the modes of a field φ) after the black hole has formed will be then to evaluate the shape of our probe rays after they hit I + . For sufficient late times, there will be a last ray which is able to enter the shell and exit it just before the horizon forms. The shape of this ray will be asymptotically close to that of the fundamental modes in the newly formed black hole geometry, from which we can obtain the flux of Hawking's radiation.
We will thus from now on follow the trajectory of a free-falling ray from I − to I + , described in the exterior of the shell at early times by the simple function where v = t + r * is the null coordinate of the metric (3), with r * = r + 2µ log r 2µ − 1 being the tortoise coordinate, and Ω is the Killing frequency, defined as the energy of the wave at r → ∞.
This ray will move towards the center of the geometry and will eventually cross the shell at r = R(τ ). At that point, it will start propagating in the flat interior geometry, which will be described by a different pair of null coordinates U, V . In order to be able to further follow the trajectory we thus need the relation between inner and outer coordinates, β(V ) = v(V ), so that the ray becomes Let us postpone for now the discussion of the explicit form of β(V ) and assume that we know it for the moment.
The inner null coordinates are defined by dU = dt − dr, dV = dt + dr, (32) and therefore, by integrating and imposing boundary conditions on the shell we have This allows us to place the center of the geometry at V = U − 2R 0 .
The ray will continue moving through the flat region until it arrives to the center of the shell at r = 0, where it will get re-emitted. By imposing reflecting boundary conditions, the out-going wave reads The ray will now travel outwards and cross the shell to get out to the Schwarzschild geometry again. Defining α(u) = U (u), we finally have which can be propagated to I + without any additional modification. Note that u corresponds to the second null coordinate in the exterior region, defined by u = t − r * .
As we have seen, in order to obtain the form of the wave at I + we need to compute the matching functions α(u) and β(V ). By differentiating with respect to the proper time of the shell, we have Luckily, we are only interested in the form of the rays after the horizon forms. Thus, we can just focus on their late time behavior, when the shell gets close to the Schwarzschild radius. We thus expand its trajectory around this point, where τ • and ξ are the time and (absolute value of the) speed of the shell at horizon crossing. Note that the actual speed ∂R ∂τ = −ξ must be negative, indicating that the shell is collapsing. Thus ξ > 0.
Using the metrics to obtain the differentials in (36), the two equations can be solved in the limit r → 2µ by keeping only the leading order contributions. We get where the constant in front of V is determined in terms of ξ, µ, R 0 and τ • but its specific form is irrelevant. C 1 and C 2 are integration constants. When both solutions are combined we find that the reflected wave takes the form where we have redefined the integration constants to absorb the constant in front of V . As it can be seen, the wave gets an exponential red-shift, controlled by the size of the Schwarzschild radius 2µ. We thus observe that while an observer in I − would see vacuum -because the mode that we have sourced in the geometry is a plane wave of positive energy -, an equivalent observer sitting at I + will instead observe a flux of particles coming out from the black hole. Hence, we have identified Hawking's radiation.
What is left is to determine the shape of the spectrum of the flux of particles in the radiation observed at I + . This can be done by decomposing the reflected wave onto the plane modes which serve as a basis in this region of space-time where α ΩW and β ΩW are the corresponding Bogolubov coefficients.
In order to compute them, however, it is more convenient to invert the relationship and rewrite the modes defined in I + in terms of the coordinates in I − by using the inverse of the functions α(u) and β(V ) -so that the logarithm in α −1 can be cancelled against an exponential-and using the inner product in I − Things are straightforward from here and more details can be found in [43]. It might happen, however, in more complicated settings, that the coefficients α ΩW and β ΩW contain a divergence related to the number of degrees of freedom in the radiation quanta 6 . This can be eased by noting that the divergence is shared between both coefficients and thus its ratio must be finite. In particular, we have This, together with the normalization condition for the modes dΩ|α ΩΩ | 2 − dΩ|β ΩΩ | 2 = 1, implies that the number of particles per mode in the radiation which corresponds to a Bose-Einstein distribution with temperature

IV. HAWKING'S RADIATION FROM A UNIVERSAL HORIZON
We return now to the situation of a black hole in Einstein-Aether Gravity, where we follow the motion of a field with an anisotropic dispersion relation (13) in the aether frame. Our aim in the following will be to reproduce the previous computation of the horizon temperature in this case.
We will consider, as in the previous section, an infinitely thin shell following a time-like trajectory r = R(τ ) towards the center of the geometry 7 . However, in this case, and unlike in the relativistic case, we will not stop tracking the shell after it reaches the Killing horizon. The reason is that now observers and rays travelling at speeds larger than the speed of light are allowed. While in a relativistic setting any ray that crosses the Killing horizon is forever trapped within its interior, that is not the case anymore. Here we have rays that move at arbitrary speed, so the process of crossing the Killing horizon and escaping from its interior is possible. Instead, it is the universal horizon which serves as an event horizon trapping any ray in its interior. Therefore this will be the ultimate region down to which we can follow the evolution of the shell. After the shell crosses the universal horizon at τ = τ • , we will consider that the black hole is formed. Note also that the aether field is completely regular at the universal horizon, when r = µ, laying along the timelike direction of space-time. Thus, we can continuously glue a flat interior to the shell without any discontinuity. A deeper analysis of this junction can be found in [24] 8 .
Another important difference to discuss here are the features of the asymptotic null infinity regions for incoming and out-coming rays, before and after crossing the shell. Indeed, high energy rays can have superluminal group velocities and in principle start in the past, and reach in the future, not only the low energy null asymptotics I ∓ respectively, but they can also start and reach spatial infinity i 0 . We shall hence construct our modes using the previously introduced effective causal coordinates from (18) which include the possibility of superluminal motion and hence of incoming rays from spatial infinity. For what regards the outgoing ones we shall see a posteriori that the spectrum of rays that arrive at large radii after crossing the shell is exponentially red-shifted at late times, thus peaked at low energies as usual. The propagation of these rays will then be close to the general relativistic one, and therefore I + will be approximately equal for all rays.
We start by describing the incoming ray at I − . In this region, the energy of the wave corresponds to the Killing frequency. Using the effective causal coordinates from (18), the incoming wave reads Again, in the asymptotic region r → ∞ space-time is flat and we can construct the usual Minkowsky vacuum there |0 M .
As in the relativistic case, this wave travels inwards towards the shell following the effective causal coordinate, i.e. a free-falling trajectory. Once it hits the surface of the shell, however, it will start propagating along the corresponding coordinate in the inner flat space. Afterwards, it will bounce back in the center of the geometry and, after crossing another portion of flat space, it will escape the shell and travel up to I + . In an identical way to the relativistic case, the reflected wave, once it arrives to I + , will read The coefficientc denotes the fact that the waves propagate at a speedc different from the speed of light within the inner flat geometry. The corresponding causal coordinates are thus now and the center of coordinates corresponds toV =Ū − 2cR 0 . The concrete value ofc could be determined by the speedc = c(q)| r=R(τ ) of the wave at shell-crossing but, being constant, its particular value is irrelevant for our computation here.
We are thus left with the task of computing the matching functions α(ū) and β(V ) between the inner and outer coordinates. As in the relativistic case, we can differentiate with respect to the proper time of the shell τ , getting Let us first focus on α(u). From the metrics and the explicit form of the inner null coordinates we get so that whereṘ = dR dτ .
Let us note a striking difference with respect to the relativistic case. Since the ray is accelerating and enters the shell with different speeds at different time snapshots, we find an explicit dependence in the group velocity, contained within ζ + , which changes with r. Since we are interested in the late time epoch, when the shell is asymptotically close to the universal horizon, we will only need to evaluate the behavior of ζ + in this region. In order to do that, we need the relation between the momentum q(r) and the Killing frequency Ω in (19) for the particular mode that we are tracking. Since α(ū) is the matching function for the coordinates following the trajectory of the out-going ray, after being reflected at the center of coordinates, we need to insert here the group velocity c(r) = dω dq for out-going rays, that we can evaluate in the overlap region when the reflected ray approaches the shell at r ∼ µ + 0 − . Thus, from (26) we get Going back to (52) and expanding the trajectory of the shell close to the universal horizon at leading order in (r − µ), which corresponds to leading order in (τ − τ • ), we get so that we finally have where we have evaluatedṘ = −ξ and we remind thatc is a constant.
The last step is to find a relation between R(τ ) andŪ in order to have a differential equation that we can solve. This can be done easily by using the shell trajectory and the metric with C 0 an integration constant. With this we finally get This can now be easily integrated to yield where C 1 is another integration constant and we have redefined C 0 to absorb constant terms into it.
We turn now to computing β(V ). Using again the metrics, we find that the differential equation governing the behavior of β(V ) can be written as In this case, we must take into account that β(V ) is the matching function for the coordinates of the ray incoming into the shell. Therefore, we must compute ζ − by evaluating the mode q 0 + in (21), which gives We plug this together with the expansion of the trajectory of the shell close to the universal horizon and we find which trivially integrates to where C 2 is an integration constant.
Collecting both solutions together we thus find that the reflected wave takes the final form where we have redefined the integration constants, absorbing factors of ξ,c, Ω and Λ into them.
As in the relativistic case, we find that the reflected wave gets exponentially red-shifted by a quantity controlled by the position of the universal horizon r uh = µ with an inverse power law. However, the numerical factor does not coincide with the one derived in General Relativity, although both are O(1). In the next section, we will discuss this numerical difference, tracing it back to the form of the high energy correction dominating the dispersion relation for large q.
The final step is to project this reflected wave onto the eigenstates of positive and negative energy in I + .
As shown in figure 2, for large radii we only have two possible solutions for the dispersion relation, which take a form analogous to the relativistic case Owing to this fact, we can straightforwardly borrow the result from the relativistic case and read the universal horizon temperature from our result to be Note however that in this caseκ = 2 3µ does not correspond to the surface gravity at the universal horizon, unlike what happens in a relativistic setting, where the surface gravity of the Killing horizon controls the temperature. This can be easily confirmed by computing κ uh following the prescription in [14,25] (appropriately modified for our different signature here) With this knowledge we can write our result for theκ governing Hawking's radiation as κ = 4κ uh 3 . (67)

V. OVERLAPPING OF WKB MODES NEARBY THE UNIVERSAL HORIZON
We offer here an alternative derivation of the temperature of the Hawking radiation just obtained in the previous section. We will follow the steps of [24], aiming to describe the WKB modes of out-going radiation close to the collapsing shell. These modes correspond to the natural plane waves in this region of the space-time. Afterwards, we will compute their inner product with a plane wave annihilating the Minkowsky vacuum inside the shell. By taking the late-time limit, as before, this will give an asymptotically good approximation to the overlapping of modes once the universal horizon has formed, allowing us to derive the thermal properties of the radiation, if any.
The main difference between [24] and our work here will be the choice of frame to construct WKB modes of the radiation close to the universal horizon. Since WKB modes are not diffeomorphism invariant, and this theory is not either, there is no reason as to why modes constructed in different frames must lead to the same Bogoliubov coefficients once projected onto flat waves. This of course opens a conundrum, the choice of frame for describing the modes corresponds to different observers moving close to the universal horizon with different kinematics. Since, as we have seen, this property is crucial for determining the thermal flux of the universal horizon, one can a priori expect the result to depend strongly on the chart of coordinates used when approaching the universal horizon. In [24] modes described in the aether frame were considered, thus depending on the preferred coordinates in this frame, which can be extended towards the inner region. Here instead, we choose to follow the generalization of free-falling observers to the case of a modified dispersion relation (13). Thus, we need to describe our WKB modes in the proper frame of these observers.
We start by performing a change of coordinates from dt to dū = dt − ζ + dr in order to adapt the coordinate system to the propagation of modes outgoing from the universal horizon. Although the expression forū might be complicated enough that it cannot be obtained analytically, we only need to perform the change on the metric and the aether vector, which depend only on r and the differentials.
We now go back to the decomposition of the momentum vector (12) and contract it with the generator of the spatial radial coordinate ∂ r in order to get the spatial momentum P in this frame. We thus have As previously discussed, we are interested in the shape of outgoing radiation. Note however that in this case we also need to take into account the branch of the solutions with negative energy, since we need to consider the whole vector space that the four modes span.
We now construct WKB modes in this region from the expression where c(r) and P (r) must be plugged case by case for every mode. Here we use Ω to denote the killing energy of these modes since this will be the base onto which we are projecting.
Going back to figure 2 we remind that there exist four possible modes propagating close to the universal horizon, which correspond to the soft modes, say η ± , and the hard modes, which we shall call Ψ ± . The ± denotes whether they correspond to the blue (orange) curve, identifying positive (negative) energy in the aether frame. Thus, in the region outside the universal horizon, any perturbation φ(u, r) with Killing frequency Ω can be decomposed in a basis of these modes where the ± denotes if they correspond to the branch with positive or negative energy.
As before, for those modes annihilating the Minkowsky vacuum, we take plane waves with Killing frequency Ω which we interpret here as the modes in the flat interior of the shell. These are the modes that we want to compare with those on the other side of the shell.
In order to evaluate the late-time behavior of the Bogoluvob coefficients, we can take into account that the shell is a closed surface, which implies that only modes with a certain frequency, inversely proportional to the radius R(s) of the shell, are possible. Thus, asymptotically late times can be identified onto the limit Ω → ∞. Finally, since we want to evaluate them in the exterior of the shell, we need to perform a change of variables involving the matching function α(u) = U (u) in (58).
The causal inner product of φ M (u) with η ± can be easily shown to vanish. Close to the universal horizon, the modes are regular and a standard dispersion relation |q| = Ω is recovered. Therefore, only soft positive modes overlap between themselves, leading to no radiation. We are thus left with the task of evaluating the overlap with the hard modes Ψ ± . Their momentum P ± (r) can be easily computed perturbatively in Λ to get where the sub-leading term distinguishes the two ± roots and we need to keep it in order for the derivative ∂ Ω P ± not to vanish. Note that the pre-factor that controlsκ appears in front of the term (r − µ) −1 which becomes a logarithm upon integration to construct the WKB mode. Indeed, if we compare this with q Λ in (26) we see that no term with a (r − µ) −1 arises for the momentum in the aether frame. This hints that the logarithm must take an important role in providing the thermal spectrum here.
We thus have two WKB modes and we focus on computing their inner product with the flat modes of positive energy φ M in the proximity of the universal horizon where Π ± = U µ ∂ µ Ψ ± and π M = U µ ∂ µ φ M are the momenta associated to the modes, derived from the action (10).
At the leading order the integration measure is just so that we get γϵ γR γI γ∞ Re x Im x Figure 4: Integration contour γ in the complex plane. The outer circumference has radius |x| → ∞, while the inner one is described by the curve x = e iθ , with 0 < θ < π/2. The branch cut is shown in blue.
where we have enclosed the different pre-factors into the constant A.
The computation of this integral is not trivial, due to the branch cut of both the square root and the logarithm. One of the integration limits lays on the branch point and therefore the exponent within the integrand is not a holomorphic function along the integration regime, which forbids us from applying a saddle point approximation. One could in principle try to regularize this by shifting the integration regime to {µ+ , ∞} and then taking the limit of to vanish. However, this fails, since the integral does not converge for real values of Λ and µ. Instead, an approximate formula which is enough for our purposes here can be obtained by deforming the integration contour.
First, we shift the integration variable to x = r − µ so that we attach the branch point to x = 0 and we choose to define the integrand through its principal value, with the branch cut running along the negative real axis in the complex plane. We then consider the integral where the contour γ = γ R + γ ∞ + γ I + γ is depicted in figure 4.
Using Cauchy's integral theorem, we find that ∆ ± = 0. On the other hand, the integral over γ ∞ can be shown easily to vanish, while that over γ does too if we close the contour through the upper half of the complex plane, as shown in figure 4. We are thus left with which implies We now perform another change of variables x = iy, so that we get Note that now this integral is more convergent that the original one when → 0, due to the real contribution to the exponential from the square root. Indeed, it now converges for positive as long as both Λ and µ are positive. Using the analytic extension of the Gamma function we get where the proportionality factor is the same for both Ψ ± .
Therefore we have that the ratio between the Bogoluvob coefficients becomes where the dots stand for sub-leading terms in r − µ.
So finally, in strict analogy with the derivation of Eq. (44), we can use the above result and the normalization condition for the Bogoliubov coefficients to deduce that the spectrum observed at I + will be Planckian with temperature T uh = 1/(3πµ) =κ/(2π). Noticeably, this result agrees with expression (65) obtained from the alternative derivation based on mode-matching for a collapsing shell.
An obvious question arises here. Why are we finding a thermal spectrum coming from the universal horizon while the authors of [24] did not find such an effect? In our language here the reason is clear. The temperature of the radiation is controlled by the (r − µ) −1 term in P ± (r), which integrates to a logarithm in the exponent of the WKB mode. This term only appears once we change frame and construct these WKB modes in the frame co-moving with free-falling observers or physical rays of the field φ. This choice was motivated by the fact that the quantum modes of the field correspond indeed to avatars of physical rays and therefore this seems to be the natural frame to follow their dynamics on this space-time. However, this choice is of course not the unique one. Instead, one could decide to stick to the preferred frame given by the aether time and its orthogonal hyper-surfaces, which corresponds to the choice made by [24]. In that case 9 , however, note that the momenta that we need to insert in order to construct the WKB modes are the q Λ ± , which do not have any (r − µ) −1 term. Instead, around the universal horizon, these WKB modes would collapse onto flat modes, since from which is trivial to check that no radiation is obtained.
In summary, we can say that the difference between our framework and that of [24] is tantamount to a different choice of the vacuum state at the universal horizon which, due to the non-relativistic causality structure of the theory, is not fully fixed by the choice of vacuum state just on I − . In our case, the state is supposed to be the vacuum for physical rays of the field, while in the case of [24] it was chosen to be that for the aether frame. Of course, we expect that if one would be able to set up a well posed Cauchy problem for a collapsing geometry -with the vacuum defined on some early times constant-khronon slice (on which static and freely follow observers basically coincide) -then the evolution of the problem will uniquely determine the global vacuum state defined on this geometry.

VI. UV-COMPLETING THE DISPERSION RELATION
A word of caution must be shed here with respect to the derivation of Hawking radiation that we have performed in this work. Going back to our result, we found that the thermal spectrum of Hawking radiation emitted by a universal horizon was controlled by the peeling of rays close to it through where κ U H = 1 2µ is the surface gravity of the universal horizon. The numerical coefficient in this result can be traced back to the behavior of the blue-shifted modes of the field φ departing close to the universal horizon, whose trajectory obeys the differential equation However, one must also note that the derivation in this paper was done perturbatively in large Λ, by expanding the dispersion relation (13) to leading order where therefore Λ takes the role of an UV cut-off for our computation.
Immediately, a question thus arises, since the momenta of the out-going modes, ultimately responsible for the production of Hawking radiation, become of the order of the cut-off close to the universal horizon. Indeed, So, although the perturbative expansion of the mode is under control, with sub-leading terms decaying appropriately with inverse powers of the cut-off and not hitting strong coupling at any point close to the universal horizon (where the expansion is valid), it is reasonable to question if we could be actually out of the validity regime of our perturbative approximation. Indeed, if instead of cutting the expansion of the dispersion relation at q 3 we retain the next suppressed correction we do find that our result forκ changes and becomes with a different numerical pre-factor. Actually, if we cut the series at an order q N Λ N −1 we find This agrees with the results obtained in [16,20,21] 10 using other methods and confirms our suspicion that our result for the thermal spectrum actually lays outside the range of validity of the perturbative expansion on Λ, and it is actually controlled by the higher order polynomial term in the expansion for the momentum. However, there is hope.
Instead of expanding the dispersion relation, and once we have understood how the mechanism of production of Hawking radiation around a universal horizon works, let us go back to the full dispersion relation (13).
Getting an analytical result with the full formula is hard because it requires to solve quartic equations that quickly become cumbersome and intractable. However, we do know now that Hawking radiation is solely controlled by the high energy behavior of the dispersion relation. We thus take, instead of the previous perturbative expansion, the opposite limit Λ → 0, for which thus focusing on the high-energy limit of the dispersion relation. Derivation ofκ is now immediate from (89) and leads to 11κ which is the result that one would have actually expected from the beginning if the whole computation that we have performed here was never attempted! Let us thus collect what we have learned here. If we were dealing with some unknown theory with dispersion relation given in a polynomial form then the derivation that we have performed along this work applies exactly and one concludes that universal horizons emit Hawking Radiation with a spectrum controlled byκ as given by (89). However, if instead, we are dealing with a dispersion relation which is an analytic function -but not a polynomial -we have to focus on the limit of large q of such dispersion relation, controlling the high energy behavior of the quantum modes in the theory, and again read the value ofκ from our result (89). Using this we find that for the case of the field with an action quartic in space derivatives (i.e. N = 2),κ exactly corresponds to the surface gravity of the universal horizon, eq. (91).
The fact that Hawking's temperature depends on the specific form of the dispersion relation for different matter fields opens up an intriguing possibility, as discussed in detail in [7]. Imagine that we have two different field species, that we denote as φ 1 and φ 2 , whose dispersion relation at very high energies takes the form with N 2 > N 1 . Then, the corresponding Hawking radiation emitted by a black hole of the kind that we have studied in quanta of the fields φ 1 and φ 2 will have different temperature which implies T 2 > T 1 . Now let us add two thin shells surrounding statically the black hole 12 with temperatures T A and T B satisfying The key property of these shells will be that the shell B only interacts with the matter particles φ 2 , while the shell A only does so with particles of φ 1 . Then, this implies that the system will try to move towards thermodynamic equilibrium by radiating from A onto the black hole using quanta of φ 1 -since T A > T 1and, at the same time, from the black hole onto B using quanta of φ 2 . For an observer sitting at infinity and observing the system, the net effect of this heat transfer will be a flux from A to B mediated by the black hole which, since T B > T A , violates the second law of thermodynamics and allows for the construction of a mobile of the second kind.
Violation of the second law of thermodynamics is, of course, an undesirable property which in our case seems to hint towards the need for a universal dispersion relation for all matter species coupled to Einstein-Aether gravity. Indeed, if in the previous case N 1 = N 2 , then all kinds of radiation have the same Hawking temperature and the violation of the second law becomes impossible. This is what happens, for instance, in Hořava Gravity, where all matter fields inherit the same universal momentum dependence on their dispersion relation due to matter-gravity interactions so that where the coefficients c i are species-dependent and d is the number of spatial dimensions, thus d = 3 for a 3 + 1 dimensional space-time. This implies that generically we should expect a universal spectrum with κ = 4/3κ uh (species-dependent constants like c 2d being always removable by a simple redefinition of the cut-off scale Λ which does not appear in the final expression for the temperature).

VII. DISCUSSION AND CONCLUSIONS
In summary, we have shown in two alternative ways that the Hawking radiation characterizing a black hole endowed with a universal horizon is set by the effective surface gravity characterizing the peeling of physical null rays in its proximity Eq. (89). This result implies that the temperature generically is not just set by the surface gravity of the universal horizon Eq. (66) but depends on the high energy behaviour of the dispersion relation of the specific field considered. Robustness of the black hole thermodynamics requires then that such UV behaviour of the dispersion relation should be universal (to avoid different temperatures of different fields on the same black hole space-time), like for example the one expected in Hořava gravity. While we deem these results a definite progress in understanding the internal consistency of these theories and the basic building blocks of black hole thermodynamics, there are nonetheless some open issues that still need to be addressed.
Although we have explicitly shown that Hawking temperature is determined by the UV behaviour of the field dispersion relation, one has to be aware that the black hole solutions with universal horizon have been so far derived only as solutions of the low energy gravitational Lagrangian (for example in Hořava gravity the so called L 2 part, providing just second order field equations for the gravitational field). It is not known if the UV completion of the gravitational Lagrangian will still admit such black holes as solutions, nor if the modified solution will still admit universal horizons. However, an examination of curvature invariants on the universal horizon, which can be arbitrarily small, seems to indicate that one can always find a situation in which the UV-completing terms -given in terms of higher powers of curvature invariants, at least in Hořava Gravity -can be always taken as small perturbations. If this is the case, it seems plausible that the universal horizon will survive in the full theory with similar properties to the ones described here. It is true nonetheless that a full account of its dynamics requires the knowledge of the UV-complete theory up to arbitrary energy.
Another point worth further investigating is related to the difference between our computation based on the WKB approximation and the one reported in [24]. As we have seen, the two derivations lead to different conclusions just because of a different choice of the vacuum state. In our derivation, we assumed that the regular state on the black hole space-time is the analogue of the standard Unruh state for physical free falling observers (which do not follow the geodesics of the metric but rather those determined with the metric-aether background). In the case of [24] a vacuum state defined with respect to the aether frame was chosen instead. Which of the two choices is the physically relevant one, cannot be settled just by the concordance of the second derivation with the first one, based on matching regular sets of coordinates along a physical ray trajectory, as the latter implicitly assumes the regularity of the vacuum state all along these worldlines.
In relativistic derivations, it can be shown exactly that the Unruh state is the only regular state compatible at late times with a collapsing geometry leading to the formation of a horizon (see e.g. [47,48]) 13 . However, some ambiguity arises in non-relativistic theories as no infinite accumulation of physical rays happens anymore at the Killing horizon (due to the modified dispersion relation) and hence no infinities are expected there in the energy distributions in any frame. Considering non-relativistic theories of gravity as well, led to the discovery of universal horizons which may remove such ambiguity by restoring the peeling behaviour typical of relativistic theories and thus reintroducing the key element for the singular behaviour of global quantum states which are not vacuum at (universal) horizon crossing. Nonetheless, an ambiguity remains due to the fact that such regularity can be defined w.r.t. alternative frames: the one for physical freely falling observers and the one for observers falling with the aether frame. In this work, we have deemed the first one as the physically relevant frame, but only an eventual calculation of the renormalized stress energy tensor analogue to the one carried out in relativistic physics will be able to address the uniqueness and regularity of the vacuum state on a collapsing geometry in a non-relativistic framework.
In order to do this, several road blocks will have to be removed, for example, the study of the Hadamard condition for Green functions (see e.g. [43]) in non-relativistic settings, as well as the generalization to spacetimes with universal horizons of the theorems concerning vacuum state regularity in globally hyperbolic space-times [47,49,50]. Let us stress that such generalizations are far from clear in the quite contrived causality structure of such geometries [51] and that the possibility to fully characterize the vacuum state is not guaranteed due to the peculiar nature of universal horizons that make them similar to Cauchy horizon in non-relativistic frameworks [10,51]. 14 We hope to be able to tackle these issues and more in forthcoming investigations.