Hawking Radiation Near the Horizon

This repository contains code for calculating Hawking radiation measurements by an Unruh-DeWitt particle detector near the horizon of a Schwarzschild black hole.

Results and a discussion of the implementation can be found in the paper:

C.J. Shallue and S.M. Carroll, What Hawking Radiation Looks Like as You Fall into a Black Hole, 2024.

Installation

Required packages and their versions can be found in the hawkrad.yml Conda environment file. The easiest way to install the pre-requisites is to make a new Conda environment from that file:

conda env create -f hawkrad.yml

The code in this repository does not need any special installation because it is pure Python. The easiest way to use it is to add the top-level repository directory (i.e., the one containing the hawkrad directory) to the PYTHONPATH environment variable.

Data

As described in more detail below, the code makes use of asymptotic expansions of solutions to the radial wave equation. Tables of coefficients of these asymptotic expansions have been pre-computed and are archived at https://zenodo.org/records/14205556.

Usage

Measurements by an infalling detector

Measurements by an infalling detector can be generated by the script scripts/radial_response.py .

Required command-line flags:

• --Omega: Detector energy gap ($E$ in the paper) in units of $R_S^{-1}$. Can be a list.
• --omega: List of $\omega$ values to evaluate in units of $R_S^{-1}$.
• --ell: Which ells to evaluate.
• --r_mid or --tau_mid: $r_\text{mid}$ or $\tau_\text{mid}$ in units of $R_S$. Can be a list.
• --delta_tau: $\Delta \tau$ in units of $R_S$. Can be a list.
• --R_rest: Schwarzschild radial coordinate in units of $R_S$ at which the detector is instantaneously at rest ($R$ in the paper). Default is "inf" (meaning $R=\infty$).
• --asym_coeffs_dir: Directory containing tables of coefficients of the asymptotic forms of the radial modes (from this data archive or computed as described in more detail below).
• --output_dir: Output directory.

Additional options can be found by running radial_response.py --help. These control the precision of the integrals and ODE solvers, the convergence tests, etc., and can be used to optimize the runtime.

The main output files have the form terms_ell*.npy. Each of these corresponds to a particular value of $\ell$ and contains a numpy array of shape (n_omega, n_Omega, n_chi, 6), where n_omega is the number of $\omega$ values sampled, n_Omega is the number of detector energy gaps, n_chi is the number of switching functions each defined by a pair of values $(\tau_\text{mid}, \Delta\tau)$, and the remaining axis corresponds to the following quantities from the paper:

1. $F_\text{B}^\text{in}$
2. $F_\text{B}^\text{up}$
3. $F_\text{H}^\text{in} - F_\text{B}^\text{in}$
4. $F_\text{H}^\text{up} - F_\text{B}^\text{up}$
5. $F_\text{H}^{\overline{\text{in}}}$
6. $F_\text{H}^{\overline{\text{up}}}$

Each of these quantities needs to be integrated over $\omega$ (e.g. using the trapezoidal rule) and summed over $\ell$. Once integrated and summed, the detector response in the Boulware, Hartle-Hawking, and Unruh vacuum states can be calculated by:

$$\begin{align} \mathcal{F}_\text{B} &= F_\text{B}^\text{in} + F_\text{B}^\text{up}, \\\ \mathcal{F}_\text{H} &= F_\text{H}^\text{in} + F_\text{H}^{\overline{\text{in}}} + F_\text{H}^\text{up} + F_\text{H}^{\overline{\text{up}}}, \\\ \mathcal{F}_\text{U} &= F_\text{B}^\text{in} + F_\text{H}^\text{up} + F_\text{H}^{\overline{\text{up}}}. \end{align}$$

Measurements by a static detector

The instantaneous response rate of a static detector (i.e., with fixed Schwarzschild coordinates) near the horizon can be calculated using the script scripts/radial_static.py.

The required flags are

• --R: Schwarzschild radial coordinate in units of $R_S$. Can be a list.
• --ell_max: Maximum value of $\ell$.
• --omega: Same as in radial_response.py.
• --asym_coeffs_dir:Same as in radial_response.py.
• --output_dir: Output directory.

The output is saved in a file F_in_up.npy, which contains a numpy array of shape (n_ell, n_omega, n_R, 2), where n_ell is the number of $\ell$ values, n_omega is the number of $\omega$ values sampled per $\ell$, n_R is the number of $R$ values, and the remaining axis corresponds to the following quantities from the paper:

1. $|\Phi_{\omega\ell}^\text{in}(R)|^2$
2. $|\Phi_{\omega\ell}^\text{up}(R)|^2$

These two values can be combined to compute the instantaneous response of a static detector using the formulas in Hodgkinson et al. [2].

Implementation details

Units

The code uses units in which the fundamental constants $c, \hbar, k_B$ are 1. In addition, the Schwarzschild radius $R_S = 2 G M / c^2$, where $M$ is the mass of the black hole, is set to 1.

This means, for example, that lengths have units of $R_S$ and times have units of $R_S / c$.

In these units, the Hawking temperature of the black hole is

$$T_\mathrm{BH} = \frac{\hbar c^3}{8 \pi G M k_B} = \frac{\hbar c}{4 \pi R_S k_B} = \frac{1}{4\pi}.$$

Another choice of units, which is used in the paper, is to set $c = \hbar = k_B = 1$, and express quantities in units of the Hawking temperature. We can convert from the units used in the code to this convention by multiplying by the appropriate power of $4 \pi$. For example, a proper time interval $\Delta \tau$ in the code has units of $R_S$, which can be converted to units of $T_{BH}^{-1}$ by dividing by $4 \pi$.

Radial coordinates

Depending on the context, the Schwarzschild radial coordinate is encoded in 3 different ways in the code:

1. $r$, the Schwarzschild radial coordinate in units of $R_S$.
2. $x \equiv r - 1$. We use this representation because when $r$ is very close the horizon (within around 1e-16 for 64-bit floating point numbers), $r$ is effectively truncated to 1.0, whereas $x$ can represent values much closer to the horizon (within around 1e-308).
3. $r_*$, the Regge-Wheeler tortoise coordinate, defined by

$$r_* \equiv r + \log|r-1|.$$

Klein-Gordon modes

Notation

In the code (and in the remainder of this README), rhoIn and rhoUp correspond to the radial functions $\rho^\text{in}$ and $\rho^\text{up}$ in the paper, but re-scaled by a constant factor: rhoIn means $\rho^\text{in} / B^\text{in}$ and rhoUp means $\rho^\text{up} / B^\text{up}$. This is to facilitate specifying the boundary conditions numerically (see Hodgkinson [1] Section 6.1)

xiIn and xiUp in the code correspond to $\rho^\text{in}$ and $\rho^\text{up}$ with the normalization in the paper.

tildePhiIn and tildePhiUp in the code correspond to $\tilde{\Phi}^\text{in}$ and $\tilde{\Phi}^\text{up}$ in the paper.

Asymptotic solutions

We use power series expansions for the in modes asymptotically close to the horizon and for the up modes asymptotically far from the horizon, as described in Hodgkinson [1].

For the in modes, we write

$$\rho(r_*) = e^{-i \omega r_*} w(r)$$

and expand $w(r)$ in a power series

$$w(r) = \sum_{n=0}^\infty b_n x^{n},$$

for $x \ll 1$ (where $x = r - 1$ as defined above).

The $b_n$ coefficients can be calculated from an explicit recurrence relation [1].

For the up modes, we write

$$\rho(r_*) = e^{i \omega r_*} e^{v(r)}$$

and expand $v(r)$ in a power series

$$v(r) = \sum_{n=0}^\infty c_n r^{-n},$$

for $r \gg 1$.

The $c_n$ coefficients can be solved for algebraically in a recursive manner [1]. We calculated these expressions in Mathematica. The file c-coeffs-100.dat in the data archive contains the coefficients for $n \leq 100$ in a format that can be read by Mathematica's Get function. The file c-coeffs-100.csv contains the same coefficients in a format that can be read by the CoeffGenerator class.

The $b_n$ and $c_n$ coefficients can be calculated on the fly in Python using the CoeffGenerator class (the file c-coeffs-100.csv is needed for the $c_n$ coefficients). However, this is quite slow because we do it symbolically using sympy to avoid roundoff errors and keep the coefficients as precise as possible.

To avoid re-calculating the coefficients every time, we generated tables of $b_n$ and $c_n$ as functions of $\ell$ and $\omega$ for $n \leq 100$ (for the $b_n$'s) and $n \leq 50$ (for the $c_n$'s) accurate to $\sim 16$ decimal places (approximately machine precision) using Mathematica. These tables are found in the data archive discussed earlier. The parent directory of the in/ and up/ directories should be passed to the --asym_coeffs_dir flag. New tables can be generated using the script scripts/gen_asym_coeffs_tables.py.

Numerical ODEs

We use scipy to solve for the $\rho^\text{in}$ and $\rho^\text{up}$ functions using scipy.integrate.solve_ivp. We solve with $r_*$ as the independent variable. The system that we solve is

$$\frac{d}{dr_*} \begin{bmatrix} \rho \\\ d\rho / dr_* \\\ x \end{bmatrix} = \begin{bmatrix} d\rho / dr_*\\\ (V_\ell(r) - \omega^2) \rho \\\ 1 - r^{-1} \end{bmatrix},$$

where $x = r-1$, as defined above.

In order to solve for the in and up modes, we specify boundary conditions by evaluating the asymptotic solutions near the horizon (for the in modes) and far from the horizon (for the up modes).

Since solving these ODEs accounts for a significant component of the overall runtime, we made various optimizations to reduce the time spent solving them, which can be found in the file ode_system.py.

Citation

If you find this code useful, please consider citing the paper

C.J. Shallue and S.M. Carroll, What Hawking Radiation Looks Like as You Fall into a Black Hole, 2024.

References

[1] L. Hodgkinson, Particle Detectors in Curved Spacetime Quantum Field Theory (2013) https://arxiv.org/abs/1309.7281.

[2] L. Hodgkinson, J. Louko, and A. Ottewill, Static detectors and circular-geodesic detectors on the Schwarzschild black hole, Physical Review D (2014) https://doi.org/10.1103/PhysRevD.89.104002.