Greybody factors for Schwarzschild black holes: Path-ordered exponentials and product integrals

In recent work concerning the sparsity of the Hawking flux [arXiv:1506.03975v2], we found it necessary to re-examine what is known regarding the greybody factors of black holes, with a view to extending and expanding on some old results from the 1970s. Focussing specifically on Schwarzschild black holes, we re-calculated and re-assessed the greybody factors using a path-ordered-exponential approach, a technique which has the virtue of providing a semi-explicit formula for the relevant Bogoliubov coefficients. These path-ordered-exponentials, (being based on a"transfer matrix"formalism), are closely related to so-called"product integrals", leading to quite straightforward and direct numerical evaluation, while avoiding any need for numerically solving differential equations. Furthermore, while considerable analytic information is already available regarding both the high-frequency and low-frequency asymptotics of these greybody factors, numerical approaches seem better adapted to finding suitable"global models"for these greybody factors in the intermediate frequency regime, where most of the Hawking flux is concentrated. Working in a more general context, these path-ordered-exponential techniques are also likely to be of interest for generic barrier-penetration problems.

In recent work [1] on the sparsity of the Hawking flux the present authors, (and two others), found it necessary to numerically evaluate certain greybody factors for the transmission of massless bosons through the Regge-Wheeler potential -the combined gravitational and angular momentum potential barrier surrounding a Schwarzschild black hole. While significant work on this topic dates back to the mid-1970's, see in particular references [2][3][4][5], we felt it useful to completely re-assess the situation in terms of a radically different formalism using path-ordered matrix exponentials [6][7][8], a formalism closely related both to "transfer matrices" and to the so-called "product integral " [9][10][11][12][13]. One of the virtues of using path-ordered matrix exponentials is that it gives a semi-explicit expression for the Bogoliubov coefficients, and so gives a somewhat deeper analytic understanding of the transmission and reflection probabilities [14][15][16]. Path-ordered matrix exponentials are also relatively simple to evaluate numerically, and one never actually has to solve a numerical differential equation, one just performs a numerical integral.
In counterpoint, there has also been a lot of work done on both low-frequency and high-frequency limits for the greybody factors. (See for instance references [2,17] for low-frequency limits, and references [18][19][20] for high-frequency limits.) Information at intermediate frequencies is harder to come by, and we shall use our numerical insights to try to develop some semi-analytic understanding of the intermediate frequency regime.
We emphasize that while in the current article we are interested in black hole physics, the path-ordered matrix exponential formalism [6][7][8] is a general-purpose tool, of wide applicability to both barrier penetration and scattering processes.

Strategy
Consider a Schrödinger-like scattering problem of the form For any such problem one can define Bogoliubov coefficients relating the asymptotic left and right free-particle states, and derive the exact path-ordered-exponential result [6]: A brief sketch of the derivation is given in appendix A; see also reference [6] for details. This result, and various generalizations and modifications thereof, has been used in a number of subsequent articles to develop general bounds on transmission probabilities, both generic and black hole specific. See for instance references [6][7][8], additional formal developments in references [14][15][16], and some applications in references [21][22][23]. The key point is that the Bogoliubov coefficients satisfy |α| 2 − |β| 2 = 1, and that the transmission probability is simply In the current article, instead of looking for bounds, we shall use this technique as a basis for numerically calculating greybody factors for the Schwarzschild black hole.
The key steps are to replace the position x by the tortoise coordinate r * , and to replace V (x) by the Regge-Wheeler potential [8,28]. Some formal developments are relegated to appendix B.

Path-ordered-exponentials and the product calculus
The definition of the path ordered exponential, for a real or complex valued m × m matrix function A(x) from some initial point x i to some final point x f , is simply: Here, (closely following the usual technical construction of the Riemann integral), we shall take x * k ∈ [x k , x k+1 ] to be a tag of some partition , with individual widths ∆x k = (x k+1 − x k ), and with the mesh (D = max {∆x k }) going to zero in the limit N → ∞.
This definition of the path ordered exponential is equivalent to (a matrix-valued, noncommutative) definition of the so-called "product integral" [9][10][11][12][13]. The product calculus is exactly the idea of defining a different notion of calculus based on infinitesimal products and divisions, (recall that ordinary calculus is based on infinitesimal sums and subtractions). In the language of the product calculus equation (3.1) becomes Here the product integral is defined as [9][10][11][12][13] x f , the tag x * k , the width ∆x k , and the mesh D, again being defined as before. The equivalence of the two definitions given in equations (3.1) and (3.3) can be deduced by noting that all the second order or higher times in exp(A(x * k )∆x k ) go to zero in the limit D → 0 [12]. (Note that if the matrix A(x) happens to be nilpotent, so that A(x) 2 = 0, (which is quite often the case), then we have the exact result that exp A(x * k )∆x k = I + A(x * k ) ∆x k , so the equality holds even before the limit of vanishing mesh is taken.) See for instance references [12,13] for an overview of the product calculus.
The following results for the product integral [9][10][11][12][13] should be especially noted: • If a < b < c, and both b a (I + A dx) and c b (I + A dx) exist, then Equation (3.5) is commonly known in the mathematics community as the Peano series, and in the physics community as the Dyson series. (There are also alternative expansions available in terms of the Magnus series and/or Fer series that might sometimes be useful. We shall not explore such options here.) This Dyson series is the physics community's standard method for approximating path ordered exponentials. Using the path-ordering operator P, equation (3.5) can be re-written as where now with σ (i) being any permutation such that where Θ(x) is the Heaviside step function.
When written in the form of equation (3.2) the transfer matrix of equation (2.2) is particularly simple to numerically evaluate. Defining we note A(x k ) 2 = 0. The transfer matrix can be compactly written as Here h = (x f − x i )/N , and we have chosen the right-tagged equipartition x k = x i + kh. This naive expression is extremely easy to compute numerically; however much like a simple Riemann sum, convergence can sometimes be rather slow. In order to improve convergence, Helton and Stuckwisch [10] introduce a polynomial approximation which has error bounded by H( is a bounded interval function determined by the matrix A, and p is the order of the polynomial approximation. (This is the product integral equivalent of a higher-order Simpson rule.) A minor technical issue is that Helton and Stuckwisch use a definition of the product integral based on right hand multiplication, that is, the order of the products in equation 3) by noting that in the limit N → ∞, Thus by sending A(x i ) → −A(x i ) and using the approximation in reference [10] we obtain the inverse of the product integral we set out to calculate. Since the Bogoliubov matrix satisfies we see that there is no loss of information; this method lends itself very well to finding the transmission probabilities.
The 5th-order approximation Helton and Stuckwisch introduce is this [10]: where A j = A[x k−1 + jh/4] for j = 0, 1, 2, 3, 4 and k = 1, 2, .., N . This somewhat unwieldy expression can more compactly be rewritten as [10]: (3.14) Here This scheme will be our primary tool for actually evaluating greybody factors. We first tested our technique in two situations, (the square-barrier potential, the deltafunction potential), where exact analytic results are known. We thereby verified that the numerical integration adequately approximated the exact results, (details mercifully suppressed), before turning to the Schwarzschild geometry and its associated Regge-Wheeler potential.

Particle emission from a Schwarzschild black hole
We now turn to the Regge-Wheeler potential, which is related to particle emission from black holes. In this situation there is no exact analytic expression for the transmission probability, neither by directly analysing incoming and outgoing waves, nor via the product calculus method. Instead we must evaluate the Bogoliubov coefficients numerically. (Oddly enough for this problem the exact wavefunctions are known in terms of Heun functions, though this observation is less useful than it might seem simply because not enough is understood regarding the asymptotic behaviour of these Heun functions. See references [24][25][26][27].) We shall work in "geometric units" where both G Newton → 1 and c → 1.

Setup
The probability of emission, (the greybody factors), of massless particles from a black hole can be found by analysing the appropriate wave equation in curved spacetime. For example, for a scalar field ψ(x) one considers where ∇ a is the covariant derivative with the Christoffel connexion associated with the spacetime metric g ab . The probability of emission is related to the ratio of the amplitude of inward and outward radial travelling solutions to equation (4.1) at spatial infinity.
It can be shown, using separation of variables, that for a non charged, spherically symmetric (Schwarzschild) black hole of mass M this reduces to finding the transmission probability for equation (2.1) with the Regge-Wheeler potential (see reference [28] and, for more recent background, reference [8] and references therein): Here, r is the usual radial coordinate and r * the so-called tortoise coordinate. They are explicitly related by where W (x) is the Lambert W function [8,[29][30][31]. In addition, is the principal angular momentum number, and s is the spin of the particle. s ∈ {0, 1, 2} for scalars, photons, and gravitons respectively, and ∈ {s, s + 1, ...}. There are 2 + 1 azimuthal modes for each principle angular momentum mode, , which in the case of spherical symmetry are equiprobable. The emission rate dN s (ω)/dt dω gives the total probability for an emission, per unit time, per unit frequency, of a particle of spin s and frequency ω. It is given by the sum over the transmission probabilities T ,s (ω), (equation (2.3)), for each principal and azimuthal angular momentum mode, multiplied by the probability for a particle to be in a given mode P ,s (ω). That is [2,[32][33][34]: where for a Schwarzschild black hole and integer spin particles the probability for the particle to be in a given mode is given by the Bose-Einstein distribution, and g is the number of polarizations for a given spin s. Equation (4.4) represents the rate of the emission of particles. Each particle carries one quantum of energy, ω, and so the energy emission is given by Another physically important quantity is the cross-section, σ(ω), which represents an effective area that embodies the likelihood of a particle to be scattered, (i.e. deflected), by the black hole. This is intimately related to the probability of transmission through the potential barrier [2] σ s (ω) = πω −2 ∞ =s (2 + 1)T ,s (ω). (4.7) In the high frequency limit M ω 1 this approaches the classical geometric optics cross-section σ ∞ = 27πM 2 [2]. This can now be used to define the dimensionless cross-section, S(x), by Now the Regge-Wheeler potential is asymptotically zero at both ends, (that is, we have V (r * ) → 0 as r * → ±∞). So using equation (2.3) the transmission probabilities can be calculated from the Bogoliubov coefficients, as given in equation (2.2). In this case the transfer matrix becomes where the potential, has now been re-written in terms of the dimensionless variables u * = r * /2M and u(u * ) = r(u * )/2M . In the product calculus formalism this leads to where It is also possible to do a little pre-processing by changing variables in the path-ordered integral. This might somewhat help analytic insight, but does not seem to improve the numerics. See appendix B.

Numerics
The calculation for the transmission probabilities for the Regge-Wheeler potential was numerically implemented in Python by using the polynomial approximation of equation (3.14). The integration region [−∞, +∞] was approximated by the finite range [−50, 350], which was found to introduce an error of at worst ∼ O(10 −9 ). Numerical convergence tests for the product integrals are summarized in figure 1.  Figure 2 shows the transmission probability for each of the scalar, photon, and graviton cases. It can be seen that the s = 0 case has transmission at the lowest frequencies, and that in each case larger values require higher frequencies before there is any transmission through the barrier. Furthermore for each eventually (at high enough frequencies) there is complete transmission through the barrier. This can be interpreted physically by observing the potential V (r * ) is lowest for the = 0 = s case and as such less energy (i.e. lower frequency) is required to pass through the barrier. For larger and s values the potential is higher and so more energy is required. Eventually any particle species in any mode will have enough energy to completely pass through the barrier, i.e. T ,s (x) → 1 as x → ∞.  These greybody factors (for the Regge-Wheeler potential in Schwarzschild spacetime) exhibit very strong similarities and relatively minor differences. The greybody factors almost seem to be translations of one another, and all appear to be similar to suitably shifted hyperbolic tangent functions. See figure 3. We shall discuss the implications of this observation more fully in the subsequent section on "modelling". Figure 4 shows plots of both the number and energy emissions rates, equations (4.4) and (4.6). It can be seen that for x 0.6 the emission rate rapidly becomes negligible, exponentially decaying to zero. The emission of particles is dominated by scalars, and the rate reduces with increasing spin. This can be understood from figure 2, in which it can be seen that for lower spin s the transmission probabilities become significant at smaller x; which coincides with the peak in the probability spectrum P ,s (x). For spin s most of the interesting physics is concentrated in the range x ∈ 0, (s + 1)/ √ 27 , whereas for total emissivity most of the interesting physics is concentrated in the range  Plots of the number (left) and energy (right) emission spectra for a Schwarzschild black hole, see equations (4.4) and (4.6). Note the logarithmic scale. The emission spectrum is dominated by scalar particles, and emission rates decrease with increasing spin. The total emission rates (i.e. summing over all particle species) are bounded by the sums over the geometric optics limit, σ = 27πM 2 , for each species (shown in red), and this limit is approached as M ω → ∞.
Finally figure 5 shows a numerical plot of the dimensionless cross-section, for each species of particle. As x increases the cross-section approaches the geometric optics limit, i.e. S(x) → 1. This happens more quickly for lower spins. Each species also ex-hibits oscillations which are due to the transmission probabilities becoming appreciable for increasing values, weighted by the 2 + 1 factor.

Modelling the greybody factors and cross-section
Given the numerical data, and in view of other information we have available, to what extent can we characterize and model salient features of the greybody factors? In the low-frequency limit we know that [2,17]: (4.14) At high frequencies we know that [18][19][20]: Unfortunately this conveys relatively little information beyond the fact that T (x) → 1 exponentially rapidly as x → ∞. At intermediate frequencies rather little qualitative or quantitative information regarding the greybody factors is available. In contrast, for the dimensionless (scalar) cross section S(x) more is known. At intermediate/high frequencies Sanchez gives the equivalent of [18][19][20]: However, inspection of figure 4 shows that the bulk of the total emission spectrum is concentrated in the range x ∈ [0, 1/ √ 27], where this is not all that good an approximation. In fact the Sanchez approximation predicts a negative cross section S s=0 (0) < 0 at x = 0. Can we do better?
An extremely simple toy model that captures most (but not all) of the relevant physics is to use a simple sigmoid function and take The parameter w ,s controls the (inverse) width of the transition zone from 0 to 1, and inspection of figures 2 and 3 indicates that the transition zone is close to 1/ √ 27 wide for all spins and angular momenta. This corresponds to w ,s ≈ 27.
The parameter x * ,s controls the location of the transition zone from 0 to 1. For any x only those modes with x * ,s < x contribute appreciable to the sum in the dimensionless cross section S(x). In fact a very crude approximation is This explains why the transition zones as displayed in figures 2 and 3 are roughly equally spaced, and explains the specific value of the spacing. In fact, inspection of figures 2 and 3 indicates that a tolerable approximation for all is Thus with these specific values for the parameters our simple toy model becomes In an attempt to further improve the model, one thing we looked at was the location and width of the transition zones. For the scalar case we found that gave a slightly better estimate for the location of the transition zones. (Presumably these are the first two terms of an asymptotic expansion.) We also found it advantageous to artificially adjust the width parameter w ,s → 33. Finally we tweaked the low-x behaviour by noting that [tanh(z 1 n )] n ≈ tanh z at low z, while still approaching unity at large z, and chose n = 3 as a good fit. With these modifications in place we now have .

(4.24)
This gives a remarkably good fit to the (scalar) greybody factors and cross section. See figures 8 and 9. While it must be admitted that the model appears complicated, there are actually only three free parameters, (the exponent 3 we have used in the tanh, the width parameter w ,s ≈ 33, and the shift in the location of the switchovers). The other parameters, (the C ,s , the asymptotic location of the switchovers, the presence of the number 27) are fixed by the known asymptotic behaviour of the greybody factors and cross section. Overall, this is a quite acceptable 3 parameter fit, both to the scalar greybody factors and to the scalar cross section.
Fits to s = 1 and s = 2 could be developed along similar lines, (by tweaking the exponent n in the tanh, the width parameter w ,s , and the shift in the location of the switchovers). In the interests of brevity we restrict attention to the scalar case.

Directly modelling the cross section
In counterpoint, we have also looked at improving the Sanchez approximation directly (without explicitly worrying about the underlying greybody factors). The best we have been able to come up with is this: . Décanini, Folacci, and Rafaelli [35], and Décanini, Esposito-Farese and Folacci [36], argue on semi-analytic grounds that the 32/27, (which Sanchez obtains on purely numerical grounds), should more properly be replaced by 8πe −π . Numerically and visually these two options are indistinguishable. Fits to s = 1 and s = 2 could be developed along similar lines, (by tweaking the width and oscillations of the subdominant term).
In the interests of brevity we restrict attention to the scalar case.

Discussion
So what have we learned?
• Path-ordered exponentials (transfer matrices, product integrals) are an effective way of first analytically expressing the Bogoliubov coefficients associated with a scattering problem, and second can then be turned into an efficient algorithm for numerically calculating the Bogoliubov coefficients when the underlying problem is not analytically solvable. This observation is generic, not black-hole specific.
• The path-ordered exponential formalism is the only way we know of to write down a more or less explicit formula for the Bogoliubov coefficients (and hence the transmission and reflection amplitudes) associated with a scattering problem.
• Turning specifically to the Regge-Wheeler equation for the Schwarzschild black hole, the product integral (specifically the 5 th order Helton-Stuckwisch algorithm, effectively a higher-order Simpson rule for product integrals), allowed us to quickly and efficiently calculate numerical greybody factors (which we then compared to older extant data from the 1970's and also used in our own recent work on the sparsity of the Hawking flux). Perhaps more interestingly, once one has enough easily manipulable data on hand, it becomes feasible to undertake some semianalytic model building to explore the structural details of the greybody factors.
• The 3 parameter fit we obtained to the (scalar) greybody factors seems quite good; likewise the 2 parameter fit we obtained to the (scalar) cross section seems quite good.
• These ideas could easily be extended to Reissner-Nordström, Kerr, and Kerr-Newman black holes, and via the Teukolsky master equation to higher spins. Generic "dirty" black holes (black holes surrounded by matter fields) can also be dealt with as soon as one derives a suitably generalized Regge-Wheeler equation.
We hope that these observations may be of wider interest, both in a general scattering context, and for black hole specific calculations.
Later on we will assume that V +∞ = V −∞ but in principle this is not required. In the two asymptotic regions there are two independent solutions [6] ψ ±i The ±i refers to right (left) moving modes, e iω ±∞ x (e −iω ±∞ x ), and ω ±∞ = √ E − V ±∞ . To analyse the transmission and reflection coefficients we consider the Jost solutions, J ± (x), which are exact solutions to equation (A.1) satisfying and Here α and β are the Bogoliubov coefficients, which are related to the reflection and transmission amplitudes by These Bogoliubov coefficients are for incoming/right moving waves which are partially scattered and transmitted by the potential V (x). The reflection and transmission probabilities are then given by That is, the probability for an incident particle to be reflected off or transmitted through the potential V (x) is given by R or T respectively. Note that by definition the sum of the probabilities for a particle to be reflected and transmitted must be unity: Now the second order Schrödinger equation (A.1) can be written as a Shabat-Zakharov system of coupled first order differential equations [7]. To do this write the wave function as where a(x), b(x) are arbitrary functions, "local Bogoliubov coefficients", and the auxiliary function, ϕ(x), is chosen such that it has a non zero derivative and To reduce the number of degrees of freedom we can impose the gauge condition, d dx ]. Then substitute equation (A.10) into equation (A.1). Using the gauge condition, equation (A.12), one obtains the following system of equations: .
This has the formal solution [6], in terms of a generalized position-dependent transfer matrix, dx .
(A. 15) Here "P exp" denotes a path ordered exponential operation. In the limit x i → −∞, x f → +∞ this becomes an exact expression for the Bogoliubov coefficients: dx . (A. 16) In the case V −∞ = V +∞ there is a natural choice for the auxiliary function, simply take ϕ(x) ≡ ωx, where for simplicity we have written ω ≡ ω ±∞ . With this choice equation (A.15) can be seen to reduce to Note the formalism is extremely general and flexible, the current application to greybody factors is just one example of what can be done. See for instance references [6][7][8], related formal developments in references [14][15][16], and various applications in references [21][22][23].

B Appendix: Formal developments for Regge-Wheeler
We are specifically interested in the problem Therefore, changing variables r * → r, we have: The benefit of these transformations is that the integral now runs over a finite range; the disadvantage is the rapid oscillations near the endpoints. From a theoretician's perspective this is now probably the most explicit formulation of the problem one can realistically hope for.