Hawking radiation on a falling lattice

Hawking radiation on a falling lattice Ted Jacobsona,b∗ and David Mattinglya† arXiv:hep-th/9908099v2 12 Jan 2000 a Department of Physics, University of Maryland, College Park, MD 20742-4111, USA b Institute for Theoretical Physics, UCSB, Santa Barbara, CA 93106 [4,5]). Quantizing the sound field Unruh argued that the horizon would radiate thermal phonons, in analogy with the Hawking effect, at a temperature given by h̄/2π times the gradient of the flow velocity at the horizon. Since the sound field has an atomic cutoff, and the physics of the fluid is completely understood in principle, the fluid model may provide a source of insight into the Hawking process. The atomic cutoff produces dispersion of sound waves leading to subsonic propagation at high frequencies. It turns out that this dispersion is all that is needed to obviate the need for a trans-Bohrian reservoir. In the linearized theory, an outgoing short wavelength mode can be dragged in by the fluid flow and redshifted enough so that its increased group velocity overcomes the flow and it escapes to infinity. Investigations of a number of linear field theory models [6–11] and lattice models [12] incorporating such high frequency dispersion have demonstrated this mechanism of “mode conversion” and have shown that the Hawking radiation is recovered.2 The dispersion in these models is not locally Lorentz-invariant since a frame is picked out in which to specify which frequencies are “high”. In the fluid model this is the local rest frame of the fluid. A real black hole also defines a preferred frame but it does so in a non-local way. Perhaps quantum gravity produces a dispersion effect related to this non-local notion of a preferred frame, or perhaps microphysics is just not locally Lorentz invariant. In any case, it seems worthwhile to understand the mechanism of the outgoing modes and Hawking radiation in these models for the hints it may provide about a correct quantum gravity account. In [12] a falling lattice model was introduced to address the “stationarity puzzle” [17,8]: how can a low frequency outgoing mode arise from a high frequency ingoing mode when it propagates only in a stationary background spacetime and hence has a conserved Killing frequency? In the continuum based dispersive models there is no satisfactory resolution of this puzzle. As the outgoing modes are traced backwards in time their incoming progenitors are squeezed into the short distance regime Scalar field theory on a lattice falling freely into a 1+1 dimensional black hole is studied using both WKB and numerical approaches. The outgoing modes are shown to arise from incoming modes by a process analogous to a Bloch oscillation, with an admixture of negative frequency modes corresponding to the Hawking radiation. Numerical calculations show that the Hawking effect is reproduced to within 0.5% on a lattice whose proper spacing where the wavepacket turns around at the horizon is ∼ 0.08 in units where the surface gravity is 1. I. INTRODUCTION If nothing can emerge from a black hole then where does the Hawking radiation come from? Relativistic field theory tells us that it comes from a trans-Planckian reservoir of outgoing modes at the horizon [1]. The existence of such a trans-Planckian reservoir is doubtful on general physical grounds. It is the cause of the divergences in quantum field theory which are expected to be removed in a quantum theory of gravity, and in particular it would seem to imply an infinite black hole entropy unless canceled by an infinite negative “bare entropy”. String theory has produced an account of the Hawking radiation and black hole entropy which requires no trans-Planckian reservoir, however in that picture instead of a black hole one has a stringy object in flat spacetime that arises from a black hole as the string coupling constant tends to zero [2]. Thus gravitational redshift plays no role in the string calculations and nothing is learned about the origin of the outgoing modes in a black hole spacetime.1 If we could follow the theory as the string coupling becomes strong we would presumably learn how it is that string theory produces the outgoing black hole modes. Absent a quantum gravity description of the Hawking process, Unruh invented a fluid analog of a black hole [3] in which an inhomogeneous flow exceeds the long wavelength speed of sound creating a sonic horizon (see also 2 Dispersion can also produce superluminal propagation, which allows outgoing modes to emerge from behind the horizon. The superluminal case does not seem very healthy from a fundamental point of view, though it is surely relevant in some condensed matter analogues of black hole horizons [13–15] and is capable of producing the usual Hawking spectrum [16]. ∗ [email protected] † [email protected] 1 The fact that redshifting plays no role in producing the stringy Hawking radiation is presumably tied to the fact that the string calculations apply in the “near-extremal” limit that the surface gravity tends to zero. 1 horizon. The precise choice of worldlines for the lattice points should not be important to the leading order Hawking effect as long as the ground state evolves adiabatically in the lattice theory. Our numerical results are consistent with this expectation, in that the particular worldlines depend on the lattice spacing yet the results converge to the continuum as the lattice spacing is decreased. Moreover it was shown recently [10], in a continuum model with high frequency dispersion, that when the preferred frame is changed from the free-fall frame of the black hole to that of a conformally related metric there is no leading order change in the Hawking radiation unless the acceleration of the preferred frame is very large. Presumably all that matters is that the preferred world lines flow smoothly across the horizon. From a fundamental point of view, if there really is a cutoff in some preferred frame, it is plausible that this frame would coincide with the cosmic rest frame and would fall across the event horizons of any black holes that form by collapse. The remainder of this paper is organized as follows. In Section II the falling lattice model is set up, in Section III some results of the WKB method are shown, illustrating the frequency shift and the mode conversion at the horizon, and the role of the negative frequency branch is explained. In Section IV the results of the full numerical evolution are presented, and Section V contains a discussion or the results and directions for further work. The Appendix describes the finite differencing of the time derivates used in the numerical calculations. We use units with h̄ = c = κ = 1, where κ is the surface gravity. of the model which is not physically sensible. A lattice provides a simple model for imposing a physically sensible short distance cutoff. It might seem at first that the most natural choice would be to preserve the time translation symmetry of the spacetime with a static lattice whose points follow accelerated worldlines. On a static lattice, however, the Killing frequency is conserved, so outgoing modes arise from ingoing modes with the same frequency, and there is no Hawking radiation. The in-vacuum therefore evolves to a singular state at the horizon. If the lattice points are instead freely falling, a discrete remnant of time translation symmetry can still be preserved. On such a lattice there is still a frequency conservation law, and the unsatisfactory resolution of the stationarity puzzle is that the lattice spacing in this case goes to zero at infinity (if the lattice points are asymptotically at rest) [12]. If we insist that the lattice spacing asymptotically approaches a fixed constant and that the lattice points are at rest at infinity, we find a fascinating resolution of the stationarity puzzle [12]: the lattice cannot have even a discrete time translation symmetry, since there is a gradual spreading of the lattice points as they fall toward the horizon. The timescale of this spreading is of order 1/κ where κ is the surface gravity. This time dependence of the lattice is invisible to long wavelength modes which sense only the stationary background metric of the black hole, but it is quite apparent to modes with wavelengths of order the lattice spacing. On such a lattice the long wavelength outgoing modes come from short wavelength ingoing modes which start out with a high (lattice scale) Killing frequency which is exponentially redshifted as the wavepacket propagates to the horizon and turns around. This behavior of wavepackets in the (time-dependent) falling lattice model was determined in [12] with the help of the WKB approximation in which the frequency is used as a Hamiltonian to solve for the wavepacket trajectory. It was also argued that, since the time dependence of the lattice is adiabatic for those modes with wavelength short enough to know they are on a lattice, these modes will remain unexcited on their way to the horizon and the ground state condition for the Hawking effect will be met. However, since the wavelength is of order the lattice spacing and the wavepacket can vary wildly on the lattice near the horizon, it is not obvious that the WKB approximation is reliable nor that the adiabatic argument really applies. The primary motivation for the work reported here was to check the WKB and adiabatic assumptions by carrying out the exact calculation numerically using the lattice wave equation. We found that indeed the field behaves this way, and the Hawking radiation is recovered to within half a percent for a lattice spacing δ = 0.002/κ which corresponds to a proper spacing ∼ 0.08/κ where the wavepacket turns around at the horizon. We also studied the deviations from the Hawking effect that arise as κδ is increased, and our simulations reveal an interesting picture of how the wavepackets turn around at the II. FALLING LATTICE MODEL In this section we set up the field theory on a lattice falling into a black hole in two spacetime dimensions. A. Black hole spacetime We assume the spacetime metric is static, so [18] coordinates can be chosen (at least locally) such that the line element takes the form
2 ds2 = dt2 − dx − v(x)dt . (2.1) In these coordinates the Killing vector is given by χ = ∂t , (2.2) χ2 = 1 − v 2 (x). (2.3) with squared norm For v(x) we choose v(x) = − 2 vmax , cosh(βx) (2.4) where vmax and β are positive constants. As |x| → +∞, v(x) vanishes, so the line element becomes that of Minkowski spacetime. Provided vmax > 1 there is an ergoregion centered on x = 0 in which the Killing vector is spacelike, and the boundaries of this region are black and white hole horizons at xH = ±β −1 cosh−1 (vmax ). ∂t a/a = (2.13) whose maximum outside the √ horizon (cosh(βx) p > vmax ) 2 2 is κ (at the√horizon) if vmax ≥ 2 and κvmax /2 vmax −1 < if vmax ≤ 2. Thus |∂t a/a| ∼ O(κ) everywhere outside the horizon as long as vmax > 1 is not too close to unity. (2.5) The surface gravity of the horizons is given by p −2 . κ = v ′ (xH ) = β 1 − vmax βvmax sinh(βx) cosh2 (βx) B. Lattice (2.6) The lattice is defined by discretizing the z coordinate with spacing ∆z = δ. That is, we have only the discrete values zm = mδ. We will consider values of κδ which range from 0.002 to 0.128. To give a feeling for the lattice we plot in Fig. 1 the (t, x) coordinates of lattice point worldlines at intervals of 50 lattice points. The surface gravity sets the length scale for our spacetime, and we will keep it fixed. Hence it will be useful to employ units for which κ = 1. We now introduce a freely falling “Gaussian” coordinate z which will be discretized to define a falling lattice. The worldlines with dx = v(x)dt are geodesics of the metric (2.1) with proper time t, at rest as |x| → ∞, and they are orthogonal to the surfaces R of constant t. On these curves the quantity t − dx/v is constant, so if z is constant on these curves we must have R W (z(x, t)) = t − dx/v for some function W . We choose this function R z so that z = x at t = 0, which implies that W (z) = − dx/v = sinh(βz)/βvmax , hence sinh(βz) = sinh(βx) + βvmax t. 30 20 (2.7) In terms of z, the line element (2.1) takes the form ds2 = dt2 − a2 (z, t) dz 2 , (2.8) 10 (2.9) 0 with cosh(βz) . a(z, t) = p 1 + (sinh(βz) − βvmax t)2 We call a(z, t) the “scale factor”, even though it depends on z as well as t. It can also be expressed as v(x(z, t)) a(z, t) = . v(z) 0 v(x(z, t)) ∂z . a(z, t) 2 3 4 5 6 FIG. 1. Every 50th lattice point, t vs. x, for the case √ δ = 0.004 and vmax = 2. In the black region of the plot the points are too close to resolve. The horizon is at x ≃ 0.6. (2.10) In (t, z) coordinates the Killing vector (2.2) is given by χ = ∂t − v(z)∂z = ∂t − 1 Since dx = a(z, t)dz at constant t, the separation ∆x of the lattice points in Fig. 1 gives a direct representation of the scale factor. The proper lattice spacing is given by a(z, t)δ which, according to (2.12), grows approximately linearly with time at the horizon. Note that there is a region where the lattice spacing becomes very small. This happens because we chose the spacing to be uniform at t = 0. If the points did not bunch up outside the black hole before t = 0, they could not wind up equally spaced at t = 0, since they are all on unit energy free-fall trajectories. (The line element is symmetric under (t, x) → (−t, −x), so the points bunch up also outside the white hole horizon after t = 0.) To get a better idea of how sparse the lattice is near the horizon at the turning point of a typical wavepacket, (2.11) The scale factor is unity at t = 0, and outside the black hole it grows towards the future. It is important to know how rapidly the scale factor is changing in time, since in the lattice theory that time dependence is not merely a coordinate effect and can therefore excite the quantum field. At the horizon the scale factor takes the value q −2 a(zH , t) = 1 + 2κt + κ2 t2 /(1 − vmax ), (2.12) so a(zH , t) ∼ κt for κt > ∼ 1. Thus the fractional rate of change ∂t a/a is of order κ. In general, we have 3 we plot in Fig. 2 the velocity function v (2.4) for every lattice point near √ the horizon at time t = 27, for the case vmax = 2 and κδ = 0.004. In this case there are only 6 points in the region between v(x) = 0.5 and the√horizon, and the proper spacing at the horizon is ∼ 2κδt ∼ 0.15/κ. -1 -2 1 2 hence there is no wave scattering. The discretization breaks this conformal invariance, however there is still essentially no scattering on the lattice providing the lattice spacing is much smaller than the radius of curvature. The reason is that modes with wavelength long compared to the lattice spacing cannot tell they are not in the continuum, while modes with wavelength comparable to the lattice spacing cannot tell that the spacetime is curved. Thus, when we compute the Hawking occupation numbers, there is no “greybody factor”. 3 -0.2 -0.4 -0.6 D. Quantization and the Hawking effect -0.8 The field is quantized as a collection of self-adjoint operators ϕ̂m (t) satisfying the equation of motion (2.16) and canonical commutation relations. For each complex solution to the equations of motion we define the operator a(f ) = (f, ϕ̂) using the inner product (2.17). The commutation relations are equivalent to the relations -1 -1.2 -1.4 FIG. √ 2. v(x) for every lattice point at t = 27, for the case vmax = 2 and δ = 0.004. [a(f ), a† (g)] = (f, g) (2.18) for all f and g. For a positive norm solution p, a(p) acts like a lowering operator, while for a negative norm solution n, the conjugate n∗ has positive norm, so a(n) = −a† (n∗ ) acts like a raising operator. The Hilbert space of ingoing modes is just the Fock space built from solutions with positive frequency with respect to t (hence positive norm), and similarly for the outgoing Hilbert space. We assume the boundary condition that the ingoing positive frequency field modes are all in their ground state, C. Scalar field We consider a massless, minimally coupled scalar field on the lattice. In the continuum the action would be Z √ 1 Scont = d2 x −gg µν ∂µ ϕ∂ν ϕ 2   Z 1 (∂z ϕ)2 . (2.14) = dt dz (a(z, t)(∂t ϕ)2 − a(z, t) where we substituted the metric (2.8) in the second line. On the lattice we adopt the same action except that the partial derivative ∂z ϕ is replaced by the finite difference Dϕm := (ϕm+1 − ϕm )/δ and am (t) is replaced by its average over the relevant lattice sites:  Z X 1 2(Dϕm (t))2 2 S= dt am (t)(∂t ϕm (t)) − 2 am+1 (t) + am (t) m a(pin )|0i = 0. (2.19) The occupation number h0|a† (p)a(p)|0i of an outgoing normalized positive frequency wavepacket p in the in-vacuum |0i is is just minus the norm of the negative frequency part of the ingoing wavepacket that gives rise to p. To see this, suppose the ingoing wavepacket q+ +q− evolves to p according to the field equation (2.16), where q+ has positive frequency and q− has negative frequency. ∗ Then a(p) = a(q+ ) + a(q− ) = a(q+ ) − a† (q− ), hence (2.15) Varying the action (2.15) gives the equation of motion for ϕm (t),   2Dϕm−1 (t) = 0. (2.16) ∂t (am (t)∂t ϕm (t)) − D am+1 (t) + am (t) h0|a† (p)a(p)|0i = −(q− , q− ), (2.20) where we used the in-vacuum condition (2.19) and the commutation relations (2.18). If p is not normalized then to obtain the occupation number we must divide by (p, p). Hawking’s resultR[19] in standard field theory yields, for a wavepacket p = dω cω exp(−iω(k)t + ikz), the thermal occupation number Z ω|cω |2 NHawking = dω ω/T , (2.21) e H −1 There is a conserved (not positive definite) inner product between complex solutions: X ∗ ∗ (ψ, ϕ) = i am (t)(ψm (t)∂t ϕm (t) − ϕm (t)∂t ψm (t)). m (2.17) In two dimensions the continuum action (2.14) is conformally invariant, and all metrics are conformally flat, where TH = κ/2π is the Hawking temperature. 4 (in the free fall frame) changes sign and the wavevector is equivalent to +π/δ. From there it continues to decrease, winding up small and positive. This reversal of group velocity by a monotonic change of the wavevector on a lattice is analogous to the Bloch oscillation of an electron in a crystal acted on by a uniform electric field. Thus the interval (−π/δ, −|kc |) is mapped onto the interval (0, π/δ). The map is onto since every outgoing wavevector arises from some ingoing wavevector. This dynamical stretching of the k-space interval of length π/δ − |kc | to one of length π/δ is not problematic in a particle phase flow, but if we think about the wavepackets that are being approximated by this WKB analysis it is disturbing since there really are more outgoing modes than ingoing modes which produced them. It seems the wave evolution is somehow not conserving the number of wave degrees of freedom. The resolution of this puzzle is the essence of the Hawking effect. The wave evolution mixes positive and negative frequency modes, so a purely positive frequency outgoing wavepacket arises from a mixture of positive and negative frequency ingoing wavepackets. WKB evolution breaks down at the turing point, hence it misses this mode mixing. An example of a positive frequency WKB trajectory is shown in Fig. 4. The graph shows several different quantities as functions of t (in units with κ = 1): the static position coordinate of the wavepacket x and the horizon xH (solid lines), the free-fall frequency ω (3.1) (dashed line), the Killing frequency ωK (3.2) (dot-dashed line), and 50 sin(kδ/2) (sparsely-dashed line). The positions are scaled up by a factor of 10 so they can be seen on the same graph. In √ this example the lattice spacing is δ = 0.004, vmax = 2, and the trajectory has z(34) = 6 and k(34) = 1. III. WKB ANALYSIS The motion of wavepackets in the falling lattice model can be studied using the WKB approximation [12]. This amounts to using Hamilton’s equations for the position z and momentum k with a Hamiltonian determined by the dispersion relation H(z, k, t) = ω = ± 2/δ sin(kδ/2) a(z, t) (3.1) which is obtained by inserting a mode of the form exp(−iωt + ikz) into the field equation and keeping only the terms with the highest derivatives or finite differences of the field [20]. A plot of this dispersion relation is shown in Fig. 3 for a(z, t) = 1. Wavevectors differing by 2π/δ are equivalent, so only the range (−π/δ, π/δ) (the “Brillouin zone”) is plotted. 2 1 -3 -2 -1 1 2 3 -1 -2 FIG. 3. Dispersion relation ωδ = ±2 sin(kδ/2) plotted vs. kδ. The frequency ω is not conserved, even when kδ ≪ 1, simply because it is not the Killing frequency. Using (2.11) we see that in the continuum the Killing frequency would be given by 50 40 ωK v(x(z, t)) =ω+ k, a(z, t) (3.2) sin(kδ/2) 30 with which the dispersion relation can be expressed as a(z, t)ωK − v(x)k = (2/δ) sin(kδ/2). 20 free-fall frequency Killing frequency packet trajectory (3.3) 10 Note however that this relation is only valid when kδ is not too large, otherwise the fact that the Killing frequency involves a partial derivative in a direction that is not along a lattice point worldline renders the relation meaningless. It was shown in [12] that the WKB trajectory of an ingoing positive frequency wavepacket bounces off the horizon and comes back out provided the ingoing wavevector is greater in magnitude than some critical value |kc |. (This critical value depends on the time and place from which the packet is launched.) The wavevector evolves from negative values to positive values by decreasing until it is less than −π/δ, at which point the group velocity horizon 20 22.5 25 27.5 30 32.5 Time FIG. 4. WKB evolution. Notice first the behavior of sin(kδ/2). Backwards in time this starts out very small, rises to the maximum, and comes back down but not to zero. The group velocity in the frame falling with the lattice vanishes at the top, when kδ/2 = π/2, corresponding to a wavelength equal to two lattice units. This rise and fall happens via a monotonic change of k, and is the Bloch oscillation 5 mentioned above. The wavepacket that ends up with ω = ωK = 1 at late times starts out at early times with ω = ωK ∼ 43. The maximum possible lattice frequency is 2/δ = 500 in this case, so the ingoing frequency is certainly “high”. Whereas the free-fall frequency continues to redshift (forward in time) as the wavepacket climbs away from the horizon, the redshifting of the Killing frequency is essentially complete by the time the wavepacket reaches the turning point. Thus, past the turning point, each Killing frequency component propagates autonomously, more or less as it would in the continuum. As described in the Introduction, this, together with the fact that the evolution of ω at the lattice scale is adiabatic, is the reason why the usual Hawking effect is reproduced. B. Wavepackets and parameters If the Killing frequency were conserved as it is in the dispersive continuum models the entire Hawking spectrum could be probed just by propagating one wavepacket and decomposing it into is Fourier components, each of which would pick up a negative frequency part appropriate for the corresponding frequency. This is how Unruh checked the spectrum in his model [6]. Since the Killing frequency is not conserved on our lattice we cannot avail ourselves of this option.3 Thus instead we compute the occupation numbers for a sequence of wavepackets ϕ es (k) with different profiles: ϕ es (k) = k exp[−(k − 0.05s)2 ], s = 1, . . . , 9 (4.1) (in units with κ = 1, as usual). The form of these wavepackets was dictated by the need to have them well enough localized to be contained in the flat region of the spacetime (so we could construct the corresponding positive frequency initial data) while at the same time containing enough power at low frequencies of order κ to have a measurable Hawking occupation number. To shift the wavepacket to the desired starting position the Fourier transform (4.1) was multiplied by a phase factor exp(−ikzinitial). The positive frequency initial data was constructed by evolving each Fourier component one time step ∆t by multiplying with the phase factor exp(−iω(k)∆t) and Fourier transforming the result back into position space. The time dependence of the lattice (illustrated in Fig. 1) means that the results will not be independent of the “launch time” of the wavepacket from a given location. If we start our backwards evolution too far in the future then when the wavepacket nears the horizon the lattice points will be pathologically sparse. On the other hand if we launch backwards at too early a time then the wavepacket can become partially trapped in the region where the lattice spacing gets very small. This happens because a sufficiently high frequency wave is turned around when it tries to propagate into a region of larger lattice spacing. We encountered this “channeling” phenomenon in our simulations, and avoided it just by launching later. The runs reported here were all done with wavepackets that reached the horizon around t = 27, and were launched from far enough away to be contained in the flat region of the spacetime. Typically, IV. LATTICE WAVE EQUATION ANALYSIS In this section we discuss the results of our numerical calculations. The basic method is to take an outgoing positive frequency wavepacket and evolve it backwards in time. After it has “bounced” from the horizon and propagated back into the flat region of the spacetime we decompose it into positive and negative frequency parts. Since it is a purely left moving wavepacket at this stage the negative frequency part is just the part composed of positive wavevectors. The occupation number (2.20) of the original outgoing packet is then given by minus the norm of this negative frequency part divided by the total norm. We typically computed the total norm in position space using (2.17). Since the negative frequency part ϕ(−) was projected out in the wavevector space we usedPan alternate expression for its norm, (ϕ(−) , ϕ(−) ) = −2 k>0 |ω(k)||ϕ e(−) (k)|2 , where the frequency ω(k) is given by (3.1). A. Numerical issues Our system of equations (2.16) is a coupled set of ordinary differential equations, which is to be solved in the limit that time is continuous but the spatial lattice is fixed. In order to satisfy the Courant stability criterion the time step ∆t must satisfy ∆t < a(z, t)δ. (In this 1+1 dimensional setting the Courant condition coincides with the condition that causality be respected by the differencing scheme.) In the asymptotic region a(z, t) = 1 so the stability condition there is ∆t < δ. For a fixed ∆t/δ the general condition would be violated wherever a(z, t) < ∆t/δ. Since this only happens deep behind the horizon (for t > 0) we just modified a(z, t) inside the black hole so that it is never less than ∆t/δ. We used a time discretization scheme (written out in the Appendix) with errors of order (∆t)2 , and found that ∆t = 0.4δ was adequately small to obtain very good accuracy. 3 We studied the possibility that the WKB evolution could be used to establish a mapping between in and out frequencies, thus allowing us to check the entire spectrum of particle creation just by propagating a single wavepacket. This turns out to be a flawed idea however since the modes are mixed by the evolution and, moreover, such a map misses the role of the negative frequency piece. 6 the launch time was around t = 56 when the wavepacket was centered at x = 28. All of the calculations reported here were done using √ the line element (2.1), with vmax (2.4) ranging from 2 to 30. A sample of the comparisons for κδ = 0.002 is listed in Table I. The typical relative difference in this case is of order half a percent. 3. Lattice dependence C. Results To study the dependence on the lattice we ran the s = 1 wavepacket backward at several different values of the lattice spacing and for several values of vmax (2.4) √ from 2 to 30. The change of vmax produces only a rather small change in the shape of the function v(x) outside the horizon, however it produces a shift of the position of the horizon (2.5). For large vmax the horizon is given approximately by xH ≃ ln(2vmax ).√The horizon position ranges from about 0.6 for vmax = 2 to about 4 for vmax = 30. Thus a wavepacket launched backward in time from a given position at a given time will reach the horizon later in time for larger vmax , so that the proper lattice spacing at the horizon will be larger, though not by a very large amount for the wavepackets considered here. In Fig. 5 we plot the occupation numbers for lattice spacings δ ranging from 0.002√to 0.128, and for several values of vmax (2.4), vmax = 2, 2, 2.5, 3, 4, 5, 10, 20, 30. In all cases the late time wavepacket was launched from the same position and time. All wavepackets had the same envelope for their Fourier transform, however for different lattice spacings the wavepacket is necessarily different. 1. Generalities The behavior of a typical wavepacket throughout the process of bouncing off the horizon is illustrated in Fig. 6. The real part of the wavepacket is plotted vs. the static coordinate x at several different times. Following backwards in time, the wave starts to squeeze up against the horizon and then a trailing dip freezes and develops oscillations that grow until they balloon out, forming into a compact high frequency wavepacket that propagates neatly away from the horizon backwards in time. This ingoing wavepacket contains both positive and negative frequency components, in just the right combination to produce only an outgoing wave when sent in towards the black hole, with no wave propagating across the horizon. To illustrate this we have taken one of these ingoing wavepackets, decomposed it into its positive and negative frequency parts, and sent them separately back toward the horizon forward in time (Fig. 7). The negative frequency part (which is the smaller of the two since the surface gravity is fairly low compared to the typical frequencies in the outgoing wavepacket), mostly crosses the horizon, with just a bit bouncing back out. The positive frequency part mostly bounces, with just a little bit going into the black hole. Occupation # 0.025 2. Spectrum The “observed” occupation numbers for the wavepackets (4.1) all agree well with the integrated thermal predictions (2.21) of the Hawking effect provided the lattice spacing is not too large. Thermal 0.01563 0.01409 0.01266 0.01135 0.01014 0.00903 0.00801 0.00709 0.00625 Lattice 0.01557 0.01404 0.01262 0.01130 0.01009 0.00899 0.00797 0.00705 0.00622 Thermal Value 0.015 0.01 0.005 0 0 0.05 0.1 0.15 Delta TABLE I. Comparison of the black hole radiation on the lattice to the thermal Hawking occupation√numbers (2.21) for the wavepackets (4.1) for the case vmax = 2 and κδ = 0.002. s 1 2 3 4 5 6 7 8 9 0.02 FIG. 5. Occupation numbers for s = 1 wavepacket for various values of δ and vmax . Higher occupation numbers correspond to larger values of vmax . Rel. Diff. 0.0038 0.0035 0.0034 0.0041 0.0045 0.0039 0.0051 0.0052 0.0046 For all values of vmax , the occupation number converges to the Hawking prediction as the lattice spacing δ decreases. The pattern of deviations from the thermal prediction appears to depend quite strongly on the value of vmax , a result which we do not understand at present. The best we can do is to list various effects which would contribute to the deviations: (i) the WKB turning point is moving away from the horizon as δ grows, so one might ex- 7 pect that the effective surface gravity v ′ (xt.p. ) and hence the effective Hawking temperature felt by the wavepacket would decrease; (ii) at larger vmax the proper lattice spacing at the horizon is larger at the time this particular wavepacket reaches the horizon; (iii) the lattice points are rather sparse near the horizon for the larger values of δ, and their precise positions will be significantly different for different values of vmax . Effect (i) would lead to a decrease in the occupation number, while one might expect that effect (ii) would lead to a relative increase due to the extra particle creation associated with time dependence of the lattice. Perhaps what we are seeing is just a delicate balance between these two effects. in others, which requires at least two dimensions to be possible. At the atomic level such a volume-preserving flow involves erratic motions of individual atoms. One possible improvement of the lattice model is to make a lattice that mimics this sort of volume preserving flow. It is not clear whether the motions of the lattice points can be slow enough to be adiabatic on the timescale of the high frequency lattice modes. If they cannot, then the time dependence of this erratic lattice background will excite the quantum vacuum. In a fluid the lattice is a part of the system, not just a fixed background. The average flow could be adiabatic for the fully coupled ground state of the system but not for the field theory of the perturbations on the “background lattice”. Similar comments apply in quantum gravity: surely if in-out mode conversion is at play the incoming high frequency modes are strongly coupled to the quantum gravitational vacuum. Ideally, therefore, we should try to find a model in which the background is not decoupled from the perturbations. As a first step in this direction one could study a one dimensional model in which the lattice points are nonrelativistic point masses, coupled to each other by nearest neighbor interactions, and “falling” or propagating in a background potential (with or without periodic boundary conditions). The perturbations of such a lattice are the phonon field, and the back reaction to the Hawking radiation is included (although the background potential is fixed). In a model like this one could presumably follow in detail the nonlinear origin of the outgoing modes and the transfer of energy from the mean flow to the thermal radiation. V. DISCUSSION We have shown that linear field theory on a falling lattice reproduces the continuum Hawking effect with high fidelity, thus verifying earlier expections [12]. For the smallest lattice spacing we studied, δ = 0.002/κ, the deviations from the thermal particle creation in wavepackets composed of frequencies ω < ∼ O(κ) were of order half a percent. This is remarkable, particularly since the proper lattice spacing encountered by the wavepacket at the horizon was ∼ 0.08/κ. When κδ is pushed to higher values eventually we see significant deviations from the Hawking prediction. There are several directions in which this work could be developed. One idea is to include self-interactions for the quantum field theory. There are general arguments [21,22] that the thermal Hawking effect holds also for interacting fields. There are also perturbative calculations [23] and calculations exploiting exact conformal invariance of certain 1+1 dimensional field theories [24]. Also the “hadronization” of Hawking radiation from primordial black holes has been analyzed by matching to high energy collision data [25], but the physics of how the thermal Hawking radiation is “dressed” by interactions as it redshifts away from the horizon has not been studied. In a 1+1 dimensional model it should not be prohibitively difficult to numerically study Hawking radiation of interacting fields, now that we see it can be done on a lattice. The falling lattice model has provided a satisfactory mechanism—the Bloch oscillation—for how to get an outgoing mode from an ingoing mode in a stationary background, but there is a serious flaw in the picture: the lattice is constantly expanding. In the fluid model by contrast the lattice of atoms maintains a uniform average density. In a fundamental theory we might also expect the scale of graininess of spacetime to remain fixed at, say, the Planck scale or the string scale (since presumably the graininess would define this scale) rather than expand. Can the falling lattice model be improved to share this feature? A fluid maintains uniform density in an inhomogeneous flow by compressing in some directions and expanding ACKNOWLEDGEMENTS We are grateful to Matt Choptuik and David Garfinkle for advice on numerical matters and to Steve Corley for many helpful discussions. This work was supported in part by the National Science Foundation under grants No. PHY98-00967 at the University of Maryland and PHY94-07194 at the Institute for Theoretical Physics. 8 APPENDIX: FINITE DIFFERENCING SCHEME The equation of motion (2.16) for ϕm (t) is    ∂t am (t)∂t ϕm (t) = D 2Dϕm−1 (t) am+1 (t) + am (t)  . We used a time discretization scheme for the left hand side arising from the replacement of f˙ by h i f (t + ∆t/2) − f (t − ∆t/2) /∆t = f˙ + (1/24)f (3)(∆t)2 + O((∆t)4 ), (A1) (A2) which yields h    i ∂t (am (t)∂t ϕm (t)) → a(n + 1/2) ϕ(n + 1) − ϕ(n) − a(n − 1/2) ϕ(n) − ϕ(n − 1) /(∆t)2 , (A3) where the spatial index is omitted and the argument n stands for n∆t. We could have used this replacement, since a(z, t) isa prescribed function that can be calculated at any t, but for some reason we substituted a(n ± 1/2) by the  average a(n) + a(n ± 1) /2. This yields the approximation h     i a(n) + a(n + 1) ϕ(n + 1) − ϕ(n) − a(n) + a(n − 1) ϕ(n) − ϕ(n − 1) /2(∆t)2  1 1 1 1 aϕ(4) + ȧϕ(3) + äϕ̈ + a(3) ϕ̇ (n)(∆t)2 + O((∆t)4 ). = ∂t (a∂t ϕ)(n) + 12 6 4 4 (A4) (A5) With (A4) in place of the left hand side of (A1) we solve for ϕ(n − 1) in terms of ϕ(n) and ϕ(n + 1) to obtain an equation for evolving ϕ backwards in time by one time step given the following two steps. [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] T. Jacobson, Phys. Rev. D44, 1731 (1991). For a review see A.W. Peet, Class. Quant. Grav. 15, 3291 (1998). W.G. Unruh, Phys. Rev. Lett. 46, 1351 (1981). V. Moncrief, Astrophys. J, 235, 1038 (1980). M. Visser, Class. Quantum Grav. 15, 1767 (1998). W.G. Unruh, Phys. Rev. D 51, 2827 (1995). R. Brout, S. Massar, R. Parentani and Ph. Spindel, Phys. Rev. D 52, 4559 (1995). S. Corley and T. Jacobson, Phys. Rev. D 54,1568 (1996). S. Corley, Phys. Rev. D 57, 6280 (1998). Y. Himemoto and T. Tanaka, “A generalization of the model of Hawking radiation with modified high frequency dispersion relation”, gr-qc/9904076. H. Saida and M. Sakagami, “Black hole radiation with high frequency dispersion”, gr-qc/9905034. S. Corley and T. Jacobson, Phys. Rev. D 57, 6269 (1998). T.A. Jacobson and G.E. Volovik, Phys. Rev. D 58, 064021 (1998); JETP Lett. 68, 874 (1998). N.B. Kopnin and G.E. Volovik, Phys. Rev. B 57, 8526 (1998); Phys. Rev. Lett. 79, 1377 (1997). G.E. Volovik, Pis’ma Zh. Eksp. Teor. Fiz. 69, 662 (1999) [JETP Lett. 69, 705 (1999).] S. Corley, Phys. Rev. D 57, 6280 (1998). T. Jacobson, Phys. Rev. D 53 (1996) 7082. See e.g. Appendix A of ref. [12]. S.W. Hawking, Nature (London) 248 (1974) 30; Commun. Math. Phys. 43 (1975) 199. The Hamiltonian technique was first applied to dispersive models of the Hawking effect in ref. [7]. G.W. Gibbons and M.J. Perry, Proc. Roy. Soc. Lond., Series A 358, 467-94 (1977). W.G. Unruh and N. Weiss, Phys. Rev. D 29, 1656 (1984). D.A. Leahy and W.G. Unruh, Phys. Rev. D 28, 694 (1983). N.D. Birrell and P.C.W. Davies, Phys. Rev. D 18, 4408 (1978). J.H. MacGibbon and B.R. Webber, Phys. Rev. D 41, 3052 (1990); J.H. MacGibbon, Phys. Rev. D 44, 376 (1991); F. Halzen, E. Zas, J.H. MacGibbon, and T.C. Weekes, Nature 353, 807 (1991). 9 √ FIG. 6. A typical wavepacket evolution. Here vmax = 2, κδ = 0.004, and the wavepacket (4.1) has s = 1. The ocsillations of the incoming wavepacket are too dense to resolve in the plots. 10 FIG. 7. The positive (top row) and negative (bottom row) frequency parts of the ancestor of an outgoing positive frequency wavepacket, sent back towards the black hole. Note the vertical scale is different in the two plots. The wavepackets crossing the horizon in the two cases are negatives of each other since they must cancel to produce only the outgoing wavepacket we started with. The parts that go through the horizon appear on the right due to periodic boundary conditions. The high frequency ocsillations are too dense to resolve in the plots. 11