# Lensing of Vacuum Entanglement near Schwarzschild
Black Holes

João G.A Caribé, Robert H Jonsson, Marc Casals, Achim Kempf, Eduardo Martín-Martínez

► To cite this version:

João G.A Caribé, Robert H Jonsson, Marc Casals, Achim Kempf, Eduardo Martín-Martínez. Lensing of Vacuum Entanglement near Schwarzschild Black Holes. Physical Review D, 2023, 108 (2), pp.025016. 10.1103/PhysRevD.108.025016. hal-04038851

HAL Id: hal-04038851

https://hal.science/hal-04038851v1

Submitted on 1 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Lensing of vacuum entanglement near Schwarzschild black holes

João G. A. Caribé,¹ Robert H. Jonsson,²,³ Marc Casals,⁴,¹,⁵,⁶ Achim Kempf,⁷,⁸,⁹,¹⁰ and Eduardo Martín-Martínez⁷,⁸,⁹,¹⁰

¹Centro Brasileiro de Pesquisas Físicas (CBPF), Rio de Janeiro, CEP 22290-180, Brazil

²Max Planck Institute of Quantum Optics, Hans-Kopfermann-Str. 1, 85748 Garching, Germany

³Nordita, KTH Royal Institute of Technology and Stockholm University,
Hannes Alfvénsväg 12, SE-106 91 Stockholm, Sweden

⁴Institut für Theoretische Physik, Universität Leipzig, Brüderstraße 16, 04103 Leipzig, Germany

⁵School of Mathematics and Statistics, University College Dublin, Belfield, Dublin 4, Ireland

⁶Laboratoire Univers et Théories, Observatoire de Paris, CNRS, Université PSL,
Université de Paris, 92190 Meudon, France

⁷Institute for Quantum Computing, University of Waterloo, Waterloo, Ontario, N2L 3G1, Canada

⁸Department of Physics and Astronomy, University of Waterloo, Waterloo, Ontario, N2L 3G1 Canada

⁹Perimeter Institute for Theoretical Physics, 31 Caroline St N, Waterloo, Ontario, N2L 2Y5, Canada

¹⁰Department of Applied Mathematics, University of Waterloo, Waterloo, Ontario, N2L 3G1, Canada

<figure><img src="image_2.png" /></figure>

(Received 21 March 2023; accepted 27 April 2023; published 24 July 2023)

An important feature of Schwarzschild spacetime is the presence of orbiting null geodesics and caustics. Their presence implies strong gravitational lensing effects for matter and radiation, i.e., for excitations of quantum fields. Here, we raise the question whether the lensing manifests itself also in the vacuum of quantum fields, namely by lensing the distribution of vacuum entanglement. To explore this possibility, we use the method of entanglement harvesting, where initially unentangled localized quantum systems are temporarily coupled to the field at different locations. We find that for the Boulware, Hartle-Hawking and Unruh vacua in 3 + 1-dimensional Schwarzschild spacetime, the harvesting of vacuum entanglement is indeed greatly amplified near caustics. In particular, we establish that preexisting vacuum entanglement can also be harvested for lightlike separations.

DOI: 10.1103/PhysRevD.108.025016

## I. INTRODUCTION

The presence of entanglement between spatially separated degrees of freedom of a quantum field is a basic phenomenon that occurs even for free fields in the vacuum state on a flat background spacetime [1,2]. The origin of this vacuum entanglement can be traced back to the fact that, in wave equations, neighboring field oscillators must be coupled to each other in order to describe the propagation of waves. The coupling between the neighboring field oscillators is through spatial derivatives in the wave operators, such as the d'Alembertian and Dirac operator. It is this coupling between neighboring field oscillators that also causes the ground state of the local field oscillators to be an entangled state.

In curved spacetimes, curvature impacts the derivatives in the wave operators which then impact the entanglement

in the field. Therefore, curvature also impacts the field correlations. Conversely, it has been shown that the imprint that curvature leaves in the field correlators is actually complete in the sense that the metric can be reconstructed from the field correlators [3,4]. As was shown in Ref. [5], the metric can also be reconstructed from the correlations between local measurements of the field. The entanglement structure of quantum fields plays a fundamental role in investigations of phenomena from holography to Hawking radiation and the black hole information loss problem [6–16].

To probe the spacetime distribution of entanglement in a quantum field, a versatile method is to couple initially unentangled localized quantum systems to the field at different spacetime regions. The amount of entanglement that the localized systems acquired can be determined by standard methods [17,18].

When the entanglement acquired from the field by localized quantum systems is extracted from the entanglement that was preexisting in the field, this protocol has become known as entanglement harvesting. Entanglement harvesting has been investigated in a number of scenarios since first hinted at in Refs. [19,20], in both flat and curved spacetime. It has been proven that entanglement harvesting

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI. Open access publication funded by the Max Planck Society.
can capture the geometry [21] and topology [22] of the underlying spacetime. So far, the scenario where entanglement is harvested from the field in the presence of black holes has only been studied in very idealized scenarios such as 2 + 1-dimensional (Bañados-Teitelboim-Zanelli) black holes [23] and 1 + 1-dimensional spacetimes with horizons [24]. The question of entanglement harvesting near a black hole in 3 + 1 spacetime dimensions has remained open. However, this case is of particular interest since one can expect new phenomenology, for example, due to lensing (see Ref. [25]) and due to the fact that orbiting null geodesics and caustics can connect the two localized harvesting systems.

The phenomenology arising from orbiting null geodesics and caustics on communication through quantum fields close to a Schwarzschild black hole was addressed in Ref. [26]. Here, we investigate in detail entanglement harvesting in a four-dimensional Schwarzschild spacetime, for the cases of the Boulware, Hartle-Hawking and Unruh vacua, using tools similar to those applied in Ref. [26]. We take the localized quantum systems to be static localized two-level quantum systems with a monopole coupling to a Klein-Gordon field, i.e., so-called Unruh-DeWitt (UDW) detectors, or detectors for short. Within this setup, we analyze the impact of the presence of caustics and of the fact that the detectors can be connected by secondary null geodesics, including the case where detectors are placed close to the event horizon.

We find, that the presence of caustics alters the characteristics of entanglement harvesting in two particular ways in comparison to flat spacetime. First, due to a lensing-like effect caused by the focusing of null geodesics, the final entanglement between detectors can become greatly amplified near the caustics in comparison to comparably placed detectors away from caustics. Second, due to changes in the singularity structure of the Wightman function which happen when the field waves cross through caustics, we observe that timelike-separated detectors can harvest pre-existing entanglement from the field if they are aligned along secondary null geodesics, i.e. null geodesics which orbit half of the black hole and so they have passed through one caustic. To the best of our knowledge, such harvesting of preexisting entanglement for lightlike separations has never been observed before in any spacetime.

Section II B begins by discussing the treatment and properties of the Wightman function of a massless scalar field on Schwarzschild spacetime and, in particular, its global singularity structure. Section III introduces the detector model, its perturbative treatment and negativity

as entanglement measure for the detectors' final state. Section III B discusses when the entanglement between detectors can be attributed to the harvesting of preexisting entanglement from the field. Section IV presents our actual results for specific detectors and the article closes in Sec. V with a discussion and outlook. The appendices collect supplemental figures. Furthermore, Appendix A gives the calculations of the perturbative contributions to the detectors' state and Appendix B discusses the numerical techniques for the evaluation of the Wightman functions and integrals evolving the Wightman function.

We use the natural units in which $c = G = \hbar = 1$.

## II. WIGHTMAN FUNCTION IN SCHWARZSCHILD SPACETIME

In this section we introduce a quantum scalar field in Schwarzschild spacetime as well as the Wightman function when the scalar field is in three quantum states of interest (namely, Boulware, Unruh and Hartle-Hawking). The reader familiar with these may wish to skip to Sec. II B, reviewing literature results on the global singularity structure of the Wightman function as the field wave front passes through caustics of Schwarzschild spacetime which motivate our search for gravitational lensing of entanglement harvesting.

### A. Klein-Gordon quantum field in Schwarzschild spacetime

In this section, we briefly review the treatment of a massless Klein-Gordon field in Schwarzschild spacetime, and the expressions for the Wightman function of the three field states we consider, which are needed for the perturbative treatment of the field-detector interaction.

The line element of the outer region of Schwarzschild spacetime in Schwarzschild coordinates {$t \in \mathbb{R}$, $r \in (2M, \infty)$, $\theta \in [0, \pi]$, $\varphi \in [0, 2\pi)$}, is given by

$$ds^2 = -f(r)dt^2 + f(r)^{-1}dr^2 + r^2(d\theta^2 + \sin^2\theta d\varphi^2), \quad (1)$$

where $f(r) := 1 - 2M/r$ $M$ is the mass of the black hole and $r = 2M$ is the radius of the event horizon. In this outer region, we consider a scalar quantum field $\hat{\phi}$ obeying the Klein-Gordon (K-G) equation $\Box\hat{\phi} = 0$. Considering the Schwarzschild spacetime and availing of its spherical symmetry, a general real-valued solution for that equation can be written as

$$\hat{\phi}(\mathbf{x}) = \sum_{\ell=0}^{\infty} \sum_{m=-\ell}^{\ell} \int_{0}^{\infty} d\omega \left( \hat{\alpha}_{\ell m \omega}^{\text{in}} \phi_{\ell m \omega}^{\text{in}}(\mathbf{x}) + \hat{\alpha}_{\ell m \omega}^{\text{up}} \phi_{\ell m \omega}^{\text{up}}(\mathbf{x}) + \hat{\alpha}_{\ell m \omega}^{\text{in}\dagger} \phi_{\ell m \omega}^{\text{in}*}(\mathbf{x}) + \hat{\alpha}_{\ell m \omega}^{\text{up}\dagger} \phi_{\ell m \omega}^{\text{up}*}(\mathbf{x}) \right), \quad (2)$$
where $x$ denotes a spacetime point, $\hat{\alpha}_{\ell m\omega}^{\text{in/up}\dagger}$ are creation and annihilation operators and

$$ \phi_{\ell m \omega}^{\text{in/up}}(t, r, \theta, \varphi) = \frac{1}{\sqrt{4\pi\omega}} e^{-i\omega t} Y_{\ell m}(\theta, \varphi) \frac{R_{\ell\omega}^{\text{in/up}}(r)}{r}, \quad (3) $$

with $Y_{\ell m}$ being the spherical harmonic of degree $\ell$ and order $m$ and $R_{\ell\omega}^{\text{in/up}}$ are radial factors.

Substituting the field modes (3) into the K-G equation leads to the conclusion that $R_{\ell\omega}^{\text{in/up}}$ obeys

$$ \frac{d^2 R_{\ell\omega}^{\text{in/up}}}{dr_*^2} + (\omega^2 - V_{\ell}(r)) R_{\ell\omega}^{\text{in/up}} = 0, \quad (4) $$

where $r_* = r + 4M \ln|\frac{r}{2M} - 1| \in (-\infty, +\infty)$ is the tor-toise coordinate and $V_{\ell}(r) = f(r) \left(\frac{2M}{r^3} + \frac{\ell(\ell+1)}{r^2}\right)$ is an effective potential.

The two linearly independent solutions $R_{\ell\omega}^{\text{in}}$ and $R_{\ell\omega}^{\text{up}}$ are defined by the asymptotic boundary conditions

$$ R_{\ell\omega}^{\text{in}} \sim \begin{cases} e^{-i\omega r_*}, & r_* \to -\infty, \\ I_{\ell\omega} e^{-i\omega r_*} + \rho_{\ell\omega}^{\text{in}} e^{i\omega r_*}, & r_* \to \infty, \end{cases} \quad (5) $$

$$ R_{\ell\omega}^{\text{up}} \sim \begin{cases} I_{\ell\omega} e^{i\omega r_*} + \rho_{\ell\omega}^{\text{up}} e^{-i\omega r_*}, & r_* \to -\infty, \\ e^{i\omega r_*}, & r_* \to \infty. \end{cases} \quad (6) $$

Here, $\rho_{\ell\omega}^{\text{in/up}} \in \mathbb{C}$ are the reflection amplitudes and $I_{\ell\omega} \in \mathbb{C}$ is the incidence amplitude. For compatibility with Ref. [27] and ease of notation, let us also define

$$ \bar{R}_{\ell\omega}^{\text{in/up}} \equiv \frac{1}{rI_{\ell\omega}} R_{\ell\omega}^{\text{in/up}}. \quad (7) $$

Note that it suffices to calculate the modes for just $\omega \ge 0$ and use the following symmetries for $\omega < 0$:

$$ \bar{R}_{\ell, -\omega}^{\text{in}} = \bar{R}_{\ell, \omega}^{\text{in}}, \quad \bar{R}_{\ell, \omega}^{\text{up}} = \bar{R}_{\ell, -\omega}^{\text{up}}. \quad (8) $$

Calculating the response of the detectors to the interaction with the quantized field $\hat{\phi}$ requires, as we shall see in the next section, the Wightman function. The Wightman function, when the quantum field is in a state $|\Psi\rangle$, is defined as

$$ W^{\Psi}(x; x') \equiv \langle \Psi | \hat{\phi}(x) \hat{\phi}(x') | \Psi \rangle. \quad (9) $$

It is thus a two-point function satisfying the homogeneous K-G equation. Henceforth we shall only consider quantum states $|\Psi\rangle$ in regions of spacetime where they satisfy the Hadamard property [namely, Eq. (14)] is satisfied for the Wightman function in such states]. The Wightman function in Schwarzschild spacetime is given by [27]

$$ W^{\Psi}(x; x') = \frac{1}{(4\pi)^2} \sum_{\ell=0}^{\infty} (2\ell + 1) P_{\ell}(\cos \gamma) \\ \times \int_{-\infty}^{\infty} \frac{d\omega}{\omega} e^{-i\omega\Delta t} G_{\ell\omega}^{\Psi}(r, r'; \Delta t), \quad (10) $$

where $\gamma$ is the angular separation between the two spacetime points $x$ and $x'$ and $\Delta t = t - t'$. The integral kernel $G_{\ell\omega}^{\Psi}$ depends on the quantum state of the field. For the Boulware [28] ($\Psi = B$), Unruh [29] ($\Psi = U$) and Hartle-Hawking [30] states ($\Psi = H$) it takes the forms

$$ G_{\ell\omega}^{B} = \theta(\omega) (\bar{R}_{\ell\omega}^{\text{up}}(r) \bar{R}_{\ell\omega}^{\text{up}*}(r') + \bar{R}_{\ell\omega}^{\text{in}}(r) \bar{R}_{\ell\omega}^{\text{in}*}(r')), \quad (11) $$

$$ G_{\ell\omega}^{U} = \frac{\bar{R}_{\ell\omega}^{\text{up}}(r) \bar{R}_{\ell\omega}^{\text{up}*}(r')}{1 - e^{-2\pi\omega/\kappa}} + \theta(\omega) \bar{R}_{\ell\omega}^{\text{in}}(r) \bar{R}_{\ell\omega}^{\text{in}*}(r'), \quad (12) $$

$$ G_{\ell\omega}^{H} = \frac{\bar{R}_{\ell\omega}^{\text{up}}(r) \bar{R}_{\ell\omega}^{\text{up}*}(r') + \bar{R}_{\ell\omega}^{\text{in}*}(r) \bar{R}_{\ell\omega}^{\text{in}}(r')}{1 - e^{-2\pi\omega/\kappa}}, \quad (13) $$

where $\kappa = 1/(4M)$ is the surface gravity, and $\bar{R}_{\ell\omega}^{\text{in/up}}$ denotes the (rescaled) radial factor of the ingoing and upgoing solutions to the wave equation, as defined in Eq. (7). Appendix B discusses the numerical techniques used for the evaluation of the Wightman function.

The expressions above convey that in the Boulware state both ingoing and upgoing modes are in their ground state, whereas in the Unruh state the upgoing modes are thermalized, and in the Hawking state both ingoing and upgoing modes are thermalized.

## B. Singularity structure of the Wightman function

The Hadamard form for the Wightman function $W^{\Psi}(x; x')$ is an analytic expression which is defined in a local neighborhood<fnref n="1" /> of the source point $x$ and explicitly shows its singularity structure. Explicitly, it is (e.g., Refs. [31–33])

$$ W^{\Psi}(x; x') = \lim_{\epsilon \to 0^+} \frac{1}{4\pi^2} \left[ \frac{u}{\sigma + i\epsilon\Delta t} - v \ln(\sigma + i\epsilon\Delta t) + w \right], \quad (14) $$

where $u = u(x, x')$, $v = v(x, x')$ and $w = w(x, x')$ are regular and real-valued biscalars. The so-called Synge's world function $\sigma = \sigma(x, x')$ is equal to one-half of the square of the geodesic distance joining $x$ and $x'$, which implies that $\sigma$ is negative/zero/positive whenever that geodesic is, respectively, timelike/null/spacelike. The biscalars $u$ and $v$ are *uniquely* determined by the geometry of

<footnote n="1">¹More precisely, the Hadamard form is only defined in a normal neighborhood of $x$: a neighbourhood $\mathcal{N}(x)$ of $x$ such that every $x' \in \mathcal{N}(x)$ is connected to $x$ by a *unique* geodesic which lies in $\mathcal{N}(x)$.</footnote>
the spacetime whereas w in principle depends on the quantum state $|\Psi\rangle$. The term in Eq. (14) with u is called the direct part and the term with v the tail part.

In order to see separately the divergences of the real and imaginary parts of the Wightman function, the following distributional limits are useful:

$$ \lim_{\epsilon \to 0^+} \frac{1}{\sigma \pm i\epsilon} = \text{PV}\left(\frac{1}{\sigma}\right) \mp i\pi\delta(\sigma), \quad (15) $$

$$ \lim_{\epsilon \to 0^+} \ln (\sigma \pm i\epsilon) = \ln |\sigma| \pm i\pi\theta(-\sigma), \quad (16) $$

where PV denotes the principal value distribution. This readily yields the anticommutator

$$ \begin{aligned} \langle \Psi | \{ \hat{\phi}(\mathbf{x}), \hat{\phi}(\mathbf{x}') \} | \Psi \rangle &= 2\text{Re}(W^{\Psi}(\mathbf{x}; \mathbf{x}')) \\ &= \frac{1}{2\pi^2} \left[ u\text{PV}\left(\frac{1}{\sigma}\right) - v \ln |\sigma| + w \right] \end{aligned} \quad (17) $$

and the commutator

$$ \begin{aligned} \langle \Psi | [ \hat{\phi}(\mathbf{x}), \hat{\phi}(\mathbf{x}') ] | \Psi \rangle &= 2i\text{Im}(W^{\Psi}(\mathbf{x}; \mathbf{x}')) \\ &= \frac{-\text{sign}(\Delta t) i}{2\pi} [u\delta(\sigma) + v\theta(-\sigma)]. \end{aligned} \quad (18) $$

The anticommutator and commutator [which, when multiplied by $i\theta(\Delta t)$ yields the classical retarded Green function] are, respectively dependent and independent of the quantum state $|\Psi\rangle$ of the field.

Equation (14) shows explicitly the singularity of the Wightman function along $\sigma = 0$, i.e., when x and a x' in a local neighborhood of x are connected by a null geodesic. It is well known [34,35] that the Wightman function $W^{\Psi}(\mathbf{x}; \mathbf{x}')$ continues to diverge when x and x' are connected by a null geodesic globally, i.e., even when x' is not in a local neighborhood of x. The global singularity structure of the Wightman function in the case of Schwarzschild spacetime was unveiled in Ref. [36]: the divergence of $W^{\Psi}(\mathbf{x}; \mathbf{x}')$ follows a fourfold pattern, with the singularity type changing every time the null wave front passes through a caustic point (i.e., a spacetime point where neighboring null geodesics are focused; in Schwarzschild spacetime, because of the spherical symmetry, caustics lie along the line $\gamma = 0$ of the point of emission of the null geodesics as well as along the antipodal line $\gamma = \pi$). Specifically, the pattern for the leading divergence in the real part of the Wightman function is<fnref n="2" />

$$ \text{PV}\left(\frac{1}{\sigma}\right) \to -\delta(\sigma) \to -\text{PV}\left(\frac{1}{\sigma}\right) \to \delta(\sigma) \to \text{PV}\left(\frac{1}{\sigma}\right) \to \dots \quad (19) $$

and that in the imaginary part of the Wightman function is

$$ -\delta(\sigma) \to -\text{PV}\left(\frac{1}{\sigma}\right) \to \delta(\sigma) \to \text{PV}\left(\frac{1}{\sigma}\right) \to -\delta(\sigma) \to \dots \quad (20) $$

where we have omitted the coefficients of the singularity factors.<fnref n="3" /> As an example, the leading singularities of the real and imaginary parts of the Wightman function along the wave front before crossing any caustics are, respectively, $\text{PV}(1/\sigma)$ and $-\delta(\sigma)$ [corresponding to the direct part in Eq. (14)], whereas after the wave front has crossed one caustic point these turn into, respectively, $-\delta(\sigma)$ and $-\text{PV}(1/\sigma)$. This change in the singularity structure will have relevant consequences in entanglement harvesting as we show below.

The above is the leading singularity structure but there is a corresponding subleading singularity structure. Reference [37] showed that, for the imaginary part of the Wightman function, its subleading fourfold structure is

$$ -\theta(-\sigma) \to \ln |\sigma| \to \theta(-\sigma) \to -\ln |\sigma| \to -\theta(-\sigma) \to \dots \quad (21) $$

Reference [37] further conjectured that the subleading structure for the real part of the Wightman function is

$$ -\ln |\sigma| \to -\theta(-\sigma) \to \ln |\sigma| \to \theta(-\sigma) \to -\ln |\sigma| \to \dots \quad (22) $$

The first terms in Eqs. (21) and (22) of course correspond to, respectively, the imaginary and real parts of the tail term in Eq. (14).

### III. ENTANGLEMENT HARVESTING

#### A. Detector model and perturbative treatment

Alice and Bob will carry detectors that can locally measure the field around them. To model their detectors we will use the conventional Unruh-DeWitt particle detector model [29,38], which consists of a nonrelativistic quantum system coupled locally to a scalar quantum field.

<footnote n="2">²An exception to the validity of Eqs. (19) and (20) is at caustic points; also, by “$\sigma$” in these equations we mean a well-defined extension of the world function outside normal neighborhoods (see Ref. [37] for more about both of these points).</footnote>

<footnote n="3">³Reference [36] showed what the global singularity structure of the Feynman Green function $G_F(x', x) = i(\theta(\Delta t)W^{\Psi}(x; x') + \theta(-\Delta t)W^{\Psi}(x'; x))$ is, from which that of the Wightman function readily follows.</footnote>
The Unruh-DeWitt model is covariant and causal<fnref n="4" /> [39–42] and captures the main features of the light-matter interaction (e.g., atoms coupled to the electromagnetic field) when exchange of angular momentum between the field and the internal degrees of freedom of the detector is not relevant [43]. The covariant treatment and use of time-dependent perturbation theory to calculate the joint time evolution of the detectors and the field in arbitrary curved spacetimes can be found in Refs. [39,42]. In particular, in the context of entanglement harvesting this model has been extensively used throughout the literature and the particular perturbative approach employed here will use the same notation and conventions used in, among many others, Ref. [44]. Hence, here we only give a brief summary stating the most relevant expressions for the present work.

We model the particle detectors D = A, B as two-level systems with energy eigenstates $|g\rangle_D$ (ground state) and $|e\rangle_D$ (excited state) which are separated by an energy gap $\Omega_D$. The detectors couple to the field amplitude $\hat{\phi}(x_D)$ along their worldline $x_D(t)$ through the interaction Hamiltonian<fnref n="5" />

$$H_{\text{int,D}}^{t} = \lambda_{D} \eta_{D}(t) \frac{d\tau_{D}}{dt} \mu_{D}(t) \otimes \hat{\phi}(x_{D}(t)), \quad (23)$$

where $\lambda_D$ is a coupling constant which is dimensionless in (3 + 1)-dimensional spacetime, $0 \le \eta(t) \le 1$ is a real-valued switching function, $\tau_D$ is the detector's proper time and $\mu_D(t) = e^{i\Omega_D\tau_D(t)}|e\rangle\langle g|_D + e^{-i\Omega_D\tau_D(t)}|g\rangle\langle e|_D$ is the monopole operator. Note that the Hamiltonian of Eq. (23) generates time translation with respect to coordinate time $t$, but not the detector's proper time $\tau_D$. (For a detailed discussion of this point, see Refs. [42,45].)

In the scope of this work, we consider static detectors in Schwarzschild spacetime, i.e., detectors with constant spatial coordinates $x_D = (r_D, \theta_D, \phi_D)$. For such detectors

$$\frac{d\tau_D}{dt} = \sqrt{f(r_D)} = \sqrt{1 - 2M/r_D}, \quad (24)$$

and we choose the relation between coordinate time and detector proper time as $\tau_D(t) = \sqrt{f(r_D)t}$. As switching functions for the detector we use Gaussians which as a function of coordinate time read

$$\eta_D(t) = e^{-(t-t_{0D})/T_D^2}, \quad (25)$$

where $T_D$ denotes the switching width and $t_{0D}$ is the center of the switching function.

We assume the initial state of the system (at $t \to -\infty$) to be a product state between the ground states of the two detectors and a field state $\rho_\psi$:

$$\rho_0 = |g\rangle\langle g|_A \otimes |g\rangle\langle g|_B \otimes \rho_\psi. \quad (26)$$

Assuming that the field state has a vanishing one-point function, the expansion of the detectors' state after time evolution $t = 0...T$ in coordinate time is

$$\rho_{AB,T} = |g\rangle\langle g|_A \otimes |g\rangle\langle g|_B + \lambda_A^2 \rho_{A,T} + \lambda_B^2 \rho_{B,T} + \lambda_A \lambda_B \rho_{AB,T} + \mathcal{O}(\lambda^4). \quad (27)$$

Using the basis order $|g\rangle_A|g\rangle_B, |e\rangle_A|g\rangle_B, |g\rangle_A|e\rangle_B, |e\rangle_A|e\rangle_B$ the final state is represented by the density matrix

$$\rho_{AB,T} = \begin{pmatrix} 1 - L_{AA}^{\Psi} - L_{BB}^{\Psi} & 0 & 0 & (M^{\Psi})^{*} \\ 0 & L_{AA}^{\Psi} & L_{AB}^{\Psi} & 0 \\ 0 & L_{BA}^{\Psi} & L_{BB}^{\Psi} & 0 \\ M^{\Psi} & 0 & 0 & 0 \end{pmatrix} + \mathcal{O}(\lambda^4), \quad (28)$$

whose entries are

$$M^{\Psi} = -\lambda_{A}\lambda_{B}\int_{-\infty}^{\infty} dt \int_{-\infty}^{t} dt' \left( \eta_{A}(t) \frac{d\tau_{A}}{dt} e^{i\Omega_{A}\tau_{A}(t)} \frac{d\tau_{B}}{dt} \eta_{B}(t') e^{i\Omega_{B}\tau_{B}(t')} W^{\Psi}(x_{A}(t); x_{B}(t')) \right. \\ \left. + \eta_{A}(t') \frac{d\tau_{A}}{dt} e^{i\Omega_{A}\tau_{A}(t')} \frac{d\tau_{B}}{dt} \eta_{B}(t) e^{i\Omega_{B}\tau_{B}(t')} W^{\Psi}(x_{B}(t); x_{A}(t')) \right), \quad (29)$$

$$L_{DD'}^{\Psi} = \lambda_D \lambda_{D'} \int_{-\infty}^{\infty} dt \int_{-\infty}^{\infty} dt' \eta_D(t) \frac{d\tau_D}{dt} \eta_{D'}(t') \frac{d\tau_{D'}}{dt} e^{-i\Omega_D \tau_D(t) + i\Omega_{D'} \tau_{D'}(t')} W^{\Psi}(x_D(t); x_{D'}(t')), \quad (30)$$

where $M^\Psi$ takes complex values while $L_{DD'}^\Psi \ge 0$ takes non-negative values. For our case of static detectors and Gaussian switching functions, as shown in Eqs. (A8) and (A9), using the mode expansion of the Wightman function, these

<footnote n="4">⁴For pointlike detectors this is strictly true. For smeared detectors this is true in an approximated sense. For details see Refs. [39,40].</footnote>

<footnote n="5">⁵In the fully covariant formulation of the interaction a Hamiltonian density is prescribed [42]. For pointlike detectors, integrating this density over hypersurfaces of constant $t$ yields Eq. (23).</footnote>
expressions can be solved analytically up to one integration over the frequency of the field modes which needs to be performed numerically.

To assess and quantify the entanglement of the two detectors in the final state $\rho_{AB,T}$ we use its negativity. Its perturbative expansion is $\mathcal{N}^{\Psi} = \max[\mathcal{N}^{\Psi,(2)}, 0] + \mathcal{O}(\lambda^4)$ with

$$ \mathcal{N}^{\Psi,(2)} = \frac{1}{2} \left( \sqrt{(L_{AA}^{\Psi} - L_{BB}^{\Psi})^2 + 4|M^{\Psi}|^2} - L_{AA}^{\Psi} - L_{BB}^{\Psi} \right). \quad (31) $$

We see that whether the two detectors end up in an entangled state, is determined by a competition between the size of the correlating term $M^{\Psi}$ and the local noise terms $L_{AA}^{\Psi}$ and $L_{BB}^{\Psi}$. In particular, if the noise terms are equal, $L_{DD}^{\Psi} := L_{AA}^{\Psi} = L_{BB}^{\Psi}$, as will be the case in Sec. IV; then the negativity and, thus, the entanglement between the detectors vanishes if $L_{DD}^{\Psi} \ge |M^{\Psi}|$, i.e., when the noise overcomes the correlations. Notice that, because $L_{AA}^{\Psi}$ and $L_{BB}^{\Psi}$ are local noise terms for which the Wightman function is evaluated along a single detector's worldline, they contain no information about field correlations between the two regions where the two detectors are interacting with the field.

## B. When is the entanglement extracted versus generated?

Entanglement harvesting is an interesting process because it can demonstrate the presence of entanglement in a quantum field between different spacetime regions. For example, it is clear that when two initially uncorrelated detectors become entangled through their interaction with the field while remaining spacelike separated, the entanglement they acquire comes from “extracting” preexisting entanglement in the field (see, e.g., Refs. [20,44] in flat spacetime). However, when the detectors are not spacelike separated, contributions to the correlating term $M^{\Psi}$ in the leading-order perturbative correction to $\rho_{AB,T}$ arise which are independent of the quantum state of the field, as recently studied in Ref. [46]. In fact, these contributions can be calculated solely from classical data, consisting of switching functions, detector worldlines and the field classical Green function. Hence, they tell nothing about the quantum properties of the field.

To see this, we follow Ref. [46]. First, note that the imaginary part of the Wightman function, given by the commutator of the field operators, is independent of the state of the field. Only the real part, which is given by the anticommutator of the field operators, depends on the quantum state. We can use this to split $M^{\Psi}$ into two contributions as $M^{\Psi} = M_{+}^{\Psi} + iM_{-}^{\Psi}$, where $M_{+}^{\Psi}$ is obtained by replacing $W^{\Psi}$ by its real part (which is symmetric) on the right-hand side of Eq. (29), and $M_{-}^{\Psi}$ by replacing $W^{\Psi}$ with its

imaginary part (which is antisymmetric). Note that, in general, $M_{+}^{\Psi}$ and $M_{-}^{\Psi}$ are complex-valued.

Below we will encounter generic scenarios where the detectors become entangled while $M^{\Psi}$ is dominated by $M_{-}^{\Psi}$, and $M_{+}^{\Psi}$ is (almost or exactly) vanishing. In such a situation most<fnref n="6" /> of the entanglement between the detectors is not to be attributed to any preexisting entanglement in the field, as we shall next argue. In such a scenario the entanglement between the detectors, as measured by $\mathcal{N}^{\Psi,(2)}$, would remain unchanged if we replaced the initial field state by a state that resulted in the same values for $L_{AA}^{\Psi}$ and $L_{BB}^{\Psi}$ (the value of $M_{-}^{\Psi} \approx M^{\Psi}$, would also remain unchanged since it is state independent). In particular, we could replace the original state of the field by a state that has the same Wightman function within the regions where the detectors are coupled to the field, while containing no entanglement between those two regions.

Hence, in this new state, the entanglement between the detectors appears to be generated due to their sequential interaction with the field instead of being extracted from preexisting correlations in the field. This line of reasoning was used in Ref. [46] in order to argue that in these cases where $M_{+}^{\Psi}$ is (almost or exactly) vanishing, the entanglement acquired by the detectors should not be referred to as “entanglement harvesting” from the field.

In its turn, in a scenario with spacelike-separated detectors where all entanglement between the detectors is harvested from preexisting entanglement in the field, we have that $M^{\Psi} = M_{+}^{\Psi}$ because the commutator vanishes between the detectors. That is, in this clear-cut scenario, the contributions from the field's state-dependent anticommutator are the ones that transfer the entanglement from the field to the detectors.

In addition, another observation made in Ref. [46] indicates that for timelike-separated detectors the processes captured in $M_{+}^{\Psi}$ are also due to the harvesting of preexisting entanglement in the field as opposed to generation of entanglement through the interaction. This is based on the fact that a process that creates entanglement between the detectors by having them sequentially interact with the field also allows for communication from the first to the second detector: the field carries information between the two detectors and that can get them entangled independently of any preexisting entanglement in the field. However, it is well known that at leading order in perturbation theory, with which we are concerned here, such causal influence of one detector on the other is determined by the commutator of the field and independent of the anticommutator of the field [see, e.g., Eq. (24) of Ref. [41]]. Thus, as argued in Ref. [46], it appears plausible that, generally, the detectors are harvesting preexisting entanglement from the field if

<footnote n="6">⁶A quantitative study of the relative contributions of $M_{+}^{\Psi}$ and $M_{-}^{\Psi}$ to the entanglement acquired by the detectors for the case of flat spacetime can be found in Ref. [46].</footnote>
$M^\Psi \approx M_+^\Psi$ is dominated by contributions arising from the anticommutator.

#### IV. HARVESTING OF GRAVITATIONALLY LENSED VACUUM ENTANGLEMENT

This section contains the main results which demonstrate the impact of gravitational lensing, which at caustics refocuses null geodesics emanating from a common source, on entanglement harvesting. To this end, we consider static detectors placed close to the black hole horizon and compare their behavior with the well-studied case of static detectors in flat spacetime.

As one would expect, an important parameter that decides if and to what extent two detectors become entangled is their distance. The distance between (the static worldlines of) the two detectors can be defined in terms of various meaningful measures (see also Ref. [26]). In the present context (of static observers in Schwarzschild or Minkowski spacetime), the light propagation coordinate time is an intuitive choice of measure which we use henceforth when referring to the distance between detectors. This is the minimal amount of (Schwarzschild or Minkowski) coordinate time that it takes light to propagate from the spatial position of one detector to that of the other detector along a null geodesic.

Once the spatial positions of the two detectors are chosen and, thus, their distance is fixed, the interaction of the detectors with the field can still be made to happen at spacelike, lightlike or timelike separation by introducing a switching delay; this means that the switching function of the second detector is shifted with respect to the first one by a certain amount of coordinate time. For example, in the following we will assume that both detectors couple through the same switching function (25) and introduce the switching delay $\Delta_{BA} = t_{0B} - t_{0A}$ as the (coordinate time) difference between the centers of the switching functions. The switching delay can then be used to maximize the entanglement in the final state of the two detectors.

In Minkowski spacetime, as far as the impact of the switching delay is concerned, the entanglement in the final state between two detectors at a fixed distance is maximized when the switching functions are exactly null aligned [46], i.e., when the switching delay equals the light propagation time.<fnref n="7" /> As far as the impact of the distance between the detectors is concerned, its impact on the final entanglement is straightforward: if the distance between the detectors is increased (while allowing the switching delay to be optimal and leaving other coupling parameters unchanged), the entanglement between the detectors in their final state decreases and goes to zero beyond a certain

distance. This holds even if the switching delay is readjusted to uphold null alignment. Furthermore, in 3 + 1-dimensional Minkowski spacetime, the entanglement between exactly null aligned detectors is dominated by correlations generated by communication rather than by harvesting, because the $M^\Psi$ term is dominated by $M_-^\Psi$ rather than $M_+^\Psi$, as shown in Ref. [46].

This raises two questions regarding entanglement harvesting in Schwarzschild spacetime. First, does the presence of caustics enhance the ability of detectors to become entangled? What is more, can caustics make the harvested entanglement no longer decrease monotonically with distance?

Second, does the singularity structure of the Wightman function allow for the entanglement to be dominated by harvesting rather than signaling if the detectors' switching instead of being aligned along a direct null geodesic is aligned along a secondary null geodesic, i.e., a null geodesic that has passed through a caustic? As discussed in Sec. II B, here the singularity structures of the real and imaginary parts of the Wightman function are shifted, so that $M_+^\Psi$ could dominate over $M_-^\Psi$ because the real part now carries the $\delta(\sigma)$ singularity. This would then constitute entanglement harvesting between timelike-separated detectors, since the secondary null geodesics lie inside the causal cone which is bounded by the direct null geodesics between detectors.

The following results answer both questions in the affirmative.

Figure 1 shows the gravitational lensing effect on the harvesting of entanglement by two static detectors in Schwarzschild spacetime when the field is in the Boulware state. The scenario is as follows. Two identical detectors, with detector gap $\Omega_{AA} = \Omega_{BB} = 5 \, \text{M}^{-1}$ and $\lambda = \lambda_A = \lambda_B$,<fnref n="8" /> are placed on static worldlines at the same radial coordinate $r = 6.009M$. The angular separation $\gamma$ between the two detectors varies along the horizontal axes of each panel. Both detectors are coupled to the field through the Gaussian switching function (25), but the difference between the two centers of the switching functions is shifted by an amount of coordinate time $\Delta_{BA} = t_{0B} - t_{0A}$ which varies between the panels as indicated.

An inset in each panel visualizes how far null geodesics, emanating from detector A at time $t_{0A}$, propagate within the coordinate time interval $\Delta_{BA}$: the black filled circle represents the interior of the black hole horizon at $r = 2M$, a black ring around it represents all points with radial coordinate $r = 6.009M$ and the red dot represents the spatial location of detector A. The violet curve then indicates the spatial position of the null wave front after propagation time $\Delta_{BA}$. In particular, the different insets

<footnote n="7">⁷This statement is not entirely precise, as in Minkowski spacetime the maximum is located very close to—but not exactly at—null alignment of the switching functions; see Ref. [46].</footnote>

<footnote n="8">⁸Note that because $\lambda$ is dimensionless and all terms are of order $O(\lambda^2)$ the numeric value of $\lambda$ only affects the result by an overall factor. We choose $\lambda = 1$ for improved readability.</footnote>
<figure><img src="image_3.png" /><figcaption>FIG. 1. Gravitational lensing of entanglement harvesting from the Boulware state: the plots show the entries $|M^B|$, $|M_{\pm}^B|$ and $|L_{DD}^B|$ of the final detectors' state (28) for two static detectors placed at radial coordinate $r = 6.009M$ with varying angular separation $\gamma$. Blue shading indicates where $|M^B| > L_{DD}^B$, i.e., the negativity (31) is positive and the detectors become entangled. All detector parameters are equal ($\lambda_D = 1, \Omega_D = 5M^{-1}, T_D = 1M$); only the delay $\Delta_{BA} = t_{0B} - t_{0A}$ between the switching functions (25) varies between plots. The insets indicate how far a null wave front propagates from the red point within a coordinate time interval $\Delta_{BA}$. In particular, at $\Delta_{BA} \approx 20.7386M$ at the antipodal point of the red point, i.e., at $\gamma = \pi$, the wave front intersects itself and the first caustic point forms. (See Supplementary material for an animated version of this plot [47].)</figcaption></figure>
show the formation of caustics located at angular separation $\gamma = \pi$ from the point of origin (which is plotted in red) where the null wave front intersects itself. The caustic reaches $r = 6.009M$ at $\Delta_{BA} \approx 20.7386M$.

The four different curves in each panel show the resulting contributions to the final density matrix of the two detectors (28), when the field is in the Boulware state. (Accordingly they are denoted with a superscript $\Psi = B$.) The horizontal, red and dashed line shows the single detector noise term $L_{DD}^B := L_{AA}^B = L_{BB}^B$. Because it only depends on the radial coordinate and the switching function's width $T_D$, it is equal for both detectors and constant for all positions of detectors considered here. The blue solid curve shows the absolute value $|M^B|$ resulting from the chosen value of the switching delay $\Delta_{BA}$ in each panel, as a function of the angular separation between the detectors. The two remaining dashed curves show the absolute values $|M^B|$ in green and $|M^B|$ in orange which, as discussed in Sec. III B, are obtained by replacing the Wightman function by its imaginary and real parts, respectively. In particular, as discussed in Sec. III B, if $|M^B|$ is dominated by $|M^B|$ (obtained from the imaginary part), then the detectors' entanglement is predominantly generated in the sequential interaction of the detectors with the field. On the other

hand, if $|M^B|$ is dominated by $|M_+^B|$ then final entanglement between the detectors is predominantly harvested from preexisting entanglement in the vacuum of the field.

Altogether the panels in Fig. 1 illustrate the lensing effect which appears in proximity to the caustics. The first two panels show the settings with the shortest switching delay $\Delta_{BA}$ between the switching functions. In these two plots, we see that the correlations between the detectors, as captured by $|M^B|$, are largest when the centers of the intervals during which the detectors couple to the field are connected by (direct) null geodesics. In the first panel there is a certain interval of angular separations in which the blue line for $|M^B|$ exceeds the noise term $L_D^B$ and the detectors end up entangled. However, in the second panel the detectors' final state remains separable, even for null separated detectors, because $\Delta_{BA}$ is increased. Note that here the contribution from $M_-^B$ dominates the peaks of $M^B$; hence entanglement between the detectors cannot be attributed to harvesting of preexisting entanglement from the field but to entanglement generation due to sequential interaction.

In flat spacetime no entanglement would be observed in the final detectors' state for larger switching delay between the detectors. Intuitively this is related to the growth of the light cone's surface area which dilutes the entanglement.

<figure><img src="image_4.png" /><figcaption>FIG. 2. Same data as Fig. 1 from a different perspective: the switching delay $\Delta_{BA} = t_{0B} - t_{0A}$ is plotted on the horizontal axis, while the panels show detectors with four different angular separations. (See Fig. 6 for a logarithmic plot.)</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 3. Relative differences between correlation terms in the Boulware state, and the Hartle-Hawking and Unruh states, respectively, i.e., plotting $\frac{|M_H^H - M_B^B|}{|M_B^B|}$, $\frac{|M_U^U - M_B^B|}{|M_B^B|}$, $\frac{|M_H^H - M_+^B|}{|M_B^B|}$ and $\frac{|M_+^U - M_+^B|}{|M_B^B|}$, in analogy to Fig. 1.</figcaption></figure>
In Schwarzschild spacetime, however, where the light cone refocuses at caustics the opposite can happen. As seen in the third through sixth panels (in order of increasing $\Delta_{BA}$), the detectors can become entangled at larger switching delays $\Delta_{BA}$ again, at angular separations close to $\gamma = \pi$ in the proximity of caustics. This answers in the affirmative the first of the two guiding questions raised above.

The second of the two questions raised above pertains to the Wightman function's singularity structure described in Sec. II B. Its effect is easy to recognize when comparing the last (bottom right) panel to the first (top left) panel of Fig. 1. In both panels the correlations between the two detectors, as captured by $|M^B|$, are maximal for null aligned detectors. However, in the first panel the detectors are connected by primary null geodesics, whereas in the last panel they are connected by secondary null geodesics. At secondary null geodesics, as discussed in Sec. II B, the singularity structure of the Wightman function is shifted from primary null geodesics so that now the real (anticommutator) part carries the $\delta(\sigma)$ singularity, which used to be carried instead by the imaginary (commutator) part in the case of primary null geodesics. Accordingly, the two contributions $|M^B_\pm|$ have a qualitatively similar overall shape in the first and last panels, except that the $M^B_+$ and $M^B_-$ have swapped places. For detectors aligned along a secondary null geodesic, the real part contribution $M^B_+$ dominates over the imaginary part contribution $M^B_-$. So if these correlations overcame the noise, $|M^B| > |L_{DD}|$, the detectors would become entangled by harvested entanglement. However, the correlations in the last panel are too weak for the detectors to end up entangled. To find detectors that get entangled by genuinely harvested correlations around secondary null geodesics, we need to position the detectors closer to caustics so as to make use of the overall enhancement of the Wightman function there. Such a setting is seen in the fifth panel (for $\Delta_{BA} = 21.5M$) of Fig. 1. In this and the following panels of Fig. 1, the switching delay $\Delta_{BA}$ is larger than the shortest light propagation (Schwarzschild coordinate) time between detectors with angular separation $\gamma = \pi$, and the detectors (at $r = 6.009M$) which are null aligned here are aligned along a secondary null geodesic and timelike separated. With this we have answered in the affirmative the second of the above questions, as we observe the genuine harvesting of preexisting entanglement from the field by timelike-separated detectors, aligned along secondary null geodesics. In the following, we will see further examples allowing for entanglement harvesting between timelike-separated detectors with angular separation $\gamma = \pi$.

Complementary to Fig. 1, the same scenario and effects are seen in Fig. 2 from a slightly different perspective. Here, the contributions to the detectors' final density matrix are plotted over the coordinate time switching delay $\Delta_{BA}$, while the four different panels correspond to four different angular separations. (Figure 6 provides a logarithmic plot of the same data.) Note that these plots over $\Delta_{BA}$ are directly

comparable to the plots also found in Ref. [46] in $3+1$-dimensional Minkowski spacetime.

In the first panel of Fig. 2, for angular separation $\gamma = \pi/5$, there are three peaks appearing in $|M^B|$ which correspond to an alignment of the switching functions along primary, secondary and tertiary null geodesics, respectively. The qualitative structure of the first peak and its contributions from $|M^B_+|$ and $|M^B_-|$ correspond to the structure found in flat spacetime in Ref. [46]. However, in the secondary peak the qualitative roles of $|M^B_+|$ and $|M^B_-|$ are interchanged due to the shifted singularity structure of the Wightman function, as is easy to see in the logarithmic plot in Fig. 6. For the tertiary peak, together with the singularity structure of the Wightman function, the qualitative structure of the peak also shifts back to its primary form.

Even if there are three peaks appearing at angular separation $\gamma = \pi/5$, only the first one corresponding to alignment along the direct, primary null geodesic exhibits entanglement in the detectors' final state. In an intermediate regime of angular separations, for $0.64\pi \lesssim \gamma \lesssim 0.81\pi$, not even the primary peak overcomes the noise, and the detectors' final state remains separable for all switching delays $\Delta_{BA}$. An example of this is seen in the third (bottom left) panel of Fig. 2 for angular separation $\gamma = 3\pi/4$. However, as the angular separation approaches $\gamma = \pi$ both the primary and secondary peaks increase their size again and they overcome the noise, thus, leaving the detectors in an entangled state. Gradually, as $\gamma \to \pi$, the primary and secondary peaks superpose and finally create one joint peak aligned at the caustic for $\gamma = \pi$. Here we see that for switching delays that are somewhat larger than the direct null alignment, the extracted entanglement can be dominated by $|M^B_+|$ and thus can be attributed to entanglement harvesting from the field.

The results shown in Fig. 1 are for the field in the Boulware state. As seen in the analogous Fig. 7 for the Unruh state and Fig. 8 for the Hartle-Hawking state, the same phenomena appear in these states as well. In fact, the quantitative differences between the states (which may be difficult to see with the naked eye) are mostly due to the difference in the noise term, which is given in Table I. The correlation terms, in particular in the scenarios of detectors that harvest entanglement in the proximity of caustics, agree to many digits for all three states, as can be seen in Fig. 3, which shows their relative differences.

<table><caption>TABLE I. Numerical values of the local noise terms in the three different states considered, for the scenarios in Fig. 1, Fig. 7 and Fig. 8, respectively.</caption><thead><tr><th>State</th><th>Local noise term</th></tr></thead><tbody><tr><td>Boulware state</td><td>$L_{DD}^B/(\lambda_A\lambda_B) = 9.82 \times 10^{-7}$</td></tr><tr><td>Unruh state</td><td>$L_{DD}^U/(\lambda_A\lambda_B) = 9.96 \times 10^{-7}$</td></tr><tr><td>Hartle-Hawking state</td><td>$L_{DD}^H/(\lambda_A\lambda_B) = 1.08 \times 10^{-6}$</td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>FIG. 4. Entanglement harvesting from the Boulware state by static detectors at antipodal positions close to a Schwarzschild black hole: the plots show the entries of the final detectors' state (28) for two identical detectors ($\lambda_D = 1, \Omega_D = 5/M$) which are placed at equal radial coordinates $r_D$ at angular separation $\gamma = \pi$, i.e., with the black hole exactly in the middle between them. To account for the different redshifts the width of the detector switching functions (25) is set to $T_D = M/\sqrt{1-2M/r}$, i.e., kept equal with respect to the detectors' proper time. The different panels show the results for different switching delays $\Delta_{BA} = t_{0B} - t_{0A}$. Detectors placed at $r = 3M$ are the first for which the switching functions are exactly null aligned at $\Delta_{BA} = 3\sqrt{3}\pi M \approx 16.32M$. Both for larger and lower radial coordinates the light propagation time is longer, and hence two peaks form in the later panels around the radial coordinates for which exact null alignment is achieved (see Fig. 5). The vertical purple lines indicate the radial coordinates at which null alignment of the switching function is achieved, i.e., where $\Delta_{BA} = \Delta t$ (see Fig. 5). (See Supplementary Material for an animated version of this plot [47].)</figcaption></figure>

We conclude this section with a study of entanglement harvesting from the Boulware state between detectors at antipodal locations with the black hole exactly in the middle between them, i.e., at identical radial coordinates $r = r_A = r_A$ and with an angular separation $\gamma = \pi$, for a range of

detector locations reaching down very close to the horizon at radial positions $r = 2.095M$.

At the different radial coordinates the detectors experience different gravitational redshifts. To account for this, we adjust the width of the switching functions (25) to
$T_D(r) = T_\infty / \sqrt{f(r)} = T_\infty / \sqrt{1 - 2M/r}$, such that, at all radial coordinates $T_D(r)$, corresponds to the same amount of proper time given by some value $T_\infty$. Analogous to the plots above, Fig. 4 shows the resulting contributions to the detectors' final state for $T_\infty = 1M$. It shows that at all radial coordinates considered, even close to the horizon, the two detectors can become entangled and the entanglement can be dominated by entanglement harvesting, when the switching delay $\Delta_{BA} = t_{0B} - t_{0A}$ is chosen appropriately.

To understand the structure of the plots and its dependence on $\Delta_{BA}$ it is useful to consider the light propagation coordinate time $\Delta t$ between the two detectors. This is the coordinate time that it takes for a null geodesic starting at spatial coordinates $(r, \theta = \pi/2, \phi)$ to reach the antipodal point $(r, \theta = \pi/2, \phi + \pi)$. Figure 5 shows this time for the range of radial coordinates we consider. It shows that $\Delta t$ is minimal at $r = 3M$. For this critical case, it takes the value $\Delta t_{\min} = 3\sqrt{3\pi}M$. Setting the switching delay equal to the light propagation time, i.e., setting $\Delta_{BA} = \Delta t$, yields exact null alignment of the detectors' switching functions.

Around exact null alignment between the detectors we expect a peak in the final detector entanglement, and this is, indeed, what we observe in Fig. 4. At first, in the panels showing the lower values of $\Delta_{BA}$, one peak forms around $r = 3M$ since the detectors located at this radial coordinate are the first to be exactly null aligned. As $\Delta_{BA}$ is increased in the following panels, a double-peak structure forms since for $\Delta_{BA} > 3\sqrt{3\pi}M$ there are always two radial coordinate positions at which the detectors are exactly null aligned. Regarding the relative size of the contributions $M_+^B$ and $M_-^B$, we observe again that both contributions are of comparable size when the detectors are exactly null aligned, as we already observed in the panel of Fig. 2, for $\gamma = \pi$. This is in contrast to scenarios where the detectors are null aligned along primary null geodesics and sufficiently far away from any caustics. There, just as in flat Minkowski spacetime, the correlations are dominated by the $M_-$ contribution when the detectors are exactly null aligned. However, at the caustics, where a whole envelope of null geodesics connects the detectors at once, the

<figure><img src="image_7.png" /><figcaption>FIG. 5. Light propagation coordinate time between detectors in the scenario of Fig. 4.</figcaption></figure>

singularity structure of the two-point function is altered (see Ref. [37]). This results in the $M_+$ contribution dominating before exact null alignment and the $M_-$ contribution dominating after exact null alignment of the switching functions.

## V. DISCUSSION AND OUTLOOK

Using UDW detectors, we investigated entanglement harvesting from a Klein-Gordon field in the Boulware, Hartle-Hawking and Unruh vacuum states in the background of Schwarzschild spacetime in 3+1 dimensions. We showed that the realistic 3+1-dimensional case possesses a particularly rich phenomenology due to the presence of secondary (and higher) null geodesics and caustics.

In particular, we investigated the ability of two detectors that are static at the same radial coordinate, but with different angular coordinates, to harvest entanglement from the Boulware, Hartle-Hawking and Unruh vacua as they get closer to the horizon.

We paid special attention to entanglement harvesting in the regions of spacetime where secondary null geodesics and caustics can connect the two harvesting detectors. We found that genuine entanglement harvesting can indeed be amplified through “entanglement gravitational lensing” effects when the detectors are close to regions where caustics appear. By the term “genuine entanglement harvesting” we mean the extraction of preexisting entanglement from the field, as opposed to the extraction of entanglement that is created through communication between the detectors. Interestingly, this genuine harvesting can also appear for timelike-separated detectors, for example, when the detectors are connected by secondary null geodesics. Mathematically, this results from a change in the singularity structure of the Wightman function as the field waves cross through caustics.

The formalism that we developed here can also be used to analyze further interesting black hole entanglement harvesting configurations, such as the case where the detectors are in motion and possibly also when crossing the event horizon.

Furthermore, since the effects of entanglement harvesting amplification that we found here are due to lensing which arises from the existence of secondary null geodesics connecting the two detectors, an intriguing follow-up question arises: to what extent could entanglement harvesting be engineered to be amplified even in flat spacetime, namely in the presence of suitable mirror and lensing arrangements? These scenarios will be explored in future work.

## ACKNOWLEDGMENTS

J. G. A. C. acknowledges financial support by CNPq (Brazil), Grant No. 1424412/2019-8. R. H. J. gratefully acknowledges support by the Wenner-Gren Foundations and, in part, by the Wallenberg Initiative on Networks and Quantum Information (WINQ). Nordita is supported in part
by NordForsk. M. C. acknowledges partial financial support by the Scientific Council of the Paris Observatory during a visit. E. M.-M. acknowledges support through the Discovery Grant Program of the Natural Sciences and Engineering Research Council of Canada (NSERC). E. M.-M. also acknowledges support of his Ontario Early Researcher award. A. K. acknowledges support through a Discovery Grant from the National Science and Engineering Research Council of Canada (NSERC) and a Discovery Project Grant from the Australian Research Council (ARC).

## APPENDIX A: DENSITY MATRIX CONTRIBUTIONS

In this appendix, we give a detailed derivation of the expressions used for the calculation of the perturbative contributions to the final density state of the two detectors in Eq. (28). To evaluate the leading-order contributions to the detectors’ final density matrix, we need certain Fourier-type integrals of the switching function (25), $\eta_D(t) = e^{-(t-t_{0D})/T_D}$. First, we have

$$ \int_{-\infty}^{s} dt \, e^{i\nu t} \eta_D(t) = \frac{\sqrt{\pi} T_D}{2} e^{-\frac{\nu^2 T_D^2}{4}} e^{i\nu t_{0D}} \left( 1 + \operatorname{erf} \left( -\frac{i\nu T_D}{2} + \frac{s-t_{0D}}{T_D} \right) \right), \quad (A1) $$

$$ \int_{-\infty}^{\infty} dt \, e^{i\nu t} \eta_D(t) = \sqrt{\pi} T_D e^{-\frac{1}{4}\nu^2 T_D^2} e^{i\nu t_{0D}}, \quad (A2) $$

where $\operatorname{erf}(z) = \frac{z}{\sqrt{\pi}} \int_{0}^{z} dt \, e^{-t^2}$. Furthermore, we need the integral

$$ \int_{-\infty}^{\infty} dt' \, e^{i\mu t'} \eta_D(t') \int_{-\infty}^{t'} dt \, e^{i\nu t} \eta_D(t). \quad (A3) $$

<figure><img src="image_8.png" /><figcaption>FIG. 6. Logarithmic plot version of Fig. 2. (See Supplementary Material for an animated version of this plot [47].)</figcaption></figure>
To this end, consider

$$ \frac{\sqrt{\pi} T_D}{2} e^{-\frac{\nu^2 T_D^2}{4}} e^{i\nu t_0 D} \int_{-\infty}^{\infty} dt' e^{i\mu t'} e^{-\frac{(t'-t_{0D})^2}{T_D^2}} \left( 1 + \operatorname{erf}\left(-\frac{i\nu T_D}{2} + \frac{t' - t_{0D}}{T_D}\right) \right) \\ = \frac{\sqrt{\pi} T_D}{2} e^{-\frac{\nu^2 T_D^2}{4}} e^{i\nu t_0 D} e^{-\frac{t_0^2}{T_D^2}} \int_{-\infty}^{\infty} dt' e^{-\frac{t'^2}{T_D^2} + t'(2^{(t_0')/2} + i\mu)} \operatorname{erf}\left(-\frac{i\nu T_D}{2} + \frac{t' - t_{0D}}{T_D}\right) + \frac{\pi T_D T_{D'}}{2} e^{-\frac{\mu^2 T_D^2 + i\nu^2 T_D^2}{4}} e^{i(\mu t_{0D} + i\nu t_0 D)}. \quad (A4) $$

For the $t'$ integral, use Eqs. (A3) and (A7) from Ref. [44], which is

$$ I(a, b) = \int_{-\infty}^{\infty} dy e^{-a^2 - iby - y^2} \operatorname{erf}(y - ia) = -i\sqrt{\pi}e^{-a^2 - b^2/4} \operatorname{erfi}\left(\frac{a + b/2}{\sqrt{2}}\right), \quad (A5) $$

where $\operatorname{erfi}(z) = -i\operatorname{erf}(iz)$, with $y = t'/T_{D'}$, $a = \nu T_D/2 - i t_{0D}/T_D$ and $b = 2i t_{0D}/T_D - \mu T_{D'}$; then

$$ \int_{-\infty}^{\infty} dt' e^{-\frac{t'^2}{T_D^2} + t'(2^{(t_0')/2} + i\mu)} \operatorname{erf}\left(-\frac{i\nu T_D}{2} + \frac{t' - t_{0D}}{T_D}\right) = -i\sqrt{\pi} T_{D'} e^{\frac{-\mu^2 T_D^2}{4} + \frac{t_0^2}{T_D^2} + i\mu t_{0D}'} \operatorname{erfi}\left(\frac{1}{\sqrt{2}}\left(\frac{\nu T_D - \mu T_{D'}}{2} - i\frac{t_{0D}}{T_D} + i\frac{t_{0D}'}{T_{D'}}\right)\right). \quad (A6) $$

Inserting this above yields

$$ \int_{-\infty}^{\infty} dt' e^{i\mu t'} \eta_{D'}(t') \int_{-\infty}^{t'} dt e^{i\nu t'} \eta_{D}(t) = \frac{T_D T_{D'}}{2\pi} e^{\frac{\mu^2 T_D^2}{4} + \frac{t_0^2}{T_D^2}} e^{i(\nu t_0 D + \mu t_{0D}')} \left( 1 - \operatorname{ierfi}\left[ \frac{\nu T_D - \mu T_{D'}}{2\sqrt{2}} + \frac{i}{\sqrt{2}} \left( \frac{t_{0D}'}{T_{D'}} - \frac{t_{0D}}{T_D} \right) \right] \right). \quad (A7) $$

For stationary detectors we use $\tau_D(t) = \sqrt{f(r_D)}t$ as the relation between proper time and coordinate time [see Eq. (24)]. With this, and Eqs. (10) and (A2), the single detector noise term $L_{DD'}^{\Psi}$ in Eq. (30) takes the following form [denoting $r = r_D, r' = r_{D'}, f_D = f(r_D)$ and $N_D = \sqrt{f_D}$]:

$$ \begin{aligned} L_{DD'}^{\Psi} &= \lambda_D \lambda_{D'} f_A f_B \int_{-\infty}^{\infty} dt \int_{-\infty}^{\infty} dt' \eta_D(t) \eta_{D'}(t') e^{-i\Omega_D N_D t + i\Omega_{D'} N_{D'} t'} W^{\Psi}(\mathbf{x}_D(t); \mathbf{x}_{D'}(t')) \\ &= \frac{\lambda_D \lambda_{D'} f_A f_B}{(4\pi)^2} \sum_{\ell=0}^{\infty} (2\ell + 1) P_{\ell}(\cos \gamma) \int_{-\infty}^{\infty} d\omega \frac{G_{\ell\omega}^{\Psi}(r_D, r_{D'})}{\omega} \left( \int_{-\infty}^{\infty} dt \eta_D(t) e^{-i(\Omega_D N_D + \omega)t} \right) \left( \int_{-\infty}^{\infty} dt' \eta_{D'}(t') e^{i(\Omega_{D'} N_{D'} + \omega)t'} \right) \\ &= \frac{\lambda_D \lambda_{D'} N_D N_{D'} T_D T_{D'}}{16\pi} \sum_{\ell=0}^{\infty} (2\ell + 1) P_{\ell}(\cos \gamma) \int_{-\infty}^{\infty} d\omega \frac{G_{\ell\omega}^{\Psi}(r_D, r_{D'})}{\omega} e^{i(\Omega_{D'} N_{D'} + \omega)t_{0D'} - i(\Omega_D N_D + \omega)t_{0D} - \frac{1}{4}(\Omega_D N_D + \omega)^2 T_D^2 - \frac{1}{4}(\Omega_{D'} N_{D'} + \omega)^2 T_{D'}^2} \end{aligned} \quad (A8) $$

Analogously, for $M^{\Psi}$ in Eq. (29) we obtain using Eq. (A7),

$$ \begin{aligned} M^{\Psi} &= \frac{-\lambda_{\Lambda}\lambda_{\Lambda'}f_{\Lambda}f_{\Lambda'}}{(4\pi)^2} \sum_{\ell=0}^{\infty} (2\ell + 1) P_{\ell}(\cos\gamma) \int_{-\infty}^{\infty} d\omega \frac{1}{\omega} \left( \int_{-\infty}^{\infty} dt e^{i(\Omega_{\Lambda}f_{\Lambda}-\omega)t} \eta_{\Lambda}(t) \int_{-\infty}^{t} dt' e^{i(\omega+\Omega_{\Lambda}f_{\Lambda})t'} \eta_{\Lambda'}(t') G_{\ell\omega}^{\Psi}(r_{\Lambda}, r_{\Lambda'}) \right. \\ & \qquad \left. + \int_{-\infty}^{\infty} dt \eta_{\Lambda}(t) e^{i(\Omega_{\Lambda}f_{\Lambda}-\omega)t} \int_{-\infty}^{t} dt' e^{i(\omega+\Omega_{\Lambda}f_{\Lambda})t'} \eta_{\Lambda'}(t') G_{\ell\omega}^{\Psi}(r_{\Lambda}, r_{\Lambda}) \right) \\ &= \frac{-\lambda_{\Lambda}\lambda_{\Lambda'}f_{\Lambda}f_{\Lambda'}T_{\Lambda}T_{\Lambda'}}{32\pi} \sum_{\ell=0}^{\infty} (2\ell + 1) P_{\ell}(\cos\gamma) \int_{-\infty}^{\infty} d\omega \frac{1}{\omega} \\ & \qquad \times \left( e^{-\frac{\mu_{\Lambda}^2 T_{\Lambda}^2 + i\mu_{\Lambda}^2 T_{\Lambda}^2}{4}} e^{i(\nu_{\Lambda}t_{0\Lambda} + \mu_{\Lambda}t_{0\Lambda})} \left( 1 - \operatorname{ierfi}\left[ \frac{\nu_{\Lambda}T_{\Lambda} - \mu_{\Lambda}T_{\Lambda}}{2\sqrt{2}} + \frac{i}{\sqrt{2}} \left( \frac{t_{0\Lambda}}{T_{\Lambda}} - \frac{t_{0\Lambda}}{T_{\Lambda}} \right) \right] \right) G_{\ell\omega}^{\Psi}(r_{\Lambda}, r_{\Lambda}) \right. \\ & \qquad \left. + e^{-\frac{\mu_{\Lambda}^2 T_{\Lambda}^2 + i\mu_{\Lambda}^2 T_{\Lambda}^2}{4}} e^{i(\nu_{\Lambda}t_{0\Lambda} + \mu_{\Lambda}t_{0\Lambda})} \left( 1 - \operatorname{ierfi}\left[ \frac{\nu_{\Lambda}T_{\Lambda} - \mu_{\Lambda}T_{\Lambda}}{2\sqrt{2}} + \frac{i}{\sqrt{2}} \left( \frac{t_{0\Lambda}}{T_{\Lambda}} - \frac{t_{0\Lambda}}{T_{\Lambda}} \right) \right] \right) G_{\ell\omega}^{\Psi}(r_{\Lambda}, r_{\Lambda}) \right), \end{aligned} \quad (A9) $$
<figure><img src="image_9.png" /><figcaption>FIG. 7. Gravitational lensing of entanglement harvesting from the field in the Unruh state: the entire setup is identical to the setup of Fig. 1, except that here the initial state for the field is the Unruh state. The plots show the absolute values of $M^U$, $M^U_\pm$ and $L^U_{DD}$ in Eq. (28) for two static detectors placed at radial coordinate $r = 6.009M$ with varying angular separation $\gamma$. All detector parameters are equal ($\lambda_D = 1, \Omega_D = 5M^{-1}, T_D = 1M$); only the offset $\Delta_{BA} = t_{0B} - t_{0A}$ between the two switching functions (25) varies between the plots. The inset indicates how far a null wave front propagates from the red point within a coordinate time interval $\Delta_{BA}$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 8. Gravitational lensing of entanglement harvesting from the field in the Hartle-Hawking state: the entire setup is identical to the setup of Fig. 1, except that here the initial state for the field is the Hartle-Hawking state. The plots show the absolute values of $M^H$, $M^\pm$ and $L_{DD}^H$ in Eq. (28) for two static detectors placed at radial coordinate $r = 6.009M$ with varying angular separation $\gamma$. All detector parameters are equal ($\lambda_D = 1, \Omega_D = 5M^{-1}, T_D = 1M$); only the offset $\Delta_{BA} = t_{0B} - t_{0A}$ between the two switching functions (25) varies between the plots. The inset indicates how far a null wave front propagates from the red point within a coordinate time interval $\Delta_{BA}$.</figcaption></figure>
with $\nu_D = \omega + \Omega_D N_D$, $\mu_D = \Omega_D N_D - \omega$. For practical evaluation purposes, one can manipulate the integrals in Eqs. (A9) and (A8) further so that the integration over $\omega$ only runs over $0 < \omega < \infty$.

## APPENDIX B: NUMERICAL TECHNIQUES FOR THE WIGHTMAN FUNCTION

In order to evaluate the Wightman function in the situations presented in this work, integrals over $R_{\ell\omega}^{\text{in/up}}(r)$ with respect to $r$ are needed. However, there is no known closed-form expression for $R_{\ell\omega}^{\text{in/up}}(r)$. Hence, one has to use numerical techniques. Specifically, we use numerical methods to evaluate $R_{\ell\omega}^{\text{in}}$, $R_{\ell\omega}^{\text{up}}$ and the set of coefficients $I_{\ell\omega}$, $\rho_{\ell\omega}^{\text{in/up}}$ in the region $M\omega \in [0, 10]$ in steps of $10^{-3}$, $r^*/M \in [-4, 13]$ in steps of $1/5$ and $\ell$ for all (integer) values $0 \le \ell \le 100$. Different techniques are involved in the evaluation.

For $R_{\ell\omega}^{\text{up}}$ we used the NDSolve-based numerical method implemented in the Regge-Wheeler Mathematica [48] package of the Black Hole Perturbation Toolkit [49].

For $R_{\ell\omega}^{\text{in}}$ we define two auxiliary quantities $\bar{r} = r/2M$ and $\bar{\omega} = 2M\omega$, then use the so-called Jaffé series [50] to write the $R_{\ell\omega}^{\text{in}}$ solution to (4) for $\bar{r} \in [1, \infty)$ as

$$R_{\ell\omega}^{\text{in}}(\bar{r}) = \bar{r}^{2i\bar{\omega}}(\bar{r}-1)^{-i\bar{\omega}}e^{i\bar{\omega}\bar{r}} \sum_{n=0}^{\infty} a_n \left(\frac{\bar{r}-1}{\bar{r}}\right)^n, \quad (B1)$$

where

$$\begin{align*}
& a_n\alpha_{n-1} + a_{n-1}\beta_{n-1} + a_{n-2}\gamma_{n-1} = 0, \\
& \alpha_n = (n+1)(n+1-2i\bar{\omega}), \\
& \beta_n = -1-2n(n+1)-\ell(\ell+1)+4\bar{\omega}(i+2ni+2\bar{\omega}), \\
& \gamma_n = (n-2i\bar{\omega})^2, \tag{B2}
\end{align*}$$

with $a_0 = e^{-2i\bar{\omega}}$ and $a_n = 0 \ \forall \ n < 0$. This solution is implemented in Mathematica. However some care has to be taken: in practice we have to impose a cutoff on the number of terms $n_{\text{max}}$ to include in the series. To set an adequate value for $n_{\text{max}}$ we used the following quantity:

$$\varepsilon(\ell, \bar{\omega}, \bar{r}, n_{\text{max}}) = \frac{a_{n_{\text{max}}+1}(\frac{\bar{r}-1}{\bar{r}})^{n_{\text{max}}+1}}{\sum_{n=0}^{n_{\text{max}}} a_n(\frac{\bar{r}-1}{\bar{r}})^n}. \quad (B3)$$

Since we want to have at least 16 digits of precision in the $R_{\ell\omega}^{\text{in}}$ modes, we want to find an $n_{\text{max}}$ such that

$$\begin{align}
\operatorname{Re}[\varepsilon(\ell, \bar{\omega}, \bar{r}, n_{\text{max}})] &< 10^{-16} \quad \text{and} \nonumber \\
\operatorname{Im}[\varepsilon(\ell, \bar{\omega}, \bar{r}, n_{\text{max}})] &< 10^{-16} \tag{B4}
\end{align}$$

for all values of $\bar{r}, \ell$ and $\bar{\omega}$ in the region where the solutions are evaluated. For that, we use the fact that the convergence of the series becomes slower with increasing $\bar{r}$ to conclude that the largest $\varepsilon$ happens at the highest $\bar{r}$ we intended to evaluate, which is around $\bar{r}_{\text{max}} \approx 5$. Then, fixing $\bar{r} = \bar{r}_{\text{max}}$ we test a couple of $\bar{\omega}$ and $\ell$ values in the region of interest and conclude that $n_{\text{max}} = 5000$ is enough to satisfy Eq. (B4) in all regions of interest.

For the $I_{\ell\omega}$, $\rho_{\ell\omega}^{\text{in/up}}$ coefficients we use data from Ref. [36], which is evaluated using a Mathematica implementation of the Mano-Suzuki-Takasugi method. An extensive review of that method can be found in Ref. [51].

### 1. $L_{DD}^{\Psi}$ terms

By setting $D = D'$ in Eq. (A8) and manipulating the integration range from $-\infty < \omega < \infty$ to $0 < \omega < \infty$, we obtain the following expression for the $L_{DD}^{\Psi}$ terms:

$$L_{DD}^{\Psi}(r_D) = \frac{(N_D \lambda_D T_D)^2}{16\pi} \sum_{\ell=0}^{\infty} (2\ell + 1) \int_0^{\infty} d\omega \frac{1}{\omega} \left( G_{\ell\omega}^{\Psi}(r_D, r_D) e^{-\frac{1}{2}(\Omega_D N_D + \omega)^2 T_D^2} - G_{\ell,-\omega}^{\Psi}(r_D, r_D) e^{-\frac{1}{2}(\Omega_D N_D - \omega)^2 T_D^2} \right). \quad (B5)$$

To verify the convergence of such an integral, one should check the large-$|\omega|$ behavior of the $G_{\ell\omega}^{\Psi}(r_D, r_D)$, which varies among the quantum states (11)–(13):

$$\begin{align}
& G_{\ell\omega}^{B}(r_D, r_D) \sim \theta(\omega)(|\bar{R}_{\ell\omega}^{\text{up}}(r_D)|^2 + |\bar{R}_{\ell\omega}^{\text{in}}(r_D)|^2), \nonumber \\
& G_{\ell\omega}^{U}(r_D, r_D) \sim |\bar{R}_{\ell\omega}^{\text{up}}(r_D)|^2(\theta(\omega) - \theta(-\omega)e^{2\pi\omega\kappa}) + \theta(\omega)|\bar{R}_{\ell\omega}^{\text{in}}(r_D)|^2, \nonumber \\
& G_{\ell\omega}^{H}(r_D, r_D) \sim |\bar{R}_{\ell\omega}^{\text{up}}(r_D)|^2(\theta(\omega) - \theta(-\omega)e^{2\pi\omega\kappa}) + |\bar{R}_{\ell\omega}^{\text{in}}(r_D)|^2(\theta(\omega)e^{-2\pi\omega\kappa} - \theta(-\omega)). \tag{B6}
\end{align}$$
<figure><img src="image_11.png" /></figure>

FIG. 9. (a) Integrand in Eq. (B5) as a function of $M\omega$ and $\ell$ for the Boulware state. The value of $M\omega$ where the asymptotic regime begins grows with $\ell$. More precisely, the transition to the asymptotic regime begins around $\omega r \sim \sqrt{\ell(\ell+1)}$ (for $\ell > 0$ and radii not too close to the horizon). Notice that thanks to the $\theta(\omega)$ in Eq. (11), the slowest decaying term in Eq. (B5) for the Boulware state is proportional to $\sim 2 \frac{e^{-\frac{1}{2}(\Omega_D N_D + \omega)^2 T_D^2}}{\omega}$. (b) Relative difference between the integral in Eq. (B5) for the Boulware state when integrated up to $\omega_{cut}$ and up to $\omega_{cut} + \delta$, with $\delta = 1/10$. The reason to limit this plot to $10^{-10}$ is that those integrals where evaluated to ten significant digits; hence, a relative difference of $10^{-10}$ means that the integral converged completely to all significant digits. Comparing with Fig. 9(a), we see that exponentially fast convergence happens after $M\omega_{cut}$ becomes slightly larger than the $M\omega$ where the asymptotic regime begins. (c) Summand in Eq. (B5) as a function of $\ell$ for the Boulware state. One can notice from this figure that the summand decays superexponentially with increasing $\ell$. In fact, for $\ell > 30$, it has already converged to all ten significant digits that we got when evaluating the integrals. (d) Relative difference between the $\ell$ sum in Eq. (B5) for the Boulware state when summed up to $l_{cut}$ and $l_{cut} - 1$. The reason to limit this plot to $10^{-10}$ is that the integrals in the sums were to 10 digits of precision. Hence, a relative difference that is smaller than $10^{-10}$ should represent numerical noise and not actual significant digits.

The leading asymptotic behavior of $|\tilde{R}_{\ell\omega}^{\text{in/up}}|^2$ as $|\omega| \to \infty$ is $|\tilde{R}_{\ell\omega}^{\text{in/up}}|^2 \sim 1$ [52].<fnref n="9" /> Hence, to leading order for large $|\omega|$, $G_{\ell\omega}^{\Psi}(r_D, r_D) \sim 1$. By substituting such a result into the integrand in Eq. (B5), one concludes that the slowest-decaying term of the integrand falls off as $\sim 2 \frac{e^{-\frac{1}{2}(\Omega_D N_D - \omega)^2 T_D^2}}{\omega}$. Therefore, the integral in Eq. (B5) is fast-converging. Given that, we are able to numerically evaluate it accurately enough with a frequency cutoff $M\omega_{cut} = 10$. The

asymptotic regime of the integrand and the convergence of the integral are illustrated in Figs. 9(a) and 9(b).

In order to verify the convergence of the $\ell$ sum, we begin by evaluating the integrals in Eq. (B5), multiplying by $2\ell + 1$ and plotting the result against $\ell$. The outcome is shown in Fig. 9(c), where we can see that such a quantity decays superexponentially. Given this numerical evidence, the $\ell$ sum in that equation is expected to converge with a good accuracy even with a cutoff lower than the $\ell_{cut} = 100$ we use. In fact, as presented in Fig. 9(d), $\ell_{cut} = 30$ is enough to converge to all ten significant digits used when evaluating the Boulware integrals.

<footnote n="9">⁹Such asymptotics are nonuniform and valid only when $\omega \gg V_\ell(r)$, as illustrated in Fig. 9(a).</footnote>
## 2. M^Ψ terms

Since we already know the large-|ω| asymptotics of $G_{\ell\omega}^{\Psi}(r_D, r_D)$, to analyze the behavior of the integrand in Eq. (A9) in that regime we have to study the factors that multiply each of its terms. These factors are very similar and it suffices to analyze only one of them:

$$e^{-\frac{\mu_A^2 T_A^2 + \ell_B^2 T_B^2}{4}} e^{i(\nu_B t_{0B} + \mu_A t_{0A})} \text{erfc} \left[ \frac{i(\nu_B T_B - \mu_A T_A)}{2\sqrt{2}} - \frac{1}{\sqrt{2}} \left( \frac{t_{0A}}{T_A} - \frac{t_{0B}}{T_B} \right) \right], \quad (B7)$$

where $\text{erfc}(z) = 1 - \text{erf}(z)$. Expanding the exponential functions, we obtain

$$\alpha_{AB} e^{-\frac{1}{4}(T_A^2 + T_B^2)\omega^2} e^{-\frac{1}{2}(N_B T_B^2 \Omega_B - N_A T_A^2 \Omega_A)\omega} e^{i\omega(t_{0B} - t_{0A})} \text{erfc} \left[ \frac{i(\nu_B T_B - \mu_A T_A)}{2\sqrt{2}} - \frac{1}{\sqrt{2}} \left( \frac{t_{0A}}{T_A} - \frac{t_{0B}}{T_B} \right) \right], \quad (B8)$$

where

$$\alpha_{AB} = e^{-\frac{1}{4}(N_B^2 T_B^2 \Omega_B^2 + N_A^2 T_A^2 \Omega_A^2)} e^{i(\Omega_B N_B t_{0B} + \Omega_A N_A t_{0A})}. \quad (B9)$$

Now we use an asymptotic expression for the erfc(z),

$$\text{erfc}(z) \sim \frac{e^{-z^2}}{z\sqrt{\pi}}, \quad z \to \infty, \quad |\text{ph}(z)| < 3\pi/4, \quad (B10)$$

to conclude that as $\omega \to \infty$,

$$\text{erfc} \left[ \frac{i(\nu_B T_B - \mu_A T_A)}{2\sqrt{2}} - \frac{1}{\sqrt{2}} \left( \frac{t_{0A}}{T_A} - \frac{t_{0B}}{T_B} \right) \right] \sim \beta_{AB} \frac{e^{\frac{1}{8}(T_A + T_B)^2 \omega^2} e^{\frac{1}{4}(T_A + T_B)(N_B T_B \Omega_B - N_A T_A \Omega_A)\omega} e^{-\frac{i\omega}{2}(T_A + T_B)\left(\frac{t_{0B}}{T_B} - \frac{t_{0A}}{T_A}\right)}}{\frac{i((N_B \Omega_B + \omega)T_B - (N_A \Omega_A - \omega)T_A)}{2\sqrt{2}} - \frac{1}{\sqrt{2}}\left(\frac{t_{0A}}{T_A} - \frac{t_{0B}}{T_B}\right)}, \quad (B11)$$

where

$$\beta_{AB} = (\pi\alpha_{AB})^{-1/2} e^{-\frac{1}{2}(\frac{t_{0B}}{T_B} - \frac{t_{0A}}{T_A})^2} e^{\frac{i}{2}(N_A T_A \Omega_A + N_B T_B \Omega_B)\frac{t_{0B}}{T_B} + N_A T_A \Omega_A} e^{-\frac{1}{4}(N_A T_A \Omega_A + N_B T_B \Omega_B)} e^{-\frac{1}{4}(N_A T_A \Omega_A + N_B T_B \Omega_B)}. \quad (B12)$$

Putting everything together we conclude that

$$e^{-\frac{\mu_A^2 T_A^2 + \ell_B^2 T_B^2}{4}} e^{i(\nu_B t_{0B} + \mu_A t_{0A})} \text{erfc} \left[ \frac{i(\nu_B T_B - \mu_A T_A)}{2\sqrt{2}} - \frac{1}{\sqrt{2}} \left( \frac{t_{0A}}{T_A} - \frac{t_{0B}}{T_B} \right) \right] \\ \sim \alpha_{AB} \beta_{AB} \frac{e^{-\frac{1}{8}(T_A - T_B)\omega^2} e^{\frac{1}{4}(T_A - T_B)(N_B T_B \Omega_B + N_A T_A \Omega_A)\omega} e^{\frac{i}{2}\omega(T_B - T_A)\left(\frac{t_{0A}}{T_A} + \frac{t_{0B}}{T_B}\right)}}{\frac{i((N_B \Omega_B + \omega)T_B - (N_A \Omega_A - \omega)T_A)}{2\sqrt{2}} - \frac{1}{\sqrt{2}}\left(\frac{t_{0A}}{T_A} - \frac{t_{0B}}{T_B}\right)}, \quad \omega \to \infty, \quad (B13)$$

which leads to the conclusion that, whenever $T_A \neq T_B$, the whole term decays superexponentially as $e^{-\frac{1}{8}(T_B - T_A)^2 \omega^2}$. On the other hand, if $T_A = T_B \equiv T$, which is the case we study in this work, Eq. (B13) simplifies to

$$e^{-\frac{\mu_A^2 + \ell_B^2}{4} T^2} e^{i(\nu_B t_{0B} + \mu_A t_{0A})} \text{erfc} \left[ \frac{i(\nu_B - \mu_A)T}{2\sqrt{2}} + \frac{1}{T\sqrt{2}}(t_{0B} - t_{0A}) \right] \sim \frac{\alpha_{AB} \beta_{AB}}{\frac{i(N_B \Omega_B - N_A \Omega_A + 2\omega)T}{2\sqrt{2}} - \frac{(t_{0A} - t_{0B})}{T\sqrt{2}}}, \quad \omega \to \infty. \quad (B14)$$

Then both, the real and imaginary parts of this quantity decay linearly with $\omega$, since the absolute value of the denominator grows linearly with the frequency. Given that, by considering the extra $\omega^{-1}$ in the integrand in Eq. (A9), we conclude that it decays as $\omega^{-2}$. Hence, the integrals defining the M^Ψ terms are convergent.
<figure><img src="image_12.png" /><figcaption>FIG. 10. (a) Boulware state integrand in Eq. (A9) after rewriting the integral so that it ranges from $\omega = 0$ to $M\omega = 10$, as a function of $M\omega$, for several different $\ell = 0, 10, 20, 30$. (b) Relative difference in the Boulware state integral in Eq. (A9) after rewriting the integral so that it ranges from $\omega = 0$ to $\omega = \omega_{\text{cut}}$ for several different $\ell = 0, 10, 20, 30$. To evaluate the relative difference, we integrate up to $\omega_{\text{cut}}$ and then to $\omega_{\text{cut}} + \delta$, with $\delta = 1/10$. (c) Boulware state integral in Eq. (A9), after rewriting the integral so that it ranges from $\omega = 0$ to $\omega = \omega_{\text{cut}}$ for several different $\ell = 0, 10, 20, 30$, here integrated up to $M\omega_{\text{cut}} = 10$. (d) Relative difference in the Boulware state integral from Eq. (A9), after rewriting the integral so that it ranges from $\omega = 0$ to $\omega = \omega_{\text{cut}}$ for several different $\ell = 0, 10, 20, 30$. To evaluate the relative difference, we sum up to $\ell_{\text{cut}}$ and then to $\ell_{\text{cut}} + 1$, ranging from $\ell_{\text{cut}} = 1$ up to $\ell_{\text{cut}} = 100$. For $\ell_{\text{cut}} > 45$ all significant digits exactly cancel out when evaluating the relative difference.</figcaption></figure>

We remark that such leading asymptotic behavior in Eq. (B14) is, in general, not achieved by $M\omega = 10$, but one can still perform the numerical integration to a good accuracy because, before this leading asymptotic regime, there is an intermediary regime where the integrand is exponentially decaying, as can be seen in Fig. 10(a). We can further confirm that convergence by looking at the

integral in Eq. (A9) as a function of $\omega_{\text{cut}}$: as can be seen in Fig. 10(b), up to $\ell = 20$ we have convergence up to four significant digits, while for larger values of $\ell$, the convergence becomes worse, resulting in no significant digits at all. Yet, since the resulting integral decays superexponentially with $\ell$, as presented in Fig. 10(c), one can still obtain up to four significant digits, as can be seen in Fig. 10(d).
[1] S. J. Summers and R. Werner, The vacuum violates Bell's inequalities, Phys. Lett. **110A**, 257 (1985).

[2] S. J. Summers and R. Werner, Bell's inequalities and quantum field theory. I. General setting, J. Math. Phys. (N.Y.) **28**, 2440 (1987).

[3] M. Saravani, S. Aslanbeigi, and A. Kempf, Spacetime curvature in terms of scalar field propagators, Phys. Rev. D **93**, 045026 (2016).

[4] A. Kempf, Replacing the notion of spacetime distance by the notion of correlation, Front. Phys. **9**, 655857 (2021).

[5] T. R. Perche and E. Martín-Martínez, Geometry of spacetime from quantum measurements, Phys. Rev. D **105**, 066011 (2022).

[6] J. D. Bekenstein, Black holes and entropy, Phys. Rev. D **7**, 2333 (1973).

[7] L. Bombelli, R. K. Koul, J. Lee, and R. D. Sorkin, Quantum source of entropy for black holes, Phys. Rev. D **34**, 373 (1986).

[8] M. Srednicki, Entropy and Area, Phys. Rev. Lett. **71**, 666 (1993).

[9] J. Eisert, M. Cramer, and M. B. Plenio, Colloquium: Area laws for the entanglement entropy, Rev. Mod. Phys. **82**, 277 (2010).

[10] E. Witten, APS medal for exceptional achievement in research: Invited article on entanglement properties of quantum field theory, Rev. Mod. Phys. **90**, 045003 (2018).

[11] S. W. Hawking, Black hole explosions, Nature (London) **248**, 30 (1974).

[12] S. W. Hawking, Particle creation by black holes, Commun. Math. Phys. **43**, 199 (1975).

[13] P. C. W. Davies, Scalar particle production in Schwarzschild and Rindler metrics, J. Phys. A **8**, 609 (1975).

[14] D. N. Page, Information in Black Hole Radiation, Phys. Rev. Lett. **71**, 3743 (1993).

[15] D. N. Page, Time dependence of Hawking radiation entropy, J. Cosmol. Astropart. Phys. **09** (2013) 028.

[16] G. Penington, Entanglement wedge reconstruction and the information paradox, J. High Energy Phys. **02** (2020).002

[17] B. Reznik, A. Retzker, and J. Silman, Violating Bell's inequalities in vacuum, Phys. Rev. A **71**, 042104 (2005).

[18] A. Pozas-Kerstjens and E. Martín-Martínez, Harvesting correlations from the quantum vacuum, Phys. Rev. D **92**, 064042 (2015).

[19] A. Valentini, Non-local correlations in quantum electrodynamics, Phys. Lett. A **153**, 321 (1991).

[20] B. Reznik, Entanglement from the vacuum, Found. Phys. **33**, 167 (2003).

[21] G. VerSteeg and N. C. Menicucci, Entangling power of an expanding universe, Phys. Rev. D **79**, 044027 (2009).

[22] E. Martín-Martínez, A. R. H. Smith, and D. R. Terno, Spacetime structure and vacuum entanglement, Phys. Rev. D **93**, 044001 (2016).

[23] L. J. Henderson, R. A. Hennigar, R. B. Mann, A. R. H. Smith, and J. Zhang, Harvesting entanglement from the black hole vacuum, Classical Quantum Gravity **35**, 21LT02 (2018).

[24] K. Gallock-Yoshimura, E. Tjoa, and R. B. Mann, Harvesting entanglement with detectors freely falling into a black hole, Phys. Rev. D **104**, 025001 (2021).

[25] M. Cliche and A. Kempf, Vacuum entanglement enhancement by a weak gravitational field, Phys. Rev. D **83**, 045019 (2011).

[26] R. H. Jonsson, D. Q. Aruquipa, M. Casals, A. Kempf, and E. Martín-Martínez, Communication through quantum fields near a black hole, Phys. Rev. D **101**, 125005 (2020).

[27] P. Candelas, Vacuum polarization in Schwarzschild spacetime, Phys. Rev. D **21**, 2185 (1980).

[28] D. G. Boulware, Spin-1 2 quantum field theory in Schwarzschild space, Phys. Rev. D **12**, 350 (1975).

[29] W. G. Unruh, Notes on black-hole evaporation, Phys. Rev. D **14**, 870 (1976).

[30] J. B. Hartle and S. W. Hawking, Path-integral derivation of black-hole radiance, Phys. Rev. D **13**, 2188 (1976).

[31] S. Hollands and R. M. Wald, Quantum fields in curved spacetime, Phys. Rep. **574**, 1 (2015), quantum fields in curved spacetime.

[32] B. S. DeWitt and R. W. Brehme, Radiation damping in a gravitational field, Ann. Phys. (N.Y.) **9**, 220 (1960).

[33] J. Hadamard, *Lectures on Cauchy's Problem in Linear Partial Differential Equations* (Dover Publications, New York, 1923).

[34] P. R. Garabedian, *Partial Differential Equations* (Chelsea Pub Co, New York, 1998).

[35] M. Ikawa, *Hyperbolic Partial Differential Equations and Wave Phenomena. Iwanami Series in Modern Mathematics*, Translations of Mathematical Monographs (American Mathematical Soc., Providence, 2000).

[36] C. Buss and M. Casals, Quantum correlator outside a Schwarzschild black hole, Phys. Lett. B **776**, 168 (2018).

[37] M. Casals and B. Nolan, Global Hadamard form for the green function in Schwarzschild space-time, arXiv:1606.03075.

[38] B. S. DeWitt, Quantum gravity: The new synthesis, in *General Relativity: An Einstein Centenary Survey*, edited by S. Hawking and W. Israel (Cambridge University Press, Cambridge England, New York, 1979), p. 680.

[39] E. Martín-Martínez, T. R. Perche, and B. d. S. L. Torres, Broken covariance of particle detector models in relativistic quantum information, Phys. Rev. D **103**, 025007 (2021).

[40] J. de Ramón, M. Papageorgiou, and E. Martín-Martínez, Relativistic causality in particle detector models: Faster-than-light signaling and impossible measurements, Phys. Rev. D **103**, 085002 (2021).

[41] E. Martín-Martínez, Causality issues of particle detector models in QFT and quantum optics, Phys. Rev. D **92**, 104019 (2015).

[42] E. Martín-Martínez, T. R. Perche, and B. de S. L. Torres, General relativistic quantum optics: Finite-size particle detector models in curved spacetimes, Phys. Rev. D **101**, 045017 (2020).

[43] R. Lopp and E. Martín-Martínez, Quantum delocalization, gauge, and quantum optics: Light-matter interaction in relativistic quantum information, Phys. Rev. A **103**, 013703 (2021).

[44] A. Pozas-Kerstjens and E. Martin-Martinez, Harvesting correlations from the quantum vacuum, Phys. Rev. D **92**, 064042 (2015).

[45] E. Martin-Martinez and P. Rodriguez-Lopez, Relativistic quantum optics: On the relativistic invariance of the light-matter interaction models, Phys. Rev. D **97**, 105026 (2018).
[46] E. Tjoa and E. Martín-Martínez, When entanglement harvesting is not really harvesting, Phys. Rev. D **104**, 125005 (2021).

[47] See Supplemental Material at http://link.aps.org/supplemental/10.1103/PhysRevD.108.025016 for an animated version of Fig. 1.

[48] Wolfram Research, Inc., Mathematica, Version 13.2, Campaign, IL, 2022.

[49] Black Hole Perturbation Toolkit (bhptoolkit.org).

[50] E. W. Leaver, Solutions to a generalized spheroidal wave equation: Teukolsky's equations in general relativity, and the twocenter problem in molecular quantum mechanics, J. Math. Phys. (N.Y.) **27**, 1238 (1986).

[51] M. Sasaki and H. Tagoshi, Analytic black hole perturbation approach to gravitational radiation, Living Rev. Relativity **6**, 6 (2003).

[52] A. Levi and A. Ori, Mode-sum regularization of $\phi^2$ in the angular-splitting method, Phys. Rev. D **94**, 044054 (2016).