Hawking radiation and classical tunneling: A ray phase space approach

Acoustic waves in fluids undergoing the transition from sub- to supersonic flow satisfy governing equations similar to those for light waves in the immediate vicinity of a black hole event horizon. This acoustic analogy has been used by Unruh and others as a conceptual model for “Hawking radiation.” Here, we use variational methods, originally introduced by Brizard for the study of linearized MHD, and ray phase space methods, to analyze linearized acoustics in the presence of background flows. The variational formulation endows the evolution equations with natural Hermitian and symplectic structures that prove useful for later analysis. We derive a 2 × 2 normal form governing the wave evolution in the vicinity of the “event horizon.” This shows that the acoustic model can be reduced locally (in ray phase space) to a standard (scalar) tunneling process weakly coupled to a unidirectional non-dispersive wave (the “incoming wave”). Given the normal form, the Hawking “thermal spectrum” can be derived by invoki...


I. INTRODUCTION
The tilting of light cones in the curved spacetimes of the General Theory of Relativity (GTOR) has an acoustic analog: sound waves propagating in a background with flow.Recall that an event horizon arises in spacetimes when there are regions that are casually separated, as determined by the behavior of light rays. 1 An "event horizon" appears in the acoustic model when the background flow speed transitions through the local sound speed.This acoustic analogy was first noted by Unruh, 2 and later examined by Jacobson. 3(See Ref. 5 for a more recent summary of work in this area.)There are other analogs for Hawking radiation that have been explored in the literature, for example, Weinfurtner et al. 4 have experimentally examined surface waves on water flowing in a channel of varying width.In addition to the GTOR literature on this topic (see the references in Visser's article), there has also been some interest in laboratory analogs of Hawking radiation in the Atomic Molecular and Optical (AMO) literature, where similar effects have been proposed for light wave propagation within certain types of Bose-Einstein condensates.This has led to the study of "artificial black holes." 5 The equations for acoustic wave propagation near the event horizon exhibit an avoided-crossing-type of behavior.(This avoided crossing behavior was noted, for example, by Jacobson 3 in the study of black holes using quantum field theory.)Because of the possibility of tunneling between rays, short wavelength fluctuations in the immediate vicinity of the event horizon can escape.While Jacobson cites some of the plasma literature on mode conversion, his method of approach does not pursue a ray phase space viewpoint, 6 which we think helps to clarify things by highlighting those aspects of the theory that are representation-independent, hence of most importance.
The goal of the present paper is to examine the acoustic analogy using ray phase space methods. 7These methods can be used to systematically derive the normal form for the evolution equations in the immediate vicinity of the event horizon.The normal form is the simplest possible form for the governing equations, as defined by having the smallest number of terms, arranged in the most symmetrical manner, and the normal form displays those combinations of parameters that are invariant under various transformations, hence it reveals those combinations of parameters that are most important physically.For example, it displays that combination of parameters that will appear in the S-matrix connecting the incoming and outgoing wave amplitudes.See, for example, Ref. 8.
The current problem is not of standard type.It requires an extension of the normal form analysis presented in Ref. 8  or in Appendix F of Ref. 7 in order to identify the uncoupled ray Hamiltonians.The normal form involves a (scalar) tunneling process weakly coupled to a non-dispersive uni-directional wave (the "incoming wave," also a scalar field).Given the normal form, the Hawking "thermal spectrum" can be derived by invoking standard tunneling theory, but only by ignoring the coupling to the incoming wave.
The methods presented here can be extended to higher dimensions (for simplicity, in this paper we only discuss the case of one spatial dimension), and to include more complicated physics.For example, the variational methodology used here was originally introduced by Brizard 9 to study linearized MHD waves in the presence of flows.But we postpone consideration of those matters because of the greater complexity of the problem if we include magnetic fields.
The outline of the paper is as follows: In Section II, we present a heuristic treatment of the problem, highlighting key points that might get overshadowed when we dive into the technical details in Sections III-VI.In Section III, we discuss one-dimensional acoustics using a cold fluid model.We then show how to derive the acoustic model (including a background flow) using a variational principle.This will allow us to use standard methods to derive conservation laws and will lead us to a Hamiltonian formulation of the linearized wave equation.We follow Brizard 9 for this part of the paper, and note that the symplectic inner product that plays an important role in our work is also invoked in Ref. 4, but they do not use a ray phase space approach.In Section IV, we examine eikonal solutions of the linearized wave equations and show that we recover the dispersion function (3) away from the event horizon.In Section V, we construct the local wave operator in the vicinity of the horizon using phase space methods to find the normal form.We then solve the local wave equation and compute the S-matrix.
Elsewhere, we will discuss the discretization of the phase space variational principle.This provides a means to derive numerical schemes that are symplectic and have good stability properties.

II. HEURISTIC TREATMENT
Before diving into mathematical details, it is useful to summarize the main result, which is really quite simple if we allow ourselves to ignore some technical matters we will discuss later on.The treatment here is non-relativistic, so transformations between frames are Galilean.
Start with the wave equation for acoustic waves in a uniform and stationary background.The dispersion function in that case is the familiar Here, x is the wave frequency, and k is the wavenumber.The two roots of D ¼ 0 are left-and right-moving nondispersive waves that propagate at the sound speed 6c s .
If there is flow with fluid velocity v 0 , Doppler effects must be included, and in the lab frame the dispersion function becomes Note that if v 0 ¼ 6c s , a zero-frequency root appears for D 0 ¼ 0: This corresponds to a standing (frozen) wave pattern.Figure 1 shows the curves satisfying D 0 ¼ k þ k À ¼ 0. These are plotted on the ðv 0 ; kÞ-plane, rather than ray phase space (x, k), which allows us to see the entire range of behaviors this system can exhibit without having to specify a flow profile.For Figure 1, the sound speed is assumed to be constant c 2 s ¼ 1, implying there are "event horizons" at v 0 ¼ 61.We might guess that the local dispersion function governing sound waves in a nonuniform medium with flow is of the form (we now drop the prime) The notation we use here (and throughout the paper) emphasizes that, because the background is assumed to be timestationary, we can treat single-frequency waves by Fourier analysis.In these expressions, we want to view Dðx; k; xÞ as a function of the variables (x, k) which form a conjugate pair on the ray phase space.The frequency x appears as a parameter.
In WKB theory, the dispersion function Dðx; k; xÞ is the ray Hamiltonian.In a two-dimensional ray phase space, the zero locus of the Hamiltonian (D ¼ 0) is also the set of rays.For example, in Figure 2, we show the rays near the event horizon v 0 ¼ c s for the special case where c s ¼ 1, and use a linear velocity profile, v 0 ðxÞ ¼ 1 þ x, hence the event horizon is at x ¼ 0.
The arrows on the rays in Figure 2 indicate the direction of the ray evolution given by Hamilton's equations 7 (subscripts here denote partial derivatives) The smooth central branch that passes near the origin is a right-moving wave associated with the root k þ ¼ 0. This is the "incoming" wave.A little algebra shows that the group velocity for this wave equals 2c s at the event horizon.
The other branch, associated with k À ¼ 0, forms an avoided crossing, typical of tunneling.For those branches of the dispersion surface, as x approaches zero jkj ! 1.The Hamilton equations for these rays show that disturbances propagate from large-to small-k, meaning that to find the disturbances that escape from the vicinity of the event horizon into the subsonic region (to the left in Figure 2) we must specify initial data on the rays at the event horizon.In the FIG. 1.A plot of the surface D 0 ¼ 0, where D 0 is defined in (2), using x ¼ 1 and c s ¼ 1.The plot is shown in terms of k and v 0 .To guide the eye, we have also plotted the lines v 0 ¼ 61 and k ¼ 0 (all three of these lines are dashed).The avoided crossings, where case of black holes, "Hawking radiation" corresponds to emission of massive particles and photons from the vicinity of the event horizon due to quantum field fluctuations, presumably at very high jkj.
The notations ½0; 1; s 2 ; ð1 þ s 2 Þ in Figure 2 are the ratios for the energy on rays undergoing tunneling from a disturbance starting at large positive k, just to the right of the event horizon.Part of the energy for this disturbance tunnels through and "escapes" into the subsonic region.These tunneling connection coefficients are computed in isolation from the incoming wave, as explained in the text near Eq.(11), and derived in Section V.
Notice that to the left of the event horizon at x ¼ 0 the rays propagate in opposite directions (indicating the wave equation is bi-directional in that region), while to the right of the event horizon both rays propagate to the right.
We should also note that jkj ! 1 is usually interpreted as the signature of a "resonance."However, there is nothing fundamental in a mathematical sense about this characterization, since jkj going to infinity at a particular point in x can be removed by a linear canonical transformation.For example, we can use the linear symplectic transformation given in (108), which we will discuss in Section V. A more invariant way to characterize a resonance using normal form theory is given in Chap.6 of Ref. 7, and we will use that method here.
As will be shown in Section IV, the WKB-based intuition reflected in the expression (3) is valid away from the event horizon, but it breaks down in the immediate vicinity of the horizon because WKB is invalid there.A local treatment must be developed in order to compute the S-matrix connecting incoming and outgoing WKB solutions.We summarize the approach here briefly, with details given in Sections III-VI of the paper.
The linearized equations of one-dimensional gas dynamics, Eqs. ( 23)-(25), will be derived from a variational principle (30), and a Hamiltonian field theory is then developed using standard methods.The canonical field variables are the particle displacement, nðx; tÞ, and the momentum density, pðx; tÞ.Written in terms of this canonical formulation, time evolution is governed by a 2 Â 2 Schr€ odingertype equation (44 To focus on the key points of this heuristic discussion here, we simplify by setting the density q 0 equal to unity, ignoring Moyal terms in the product kv 0 ðxÞ (defined in (49) below), and using a constant sound speed c s .Notice that the symbol matrix A is not self-adjoint, but it becomes self-adjoint (for real x, k, and x) when we multiply it by i times the symplectic matrix, J, defined in (40).This important fact is a general characteristic of Brizard's theory. 9o isolate the tunneling process as much as possible, we need to bring the symbol matrix (5) into normal form.Then, the 2 Â 2 operator associated with that symbol will provide us with the simplest local wave equation, valid in the immediate vicinity of the event horizon.
To compute the normal form, we start with the eigenvalues and eigenvectors of the non-self-adjoint A. The eigenvalues are For a fixed (but arbitrary) frequency x, choose k 0 ðxÞ by solving k þ ðx ¼ 0; k 0 ðxÞ; xÞ ¼ 0 The next step in the normal form calculation is to Taylor expand all quantities about the point in phase space where the incoming ray intersects the event horizon, ½x ¼ 0; k ¼ k 0 ðxÞ.Writing k ¼ k 0 þ j, we keep only the leading terms in j in each of the expressions.(Higher order corrections can, of course, be included if they are deemed necessary.See Ref. 10 for a discussion.)After some lengthy algebra, summarized in detail later in the paper, we find we can recast the 2 Â 2 symbol governing the local interacting waves into the self-adjoint form We identify D þ and D À as the uncoupled dispersion functions at the event horizon, and a is the coupling between them.
(Note that the eigenvalues of Aðx; k; xÞ given in ( 6), k þ and k À , depend upon k, while D þ and D À depend upon j.) A plot of the surface D ¼ 0, now restricted to the region around the resonance at v 0 ¼ þc s ¼ 1, assuming a linear velocity profile of the form where we have set c s ð0Þ ¼ 1, and L ¼ 1.Note that this plot uses x as the horizontal axis.An avoided crossing type of behavior is present at the resonance, where k À ðx; kÞ ¼ 0, but there is an additional branch where k þ ðx; kÞ ¼ 0 is present as well, so this is not tunneling of the standard type.See text for details.
The coupling is usually evaluated at the base point of the Taylor expansion in ray phase space [here ðx ¼ 0 setting it equal to a constant value. 7This simplification is often sufficient to get good results.But in the present case, setting j ¼ 0 implies zero coupling.We carry the coupling term along to remind ourselves it is not exactly zero.Neglect of the coupling is required to recover the Hawking result outlined below, and in our derivation of the S-matrix in Section V. Testing the accuracy of the neglect of the coupling will be carried out numerically and reported in another paper.
The entry D þ is the dispersion function for a rightmoving wave with group velocity 2c s at the event horizon, as can easily be verified using the Hamilton equations with D þ as ray Hamiltonian.This is the "incoming wave." The entry D À is the wave undergoing tunneling when v 0 ðxÞ % c s .We now consider this avoided crossing separately, using where g 2 is defined as in Eq. ( 96).The ray equations, using D À as the ray Hamiltonian, are This implies that disturbances that start near x % 0 at very small spatial scales (large jjj) propagate toward smaller jjj, while moving away from the origin in x.We will continue to refer to Figure 2 in later discussions of the tunneling process, but in a slight abuse of notation the shift of the origin in k-space to k 0 ðxÞ given by ( 7) should be understood.
In Section V, we compute the S-matrix coefficients relating the general outgoing (complex) amplitudes on opposite sides of the horizon [see Eq. ( 103)] given a general incoming disturbance.For the special case shown in Figure 2, where there is no disturbance at large negative j, there is a jump in the complex amplitude as we cross the event horizon of the form a À ¼ Àisa þ , where a 6 are the amplitudes of the field to the left and right of the event horizon in the x-representation, and The energy associated with the outgoing transmitted ray is proportional to the absolute square of the wave amplitudes, hence ja À j 2 ¼ s 2 ja þ j 2 .In Figure 2, we show the ratios of energies assigned to the various rays for an initial disturbance that starts with large positive j, i.e., a fluctuation that starts just to the right of the event horizon, but which partly tunnels to the subsonic region (to the left in the figure).
Because the rays have jjj ! 1 at x ¼ 0, in the x-representation solutions of the local tunneling wave equation ( 98) have essential singularities.This leads to problems for numerical simulation.In Section V, will also show that the local wave equation ( 97) can be transformed into the standard tunneling form (114), a representation which should have better numerical properties.
To this point, these results are entirely classical.The connection with "Hawking radiation" from a black hole is through the following analogy.Rewrite the exponent in the energy transmission coefficient as [from (11), we find where h is Planck's constant, k B is Boltzmann's constant, and T eff is an "effective temperature" Note that this effective temperature is inversely proportional to the length scale.(In the famous result by Hawking, 12 the effective temperature of a black hole is inversely proportional to the mass, while the Schwarzschild radius is proportional to the mass.This means that the temperature and characteristic length scale at the event horizon are in the same inverse relation as here.)If we consider laboratory acoustics and transitional flow, for example, with jet nozzles, we can choose approximate values for the transition region at the throat of the nozzle.For example, following Unruh, we choose the length scale, L % 10 À3 m, and the sound speed, c s % 300 m=s.This gives sðxÞ % expðÀ10 À6 xÞ; (14) implying that we can only observe the "thermal" character of the emission if we look in the MHz range of frequencies, and this in an extremely turbulent environment.(Unruh 2 acknowledges that this is a challenging measurement to make, though he points out that it is easier than using a laboratory-scale black hole!)Our primary interest in this paper is in the theoretical formalism, in particular, an examination of the Unruh model from the perspective of ray phase space, with a view toward generalization to include more spatial dimensions and magnetic fields.The normal form method presented here requires a modification of earlier methods, and the symplectic inner product (80) plays a key role, which is new.The extension used here should be applicable to other problems where this symplectic structure appears (e.g., all those covered by Brizard's theory of linearized MHD).
This completes our heuristic summary of the main points.We now move to the technical details.

III. ONE-DIMENSIONAL HYDRODYNAMICS
The approach followed here is a special case of the general theory presented in Ref. 9. We use a cold ideal fluid model in one spatial dimension, x, which has as dependent variables the density, qðx; tÞ, the velocity, v(x, t), and the pressure p(x, t).These quantities obey the following evolution equations: and where the subscripts denote partial derivatives, and c is the ratio of specific heats.We recognize the first equation as the statement of mass conservation, the second of momentum conservation, and the third is required for the evolution following fluid trajectories to be adiabatic with the equation of state p / q c .

A. Acoustic waves
It is important to note that the derivation of the acoustic model given below is non-physical in the sense that we start with nonlinear one-dimensional inviscid fluid flow, then linearize around a given (time stationary) spatial profile in density and background flow.This means that we are ignoring the formation of shocks, which are well-known to occur in this system of partial differential equations (PDEs).Now linearize Eqs. ( 15), ( 16), and (17) to find the evolution equation for small amplitude (acoustic) waves in a nonuniform, stationary, background with flow.Write qðx; tÞ ¼ q 0 ðxÞ þ q 1 ðx; tÞ; v 0 ðx; tÞ ¼ v 0 ðxÞ þ v 1 ðx; tÞ, and pðx; tÞ ¼ p 0 ðxÞ þ p 1 ðx; tÞ.The constant is a formal expansion parameter.Expand in powers of and collect terms.The zeroth order equilibrium must satisfy We note that the last condition, combined with the first, implies p 0 / q c 0 .Notice, further, that these equilibrium conditions imply that the pressure, density, and velocity profiles are not independent.In particular, the sound speed profile is Eq. ( 65)] This is not independent of the background velocity profile, v 0 ðxÞ.A little algebra shows that This, of course, does not preclude these two velocities from becoming equal one another, which is the situation of interest to us.At first order in , we have (N.B.The first term on the RHS of Eq. ( 24) is missing in Eq. (15) of Ref. 9.) Following Brizard, we now replace the three fields ðq 1 ; v 1 ; p 1 Þ with the first order particle displacement n.The variations in the field quantities are determined by the particle motions through the following identities: These identities are inserted into the first-order evolution equations to derive the evolution equation for the particle displacement.In particular, the first-order momentum conservation law (24) becomes, after some straightforward but lengthy algebra Note that this result requires use of the zeroth order equilibrium conditions (18) and ( 19).Note also that (29) agrees with Eq. ( 22) of Brizard, 9 when that expression is reduced to one spatial dimension, and zero magnetic field.Now introduce the following variational principle (the overall factor of 1/2 will ensure that the canonical momentum is equal to the physical momentum density.) A standard calculation shows that the evolution equation ( 29) follows when we require stationarity with respect to the variation dnðx; tÞ, assuming the stationary background obeys (18).From this point, unless otherwise noted, time integrals are from t 0 to t 1 > t 0 , and all spatial integrals from À1 to þ1.
The variational principle (30) can be used to construct a Hamiltonian formulation using the following standard algorithm.First, define the Lagrangian density where Ĝn Fn þ @ @x q 0 v 2 0 @n @x ; ¼ c @ @x p 0 @n @x : (32) The canonical momentum density is which we identify as the physical momentum density.The Hamiltonian density is constructed using the Legendre transformation (first writing H pn t ðn; pÞ À Lðn; pÞ; (34) Notice that this leads to a Hamiltonian density that is quadratic in n (through its derivatives) and p, as expected.The Hamiltonian is the functional defined by integrating the Hamiltonian density over x and t, which we write in the form where we have introduced the inner product notation hvjki ð dx v Ã ðxÞkðxÞ: A little algebra verifies that The first of this pair of evolution equations is simply a rewriting of the relationship between n t and p, while the second is seen to be a recasting of (29), after using (32).The pair of canonical evolution equations ( 38) and (39) can be shown to have a Hermitian structure.This will prove valuable in deriving the normal form.First, introduce the 2 Â 2 symplectic matrix Second, define the two-component complex field w ðn; pÞ T .Third, the complex symplectic product of w a and w b is now introduced (The values of x and t are identical here, the symplectic product concerns the two-component vector indices.)Finally, define the (degenerate) canonical inner product on the linear symplectic complex vector space The evolution equations (38) and (39) are now written in the Schr€ odinger form i @w a @t Ba c w c ; where we sum over repeated indices.Explicitly Unless otherwise noted, derivatives act on all quantities to the right.
The Weyl symbol of an operator Â is defined as 7 The Weyl symbol mapping, which is invertible, takes operators and maps them to functions on ray phase space Â $ R aðx; kÞ: This mapping is linear, and topological, meaning that it preserves neighborhood relations in the two spaces.Since operators generally do not commute, this implies that the symbol of the operator product Â1 Â2 cannot be simply the product of the related symbols a 1 ðx; kÞ and a 2 ðx; kÞ.In fact, the symbol of the product is given by the Moyal product, denoted where The reader is referred to Chap.where In the current case, it is possible to show that B is selfadjoint with respect to the canonical inner product h j i can After integration by parts, a little algebra shows that implying that B † ¼ B with respect to h j i can as claimed.Before discussing the WKB analysis of the evolution equations, we note that using the concepts we now have in hand, we can introduce the phase space variational principle Explicitly This variational principle will prove useful when we derive the 2 Â 2 normal form.Also, it is the starting point to discretize the dynamics and derive symplectic integrators, which we will discuss in a separate paper.
It is sometimes useful to write the phase space variational principle in terms of real canonical fields ðn; pÞ, this variational principle becomes We have integrated by parts and introduced an overall constant factor of i=2 in order to cast this into a more standard form, using the fact that overall factors in the variational principle do not affect the resulting evolution equations.A short calculation leads back to Hamilton's equations ( 38) and (39), as required.We have introduced the sound speed (invoking the equation of state p 0 / q c 0 ) For there to be non-trivial solutions of (62), for a given x and x, at least one of the eigenvalues of A½x; k ¼ h x ðxÞ; x must vanish.These eigenvalues are denoted k 6 ½x; k ¼ h x ðxÞ; x x À ½v 0 ðxÞ 6 c s ðxÞkðxÞ: (66) These are a pair of Hamiltonians, one for each of the two types of rays.Rays of each type live on the surfaces k 6 ðx; k; xÞ ¼ 0.
Hamilton's equations are (using Eq. ( 4), now with k 6 as ray Hamiltonians): We can find the associated polarization for each ray by assuming the relevant dispersion relation is satisfied, using it in (62), and then solving for the associated null eigenvector.For example, write x ¼ ðv 0 6c s Þk.The null eigenvectors must satisfy The associated polarization vectors are where q 0 , c s are functions of x, and k is treated as a free parameter here.These have been normalized with respect to the symplectic inner product, as per the discussion after (80), which simplifies later expressions.eikonal solution, we set kðxÞ ¼ dh=dx, after we have solved for hðxÞ, as described below.The amplitude a(x) varies in a manner governed by the action conservation law, which we can derive as follows.Choose one of the two polarizations, for example, the "-" polarization.Using the polarization (70) with k(x) unspecified as yet, insert the ansatz (61) into the variational principle (59), keeping only the leading order terms, i.e., the derivatives act only hðxÞ, but not the background quantities, or a(x) and k(x), consistent with the leading order eikonal approximation.After some algebra, this gives A½a; h ¼ Stationarity with respect to variations in the amplitude implies which gives us the local dispersion relation Solving this for h x , we can now solve for hðxÞ Variation of (71) with respect to the phase gives us the action conservation law Therefore Clearly, the eikonal approximation breaks down near the event horizon where the denominator goes to zero.
At each x, we use the instantaneous polarization by setting kðxÞ ¼ h x ðxÞ in (70), as already mentioned.This, together with the results (74) and (76), gives the eikonal solution (61).
Therefore, away from the event horizon, the heuristic derivation of the local dispersion function discussed in the first part of Section II is a useful guide.But, the WKB approximation is invalid when v 0 ¼ 6c s , so how to proceed?We will follow the strategy outlined in Chap.6 of Ref. 7 by using the insight we have gained from the WKB analysis to construct the normal form, which isolates the tunneling behavior from the other "incoming" wave in the immediate vicinity of the event horizon to the greatest extent possible.

V. THE NORMAL FORM
To deal more carefully with the region near the event horizon, we return now to the phase space variational principle (58), which we reproduce here for ease of reference Recall the canonical inner product is defined in Eq. ( 43).The normal form is developed about a fixed point in ray phase space, x 0 and k 0 .Suppose we have an event horizon where v 0 ð0Þ ¼ c s ð0Þ.Then of course, we choose x 0 ¼ 0 as our base point.
Because the background is time-stationary, we can treat each frequency separately.Choose a fixed, but arbitrary, frequency x 6 ¼ 0, and introduce the constant (in x and k) polarization vectors and The polarization vectors are mutually orthogonal with respect to the complex form of the symplectic inner product, defined as We have normalized them so Xðe 6 ; e 6 Þ ¼ 61.
We expand about that point where the incoming ray crosses the event horizon at x ¼ 0. That is, we choose k 0 ðxÞ such that k þ ðx ¼ 0; k 0 ðxÞÞ ¼ 0 [see Eq. ( 6)], implying Using k 0 ðxÞ in the polarizations (78) and (79), construct a new ansatz for the wave functions, appropriate for the local region around the event horizon Note that we have the identity This is how the shift of the origin in ray phase space to x ¼ 0 and k ¼ k 0 ðxÞ appears in terms of wave functions.Note also that w x ðx; tÞ is a two-component object, while / 6 ðxÞ are scalars.The amplitudes / 6 ðxÞ are not assumed to be eikonal in form.We have made no approximations yet.The wave function (82) simply reflects a change of dependent variable by choosing a particular polarization basis for a singlefrequency solution.This new form for w x is inserted into (77).Because we have chosen a single frequency for the wave function, we drop the integration over t, retaining only the integral over x, which gives us the variational principle for the amplitudes Here and By construction, the operator-valued matrix D is selfadjoint with respect to the standard inner product (37).
[Recall that B is self-adjoint with respect to the canonical inner product (43)].Because the polarizations ê6 ðxÞ are constant in x, the operators D6 and ĝ are clearly just linear combinations of the entries of B. Note also that in moving from (77) to (84), no approximations have been made.The variational principle (84) is still general.We have simply chosen a particular polarization basis for the wave functions, one that best isolates the incoming wave from the tunneling process.
The calculation of the normal form for the local wave equation is carried out using Weyl symbol methods, as outlined in the discussion leading to (8).For general background densities and flow profiles, the entries of D are messy because of the derivatives acting on the background quantities q 0 ðxÞ; c s ðxÞ, and v 0 ðxÞ, in addition to the action on the amplitudes / 6 .This obscures what is going on, so let's simplify things and assume that qðxÞ ¼ q 0 ¼ 1, and take c s to be constant.(We need to retain the x dependence in v 0 ðxÞ to keep the resonance local.)We will also neglect the higher order terms in the Moyal series that appear in the symbol matrix (52), replacing v 0 ðxÞ Ã k and k Ã v 0 ðxÞ by v 0 ðxÞk.This means that we are assuming the background variation is on a long spatial scale compared to that of the amplitudes / 6 ðxÞ. 7We emphasize that this is a much less restrictive assumption than that used in WKB theory, because we do not assume any special form for the amplitudes / 6 ðxÞ.
The normal form transformation is summarized most compactly if we construct the constant (in x and k) matrix QðxÞ by using the polarizations (78) and (79) as the column entries QðxÞ ½ êþ ðxÞ; êÀ ðxÞ : (88) Note that QðxÞ is not unitary, but instead satisfies Note also that QðxÞ is not defined when x ¼ 0. The limit x !0 is singular, reflected in the fact that (89) remains true for all x 6 ¼ 0. The symbol for the 2 Â 2 wave operator associated with the 2 Â 2 Schr€ odinger form (44), derived using the variational principle (77), is The next step is to use (88) to compute the new representation of the symbol matrix, the one associated with (84), by noting that the variational principle is a bilinear form, hence the transformation induced by the change of polarization basis results in the congruence transformation 11 Dðx; k; xÞ iQ † ðxÞ Á J Á Aðx; k; xÞ Á QðxÞ: (91) This new matrix is self-adjoint by construction because iJ Á A is self-adjoint for all real x, k, and x.More details regarding the appearance of congruence transformations in ray phase space methods for multicomponent wave equations can be found in Appendix C 2 of Ref. 7. Some straightforward algebra leads to the result Here, j ¼ k À k 0 ðxÞ, and we have retained all terms linear in that quantity.The normal form isolates as much as possible the uncoupled dispersion functions, D 6 ðx; j; xÞ, which are closely related to ( 6), but with a shift in the origin in k.
Given a matrix symbol like Dðx; j; xÞ, the related operator is given by applying the usual correspondence j $ Ài@ x ; (93) entry by entry.For product terms like xj, the Weyl calculus ensures that we end up with the symmetrized product Now consider the tunneling process in isolation.Let's focus on the event horizon that occurs when v 0 ðxÞ % c s .In that case, the dispersion function D À ðx; j; xÞ is the one associated with the tunneling.Let's choose our origin such that v 0 ð0Þ ¼ c s , and linearize v 0 ðxÞ where L is the length scale characteristic of the flow at the event horizon.This implies The zero locus in x and j of D À ðx; j; xÞ ¼ 0, for x 6 ¼ 0, has the characteristic avoided-crossing (hyperbolic) shape for the rays.
To find the transmission coefficient of the tunneling process, we need to solve the wave equation for the amplitude function / À ðxÞ.That is, using (94) and (96), we need to solve This has the piecewise solution From Figure 2, we know that the initial data for the rays are fixed in j-space, so we Fourier transform The calculation is similar to the mode conversion case discussed in detail in Section 6.3 of Ref. 7. Some straightforward algebra leads to the following piecewise solution in the j-representation: Writing we can now identify the elements of the S-matrix, connecting the amplitude coefficients where sðxÞ ¼ e Àpg 2 ðxÞ : We choose to define (103) as the inverse of the S-matrix so as to retain the understanding that SðxÞ is the matrix connecting incoming to outgoing amplitudes.The incoming field amplitudes are b 6 , the outgoing amplitudes are a 6 , as seen by referring to Figure 2. Using Euler's Reflection Formula (Eq.(5.5.3) of Ref. 13) with z ¼ ig 2 À 1, it is straightforward to show that S À1 ðxÞ is unitary, which implies S is unitary as well, S † S ¼ 1, hence total energy is conserved The case shown in Figure 2 corresponds to setting b À ¼ 0. This implies as quoted in (11).Therefore, the energies on the two outgoing rays are related by ja À j 2 ¼ s 2 ja þ j 2 .From (106), the energy on the incoming tunneling ray, jb þ j 2 ¼ 1 þ s 2 , as indicated in Figure 2.
The solution (98) has an essential singularity at x ¼ 0, reflected in the ray behavior in Figure 2, where the rays connect to large j.A similar behavior can be expected locally for the full 2 Â 2 problem, now compounded by the presence of the incoming wave.This causes difficulties for numerical studies.It would be useful to find a representation of the problem that has better numerical properties.Here, we simply sketch the approach we are examining, results will be reported in a separate paper as this is still work in progress.
Staying with the tunneling problem in isolation, consider now the following linear canonical transformation: This is a rotation by 45 , and it puts (96) into the tunneling normal form (see, for example, page 243, Eq. (6.41), and Figure 6.7 on page 244 of Ref. 7) Linear canonical transformations on ray phase space induce related unitary transformations in the associated Hilbert space of wave functions.This change of representation is a generalization of the Fourier transform, called the metaplectic transform.For example, the unitary transform that takes wave functions in the x-representation to wave functions in the X-representation is where The inverse of (110) is simply See Appendix E of Ref. 7 for details.Because X and K form a canonical pair, if we choose the X-representation to write our wave equation, we have the familiar association K $ Ài @ @X ; (113) implying that the equation governing the mode shape for frequency x is The solution of this equation can be written in terms of parabolic cylinder functions, and the S-matrix elements connecting incoming and outgoing rays computed.The parabolic cylinder functions have much nicer behavior in the tunneling region than / À ðxÞ the original x-representation.
As already mentioned, the linear canonical transformation (108) from ðx; jÞ !ðX; KÞ generates a unitary transformation on the Hilbert space of wave functions.Therefore, so long as we keep track of incoming and outgoing pairings (using the rays), the S-matrix elements are unchanged.This allows us to compute the connection coefficients for the original x-representation.Preliminary results suggest that this approach is promising for numerical work, even when we include the incoming wave.

VI. SUMMARY
In this paper, we have sketched a methodology for studying linear acoustics in transitional regions using ray phase space methods.The main contribution of the paper is to show how to derive the normal form (92).This isolates the tunneling phenomenon from the non-resonant "incoming" wave, while also providing the leading order coupling between the wave undergoing tunneling and the incoming wave (while the coupling is weak, it is not zero, as assumed in the computation of the S-matrix given here).
Isolation of the tunneling and casting it into normal form uncovers the coupling constant for that process (g 2 x ð Þ, as defined in Eq. ( 96)), and a standard calculation to compute the S-matrix connecting the incoming and outgoing tunneling wave fields uncovers the Boltzmann factor in the transmission coefficient, which is at the heart of the theory of Hawking radiation.
In further work, we are pursuing the use of the phase space variational principle (77) as a tool for deriving symplectic integrators.Working in ray phase space, as opposed to only x-or k-space, implies we have a much wider class of transformations (the metaplectic transformations) available to simplify the problem, and to find the best representation for numerical work.We will report on this elsewhere.
2 of Ref. 7 for a complete discussion.Here, we simply quote results unproven.Using (49), the Weyl symbol of B is a jw b i can hw a j Ôw b i can ; 8w a ; w b : êðxÞ ðe n ; e p Þ T is assumed to vary on the same length scale as the background.Use this ansatz in (38) and (39).At leading order, assuming the derivative acts only on the phase, we get IV. WKB ANALYSISNow we introduce a single-frequency eikonal ansatz