Analytical theory of Hawking radiation in dispersive media

Hawking's 1974 prediction that black holes radiate and evaporate has been hinting at a hidden connection between general relativity, quantum mechanics and thermodynamics. Recently, laboratory analogues of the event horizon have reached the level where tests of Hawking's idea are possible. In this paper we show how to go beyond Hawking's theory in such laboratory analogues in a way that is experimentally testable.

Waves in moving media experience the analogue of the event horizon when the velocity of the medium exceeds the speed of the waves. Consider waves propagating with phase velocity c in a medium moving with local velocity u. Here c depends on the wavenumber k according to the dispersion relation in the material-in a co-moving frame-while u varies in space according to the velocity profile of the medium. We can always decompose a wave into its Fourier components φ with respect to the (angular) frequency ω. Near the horizon, we can assume a one-dimensional (1D) model [4] where both φ and u depend on the spatial coordinate z only. The medium shall move to the left, i.e. with negative u, forming just one black-hole horizon, for simplicity. In a co-moving frame, the frequency is Doppler-shifted to ω − uk and equals ck for waves propagating against the current. We do not discuss waves propagating with the flow, as they are not affected by the horizon, and we assume that waves propagating with and against the current are decoupled. The wave propagation in moving media is thus described by the dispersion relation For isotropic media the phase velocity c is an even function of k, so c(k)k is an odd function. Figure 1 illustrates the dispersion relation (1) and shows how to find solutions graphically. These solutions characterize the waves in the moving medium. Consider now localized wave packets with position z, wavenumber k and carrier frequency ω. Figure 2 illustrates the typical trajectories of such wave packets. Their dynamics are governed by Hamilton's equations dz dt = ∂ω ∂k = v g + u, Here u = du/dz denotes the velocity gradient of the medium and v g = d(ck)/dk the group velocity of the wave, which, in dispersive media, differs from the phase velocity c. The horizon is formed where the flow exceeds the wave velocity, but we need to identify which velocity we mean. Here we assume a group-velocity horizon where |u| reaches v g . Yet Hawking radiation also relies on the possibility of exciting waves for which |u| exceeds c. Figure 1 shows [19] that three types of waves may exist in a region where |u| < v g : a right-moving wave φ ur that escapes from the horizon, a left-moving wave φ ul that drifts towards the horizon and a leftmoving wave φ u with negative wavenumber k u [10]. We see from the dispersion relation (1) that k = (c + u) −1 ω, so k is negative when |u| > c; for φ u the medium exceeds the phase velocity. (A fourth wave φ v propagates with the flow and is ignored in our analysis.) On the other side of the horizon, where |u| > v g , only the negative-k wave φ u continues to propagate (apart from φ v of course). In the astrophysical black hole, φ u would be on its way to the singularity; it constitutes the Hawking partner of the escaping wave. Figure 2 indicates that incident waves with positive k are partially converted at the horizon into their Hawking partners with negative k and vice versa.  (1) for waves in moving media. Waves may propagate with wavenumbers k j given by the intersections of the dispersion curve ck (blue) with the straight lines (red) of the Doppler-shifted frequency ω − uk depending on the velocity u of the medium (u < 0). The group velocity v g is given by the slope of the dispersion curve ck. For a wave with k ur the group velocity exceeds −u, the wave thus propagates to the right, away from the horizon. For k ul the wave is left moving with v g + u < 0. Where k ur and k ul coincide the flow velocity reaches the group velocity, forming a group-velocity horizon. A third solution k u may exist that corresponds to a Hawking wave with negative wavenumber that is left-moving. The fourth solution k v describes waves propagating with the current that do not contribute to the Hawking effect. We used the dispersion relation c 2 = c 2 0 (1 − k 2 /k 2 0 ) with the constants c 0 and k 0 and plotted k and ω in units of k 0 and ω 0 = c 0 k 0 , respectively. This dispersion relation is said to be subluminal (or normal), because the phase velocity decreases with increasing k; it is the simplest subluminal dispersion relation. We could also chose superluminal (or anomalous) dispersion, for example c 2 = c 2 0 (1 + k 2 /k 2 0 ), where c increases with k. It turns out that in this case we obtain similar Hawking temperatures. Figure 2 also characterizes the in-going and out-going modes. As the mode conversion is linear, an out-going mode is the linear superposition of two in-going modes: The modes are orthogonal and normalized to δ(ω 1 − ω 2 ) in a sense made precise in the scalar product (A.2) defined in appendix A. The negative-k mode φ in u turns out to have the negative norm −δ(ω 1 − ω 2 ). From this follows an important conclusion that does not rely on the details of the scalar product (φ 1 , φ 2 ) but only on the fact that (φ 1 , φ 2 ) is bi-linear in φ * 1 and φ 2 . Calculating the norm of the superposition (3) we obtain In the diagram for the in-going mode φ in ul , an incident wave with wavenumber k ul is reflected at the horizon and leaves with wavenumber k ur , while some of the incident amplitude is converted into a wave with negative k u that crosses the horizon (dashed line). The other in-going mode φ in u is incident with negative k u and continues across the horizon, sending off a k ur to the right (dashed line). The φ out ur and φ out u describe the out-going modes that leave with wavenumbers k ur and k u and originate from φ ul and φ u contributions. We used the velocity profile u(z) = (u R + u L )/2 + tanh(z/a)(u R − u L )/2 with constant asymptotic velocities u R and u L and constant scale a and solved Hamilton's equations (2) for the dispersion relation of figure 1 (u R = −0.8c 0 , u L = −1.2c 0 , a = 1/k 0 , ω = 0.035ω 0 ); z and t are plotted in units of 1/k 0 and 1/ω 0 , respectively. which implies that the amplitude of the incident mode φ in ul is multiplied by a factor |α| larger than 1; the horizon thus acts as an amplifier [30,31]. The negative-k mode φ in u represents the reservoir of the amplification noise [30].
According to quantum mechanics [32], amplifiers generate noise even when the reservoir is empty; they create real particles from the quantum vacuum in correlated pairs [30]. In the case of a horizon, these particles constitute the Hawking radiation; the particles of the φ out ur are escaping while their Hawking partners of φ out u are swept away. Correlated particles are most strongly entangled; if we only consider the escaping radiation, averaged over the Hawking partners, the reduced quantum state of each escaping mode is therefore a state of maximal entropy [33], a thermal state. Its temperature T is given by the Boltzmann law [30] β α whereh denotes Planck's constant divided by 2π and k B is Boltzmann's constant. In the following we calculate T . Here it is wise to characterize wave packets by their wavenumbers k instead of their positions z, because wavenumber is a better indicator of a horizon than position: there is a clear gap between positive and negative k, but not in z: negative-k waves may freely cross the horizon exploring the entire physical space (see φ in ul in figure 2). Therefore, we rather represent z as a function of k, identifying the present position of the wave packet by its present wavenumber. Mathematically, we simply solve the dispersion relation (1) for z with k as variable.
The Hawking effect bridges positive and negative wavenumbers. Figure 3 visualizes this connection on the complex plane where we regard k as a complex variable. Assuming a slowly varying velocity profile, we derive in appendix A a simple rule for obtaining the Hawking temperature (5) graphically, which goes as follows. The red lines of figure 3 indicate the physically allowed wavenumbers. The relative phase between the positive and negative k wave is given by the contour integral of z from positive to negative wavenumber, avoiding the branch cut between them. This phase is not a real number but contains an imaginary part. If we close the contour (black ellipse in figure 3) the real part of the integral vanishes and twice the imaginary part remains. Now, the relative amplitude |α/β| between negative and positive component is the exponential of the imaginary part of the phase, and so |α/β| 2 must be the exponential of twice the imaginary part, i.e. of the closed contour integral (black ellipse in figure 3). We thus obtain from formula (5) our result hω Note that we can also represent z dk as − k dz, but the corresponding contour in the complex z-plane can be more complicated.
Formula (6) is the main result of this paper. It extends Corley's analytical theory [23] to arbitrary velocity profiles (beyond linearization at the horizon) and generalizes the interpretation of Hawking radiation as tunnelling [34][35][36] to dispersive media. We easily reproduce Hawking's formula in the regime of weak dispersion where v g changes little with wavenumber. In this case, we see from Hamilton's equations (2) that dz/dt remains nearly zero, while the wavenumber k falls exponentially with the rate u that remains nearly constant, because z does not change much. Differentiating the dispersion relation (1) with respect to ω we get dz/dω = 1/(u k) and so d dω which gives 2π/u since u remains constant over an exponentially large range of k. After integrating over ω we arrive at Hawking's formula [1] in moving media [5] For strong dispersion, however, we cannot regard u (z) as almost constant in the contour integral (7). Figures 4 and 5 show that in this case T varies with frequency: the Hawking  1), but they are connected for complex k. Therefore it is wise to represent the solution z(k) of the dispersion relation (1) on the complex k plane; here we illustrate the imaginary part of z using different shades of blue (for the dispersion relation of figure 1 and the velocity profile in figure 2). At the red lines Im(z) = 0, so z is real there; they are drawn along the values of k that appear during the real wave propagation. Lines with real z(k) are disconnected by branch cuts where the imaginary part of z abruptly changes (the shades of blue are discontinuous here). The solid blue line connects the wavenumbers (blue dots) that are physically allowed (figure 1) for the out-going wave φ out ur on the far right side of the horizon (figure 2); the dotted blue line corresponds to a point further to the left, beyond the groupvelocity horizon, where the wave is exponentially damped. The contour integral along the blue curve gives the logarithm of the relative amplitude and the phase of the wave components of φ out ur ; the contour is freely deformable where z(k) is continuous. The closed-contour integral (black ellipse) gives the Hawking temperature according to formula (6).
temperature (8) is no longer universal and, as a consequence, the Hawking spectrum is non-Planckian. Formula (6) still predicts that T is inversely proportional to the scale a of the velocity profile. The Hawking temperature remains proportional to 1/a until for large velocity gradients  For slowly-varying velocity gradients the contour integral (6) describes Hawking radiation very well, including deviations from Hawking's formula (8) that have only been seen as an approximation [24] or numerically [27,28] so far. Appendix B describes three examples where our theory is exactly solvable, including a case closely related to the Schwarzschild black hole. The deviations from formula (8) are remarkably small, except when the velocity profile is asymmetric ( figure 5). However, as they enter the mode-conversion rate (5) exponentially, they are observable in laboratory analogues within the validity range of our analytical method, see appendix C. What lies beyond Hawking's theory is now at the horizon, coming into view.

Acknowledgments
We are grateful to Simon Horsley, Renaud Parentani and Thomas Philbin for valuable discussions. In particular we thank Renaud Parentani for pointing out that the negativewavenumber component in the out-going mode gives directly β/α. Our work is supported by EPSRC and the Royal Society.

Appendix A.
In this appendix we describe the main method for obtaining our analytical results, the phaseintegral method in k-space [19]. Corley [23] used phase integrals in position space, but horizons are better defined in k-space; Corley's method is restricted to linear velocity profiles, whereas our method works for arbitrary u(z).
We begin with an analysis of the classical wave equation in moving media in position space [18]: For solutions φ 1 and φ 2 of equation (A.1) we define the scalar product [5] ( that is a conserved quantity [18] related to the number of particles associated with a given wave. In a region of uniform flow velocity the solutions are the plane waves for the roots k j of the dispersion relation (1) shown in figure 1. We have normalized the plane waves (A.3) to a delta function in ω. Note that φ ur and φ ul have positive norm, while the norm of φ u is negative. When u(z) varies, asymptotically plane waves (A.3) are partially converted into each other. Here it is advantageous [19] to consider the spatial Fourier transformφ(k). We derive from equation (A.1) the wave equation in k-space: For slowly varying velocity profiles we approximateφ using a phase-integral method in k-space. For this, we regard u(i∂ k ) as u(i ∂ k ) where is a small parameter, we represent kφ as exp[ −1 (φ 0 + φ 1 + 2φ 2 + · · ·)] and substitute this ansatz in equation (A.4). Sorting the result into powers of produces a coupled system of equations for theφ m that we truncate at m = 1. Finally we remove , incorporating it in the scale of u(z), by formally setting = 1. In this way we obtain the square of the dispersion relation (1) with i∂ kφ0 = z and 2∂ kφ1 + ∂ k ln[u (i∂ kφ0 )c(k 2 )] = 0, which gives in k-spacẽ To deduce φ in z-space, we evaluate the inverse Fourier transformation in saddle-point approximation: where, for 0 = 1/(4π 2 ), the φ j are the expressions (A.3) for the normalized plane waves with the solutions k j of the dispersion relation (1), but now z may vary. The phases and amplitudes of the components φ j may depend on the contours of the phase integral in equation (A.6), although not on their form, but on their topology. Choosing topologically different contours [23] we obtain the various incident and out-going waves of figure 2. Figure 3 shows the contour that distinguishes the out-going mode φ out ur of the escaping Hawking radiation defined in figure 2. The contour is chosen by the following argument. On the left side of the group-velocity horizon, φ out ur decays exponentially; therefore the only physically allowed wavenumber is k ul with negative imaginary part, which in figure 3 lies on the lower half of the red 'cross bow'. The contour for φ out ur (dotted curve) comes in from ∞ in a quadrant where z(k) has a negative imaginary part such that the phase integral in expression (A.5) vanishes there. On the right side of the horizon, we have the three saddle points k ul , k ur and k u shown as blue dots. They correspond to the three partial waves (figure 2) of φ out ur . The contour for φ out ur for the right side of the horizon must be a deformation of the contour for the left side and pass through the physically relevant saddle points, which determines the contour shown in figure 3.
The out-going mode φ out ur is the superposition (3) of the in-going modes φ in ul and φ in u . On the right side of the horizon, the k ul -component comes from α φ in ul and the k u -component from β φ in u . The k ul -component can only differ from the k ur -component by a phase factor, because the z(k) integral between k ur and k ul is real on the real axis and hence real for any contour. We conclude that |β/α| corresponds to the relative weight of the k u -component in φ out ur that we can read off from our result (A.6). One advantage of this method is that we infer |β/α| from partial waves with high wavenumbers where the saddle-point approximation required for deriving expression (A.6) is sufficiently accurate. According to expression (A.6) the amplitude of the k u -component is given by exp(−Im z dk) taking the lower half of the integration contour in figure 3. Since the velocity profile is real for real z, its analytic continuation on the complex plane obeys u(z * ) = u * (z), and the same must be true for z(k). Therefore, although we integrate from the negative to the positive wavenumber, we may close the contour on the upper half of the complex k-plane. The real part of the integral then vanishes, while the imaginary part is doubled. Exponentiating gives |β/α| 2 and thus the formula from which our result (6) follows: In this appendix we consider three simple examples where our formula (6) admits exact solutions for dispersion relations where c(k) is an analytic function in k: the linear velocity profile, the z −1 profile and the z −1/2 profile. The examples show how our theory is to be applied. Furthermore, the z −1/2 profile corresponds to a prominent case: the Schwarzschild black hole [8]. Formula (6) requires z as a function of k for a given (angular) frequency ω. We obtain from the dispersion relation (1): To obtain z(k) we thus need to invert the function u(z) and substitute ω/k − c(k) for u. Consider first a linear velocity profile, where we get Assuming c(k) to be an analytic function, the only contribution to the contour integral (6) is the ω/(αk) term, which gives Hawking's result (8). We see that the Hawking radiation caused by a linear velocity profile does not depend on c(k): it is insensitive to dispersion. The wavenumber-dependance of c can only play a role in nonlinear profiles. Consider a situation where u is inversely proportional to z. Although this is an artificial case, it has the advantage that the entire wave propagation is exactly solvable for analytic c(k) [26]. Let us see whether our theory reproduces the known result [26] for Hawking radiation. Assuming This expression has single poles at the zeros of the denominator that correspond to the solutions of the dispersion relation (1) for z → ∞ where the velocity (B.5) of the medium vanishes, i.e.
As the integration contour (figure 3) encloses the k ur of the outgoing wave, we take the k ∞ that corresponds to k ur at ∞. We evaluate formula (6) using Cauchy's residue theorem where we require expression (B.6) for k near the pole k ∞ . Here where v g denotes the group velocity. We obtain from the residue theorem which is the previously known exact result for the 1/z profile [26]. It shows that the dispersion influences the spectrum of Hawking radiation; it is no longer a Planck spectrum. We can cast expression (B.9) in a form where we immediately see how it reproduces Hawking's result for zero dispersion. Let us define the phase horizon z ∞ as the position where the medium reaches the phase velocity at k ∞ : For the profile (B.5) the velocity gradient at the phase horizon is and so we obtain from the result (B.9) In the absence of dispersion the group velocity agrees with the phase velocity, which gives Hawking's result (8).
Finally, consider a velocity profile that closely corresponds to the Schwarzschild black hole. In Painlevé-Gullstrand coordinates the Schwarzschild solution appears like a medium moving with velocity profile u = −c (a/r ) −1/2 [8], where r is the radius (the distance to the singularity), a the Schwarzschild radius and c the speed of light. Considering only the wave propagation in r direction and identifying r with z we assume u = −u 0 a/z. (B.13) In this case, (B.14) Similar to the profile (B.5) the function z(k) has poles at k ∞ that are solutions (B.7) of the dispersion relation (1) at z → ∞. As before, we take the k ∞ that corresponds to k ur at ∞. However, due to the square in expression (B.14), z(k) has both a double and a single pole at k ∞ . Only the single pole contributes to the contour integral (6). We extract the single pole and obtain from Cauchy's residue theorem As in the case of the z −1 profile we calculate the velocity gradient α at the phase horizon (B.10), where we get for the z −1/2 profile (B.13) In practice [16], this is a white-hole horizon. The white hole is simply the timereverse of the black hole. In the aquatic analogue of the white hole the magnitude of the flow is rising at the horizon from sub-to super-wave speed (at the blackhole horizon the flow is falling from super-to sub-wave speed). An incident wave packet (top) is sent against the rising current. It propagates until the flow speed reaches the group velocity of the wave, whereupon it is converted into two wave packets (bottom), one with positive wavenumber (bottom, right) and one with negative wavenumber (bottom, left).
In terms of α we obtain from expression (B.15) and the dispersion relation (B.7) the result This formula describes how the Hawking radiation of the black hole is modified if space were dispersive [3]. For instance, if space were discrete at the Planck scale, this discreteness would first appear as a wavenumber dependance of the speed of light [30], i.e. as dispersion. Our formula (B.17) would describe the non-Planckian spectrum of Hawking radiation due to the Planck scale.

Appendix C.
In this appendix we apply our theory to an example of direct experimental relevance, the observation of stimulated Hawking radiation with water waves [16]. We show that deviations from Hawking's prediction of a thermal spectrum are observable in the regime where our analytical theory works best. Consider water waves in a channel of varying height h (figure C.1) where the water flows with a velocity profile u(z) that can be regulated by h and the initial speed of the water. Water waves-gravity waves-obey the dispersion relation [37] c(k) = g tanh(hk)/k, (C.1) the relative intensity |β| 2 of the negative-k waves is given by equations (4)-(6), or, in explicit form The contour integral connects positive with negative wavenumbers. We can choose any contour, as long as we do not cross branches, because the value of the integral is contour-independent. For our calculation we chose an ellipse between the positive and the negative k value for the asymptotic flow. Figure C.2 compares |β| 2 with the analogue of Hawking's prediction, formula (8) in equations (4) and (5) that gives the Planck distribution where u is taken at the horizon for the limit ω → 0. We see that the ratio |β| 2 /|β 0 | 2 can deviate significantly from 1, in particular for slow variations of the flow velocity over a large length scale a, which is the regime where our theory is valid. Therefore, deviations of the Hawking spectrum from the Planck distribution (C.6) seem observable in laboratory analogues of the event horizon.