Lensing of Vacuum Entanglement near Schwarzschild Black Holes

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 pre-existing vacuum entanglement can be harvested also for lightlike separations.


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][7][8][9][10][11][12][13][14][15][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 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 arXiv:2303.01402v3 [quant-ph] 10 Aug 2023 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 preexisting 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.
Sec. 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. Sec. III introduces the detector model, its perturbative treatment and negativity as entanglement measure for the detectors' final state. Sec. III B discusses when the entanglement between detectors can be attributed to the harvesting of preexisting entanglement from the field. Sec. 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, App. A gives the calculations of the perturbative contributions to the detectors' state and App. 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 = ℏ = 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 ∈ R, r ∈ (2M, ∞), θ ∈ [0, π], φ ∈ [0, 2π)}, is given by 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 ϕ obeying the Klein-Gordon (K-G) equation □φ = 0. Considering the Schwarzschild spacetime and availing of its spherical symmetry, a general real-valued solution for that equation can be written aŝ where x denotes a spacetime point,α in/up † ℓmω are creation and annihilation operators and with Y ℓm being the spherical harmonic of degree ℓ and order m and R in/up ℓω are radial factors. Substituting the field modes (3) into the K-G equation leads to the conclusion that R in/up ℓω where is an effective potential.
Here, ρ in/up ℓω ∈ C are the reflection amplitudes and I ℓω ∈ C is the incidence amplitude. For compatibility with Ref. [27] and ease of notation, let us also definē Note that it suffices to calculate the modes for just ω ≥ 0 and use the following symmetries for ω < 0: Calculating the response of the detectors to the interaction with the quantized fieldφ requires, as we shall see in the next section, the Wightman function. The Wightman function, when the quantum field is in a state |Ψ⟩, is defined as It is thus a two-point function satisfying the homogeneous K-G equation. Henceforth we shall only consider quantum states |Ψ⟩ in regions of spacetime where they satisfy the Hadamard property (namely, Eq. (14), is satisfied for the Wightmant function in such states). The Wightman function in Schwarzschild spacetime is given by [27] where γ is the angular separation between the two spacetime points x and x ′ and ∆t ≡ t − t ′ . The integral kernel G Ψ ℓω depends on the quantum state of the field. For the Boulware [28] (Ψ = B), Unruh [29] (Ψ = U ) and Hartle-Hawking [30] states (Ψ = H) it takes the forms where κ ≡ 1/(4M ) is the surface gravity, andR in/up ℓω denotes the (rescaled) radial factor of the ingoing and upgoing solutions to the wave equation, as defined in Eq. (7). App. 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 Ψ (x; x ′ ) is an analytic expression which is defined in a local neighborhood 1 of the source point x and explicitly shows its singularity structure. Explicitly, it is (e.g., Refs. [31][32][33]) where x ′ ) and w = w(x, x ′ ) are regular and real-valued biscalars. The so-called Synge's world function σ = σ(x, x ′ ) is equal to one-half of the square of the geodesic distance joining x and x ′ , which implies that σ is negative/zero/positive whenever that geodesic is, respectively, timelike/null/spacelike. The biscalars u and v are uniquely determined by the geometry of the spacetime whereas w in principle depends on the quantum state |Ψ⟩. 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: where PV denotes the principal value distribution. This readily yields the anticommutator and the commutator The anticommutator and commutator (which, when multiplied by 'i θ(∆t)' yields the classical retarded Green function) are, respectively dependent and independent of the quantum state |Ψ⟩ of the field.
Eq. (14) shows explicitly the singularity of the Wightman function along σ = 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 Ψ (x; 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 Ψ (x; 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 γ = 0 of the point of emission of the null geodesics as well as along the antipodal line γ = π). Specifically, the pattern for the leading divergence in the real part of the Wightman function is 2 and that in the imaginary part of the Wightman function is (20) where we have omitted the coefficients of the singularity factors. 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, PV(1/σ) and −δ(σ) (corresponding to the direct part in Eq. (14)), whereas after the wave front has crossed one caustic point these turn into, respectively, −δ(σ) and −PV(1/σ). 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. Ref. [37] showed that, for the imaginary part of the Wightman function, its subleading fourfold structure is: (21) Ref. [37] further conjectured that the subleading structure for the real part of the Wightman function is: 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). 2 An exception to the validity of Eqs. (19) and (20) is at caustic points; also, by "σ" in these equations we mean a well-defined extension of the world function outside normal neighborhoodssee [37] about both of these points. 3 Ref. [36] showed what the global singularity structure of the Feynman Green function is, from which that of the Wightman function readily follows.

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. The Unruh-DeWitt model is covariant and causal 4 [39][40][41][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⟩ d (ground state) and |e⟩ d (excited state) which are separated by an energy gap Ω d . The detectors couple to the field amplitudeφ(x d ) along their worldline x d (t) through the interaction Hamil- where λ d is a coupling constant which is dimensionless in (3+1)-dimensional spacetime, 0 ≤ η(t) ≤ 1 is a realvalued switching function, τ d is the detector's proper time and µ d (t) = e iΩdτd(t) |e⟩⟨g| d + e −iΩdτd(t) |g⟩⟨e| d is the monopole operator. Note that the Hamiltonian of (23) generates time translation with respect to coordinate time t, but not the detector's proper time τ 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 and we choose the relation between coordinate time and detector proper time as τ d (t) = f (r d )t. As switching functions for the detector we use Gaussians which as a function of coordinate time read 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 be a product state between the ground states of the two detectors and a field state ρ Ψ : Assuming that the field state has vanishing one-point function, the expansion of the detectors' state after time evolution t = 0...T in coordinate time is Using the basis order |g⟩ a |g⟩ b , |e⟩ a |g⟩ b , |g⟩ a |e⟩ b , |e⟩ a |e⟩ b the final state is represented by the density matrix whose entries are where M Ψ takes complex values while L Ψ dd ′ ≥ 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 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 ρ ab,T we use its negativity. Its perturbative expansion is 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 Ψ and the local noise terms L Ψ aa and L Ψ bb . In particular, if the noise terms are equal, L Ψ dd := L Ψ aa = L Ψ bb , as will be the case in Sec. IV, then the negativity and, thus, the entanglement between the detectors vanishes if L Ψ dd ≥ M Ψ , i.e., when the noise overcomes the correlations. Notice that, because L Ψ aa and L Ψ bb 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, then the entanglement they acquire comes from "extracting" preexisiting 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 Ψ in the leading-order perturbative correction to ρ 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 Ψ into two contributions as M Ψ = M Ψ + + iM Ψ − , where M Ψ + is obtained by replacing W Ψ by its real part (which is symmetric) on the right-hand side of Eq. (29), and M Ψ − by replacing W Ψ with its imaginary part (which is antisymmetric). Note that, in general, M Ψ + and M Ψ − are complex-valued. Below we will encounter generic scenarios where the detectors become entangled while M Ψ is dominated by M Ψ − , and M Ψ + is (almost or exactly) vanishing. In such a situation most 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 N Ψ, (2) , would remain unchanged if we replaced the initial field state by a state which resulted in the same values for L Ψ aa and L Ψ bb (the value of M Ψ − ≈ M Ψ , 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 Ψ + 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 Ψ = M Ψ + 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 Ψ + 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 M Ψ ≈ M Ψ + 6 A quantitative study of the relative contributions of M Ψ + and M Ψ − to the entanglement acquired by the detectors for the case of flat spacetime can be found in Ref. [46]. 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 ∆ 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. 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 Ψ term is dominated by M Ψ − rather than M Ψ + , 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 Ψ + could dominate over M Ψ − because the real part now carries the δ(σ)-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. Fig. 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 Ω aa = Ω bb = 5M −1 and λ = λ a = λ b , 8 are placed on static worldlines at the same radial coordinate r = 6.009M . The angular separation γ 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 ∆ 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 ∆ 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 ∆ ba . In particular, the different insets show the formation of caustics located at angular separation γ = π from the point of origin (which is plotted in red) where the null wave front intersects itself. The caustic reaches r = 6.009M at ∆ ba ≈ 20.7386M .
The four different curves in each panel show the resulting contributions to the final density matrix of the two detectors (28) 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 ∆ 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 B d and the detectors end up entangled. However, in the second panel the detectors' final state remains separable, even for null separated detectors, because ∆ 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. 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 ∆ ba ), the detectors can become entangled at larger switching delays ∆ ba again, at angular separations close to γ = π 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 δ(σ)-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 ± | have a qualitatively similar overall shape in the first and the 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 ∆ ba = 21.5M ) of Fig. 1. In this and the following panels of Fig. 1, the switching delay ∆ ba is larger than the shortest light propagation (Schwarzschild coordinate) time between detectors with angular separation γ = π, 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 γ = π.
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 ∆ ba , while the four different panels correspond to four different angular separations. (Fig. 6 provides a logarithmic plot of the same data.) Note that these plots   Fig. 1, Fig. 7 and Fig. 8, respectively.
over ∆ 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 γ = π/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 γ = π/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π ≲ γ ≲ 0.81π, not even the primary peak overcomes the noise and the detectors' final state remains separable for all switching delays ∆ ba . An example of this is seen in the third (bottom left) panel of Fig. 2 for angular separation γ = 3π/4. However, as the angular separation approaches γ = π 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 γ → π, the primary and secondary peaks superpose and finally create one joint peak aligned at the caustic for γ = π. 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 Tab. 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.
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 γ = π, 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 ∞ / f (r) = T ∞ / 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 ∞ . Analogous to the plots above, Fig. 4 shows the resulting contributions to the detectors' final state for T ∞ = 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 ∆ ba = t 0b − t 0a is chosen appropriately.
To understand the structure of the plots and its dependence on ∆ ba it is useful to consider the light propagation coordinate time ∆t between the two detectors. This is the coordinate time that it takes for a null geodesic starting at spatial coordinates (r, θ = π/2, ϕ) to reach the antipodal point (r, θ = π/2, ϕ + π). Fig. 5 shows this time for the range of radial coordinates we consider. It shows that ∆t is minimal at r = 3M . For this critical case, it takes the value ∆t min = 3 √ 3πM . Setting the switching delay equal to the light propagation time, i.e., setting ∆ ba = ∆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 ∆ ba , one peak forms around r = 3M since the detectors located at this radial coordinate are the first to be exactly null aligned. As ∆ ba is increased in the following panels, a doublepeak structure forms since for ∆ ba > 3 √ 3π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 γ = π. 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 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. We showed that the realistic 3+1-dimensional case possesses a particularly rich phenomenology due to the pres- ence 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. To this end, consider For the t ′ -integral, use Eqs. (A3) and (A7) from Ref. [44], which is For stationary detectors we use τ d (t) = 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 ′ in Eq. (30) takes the following form (denoting r = r d , Analogously, for M Ψ in Eq. (29) we obtain using Eq. (A7), For practical evaluation purposes, one can manipulate the integrals in Eqs. (A9) and (A8) further so that the integration over ω only runs over 0 < ω < ∞. (b) Relative difference between the integral in Eq. (B5) for the Boulware state when integrated up to ωcut and up to ωcut + δ, with δ = 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. 9a, we see that exponentially fast convergence happens after M ωcut becomes slightly larger than the M ω where the asymptotic regime begins. (c) Summand in Eq. (B5) as a function of ℓ for the Boulware state. One can notice from this figure that the summand decays superexponentially with increasing ℓ. In fact, for ℓ > 30, it has already converged to all 10 significant digits that we got when evaluating the integrals. (d) Relative difference between the ℓ-sum in Eq. (B5) for the Boulware state when summed up to lcut and lcut − 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.
convergence of the series becomes slower with increasinḡ r to conclude that the largest ε happens at the highest r we intended to evaluate, which is aroundr max ≈ 5. Then, fixingr =r max we test a couple ofω and ℓ values in the region of interest and conclude that n max = 5000 is enough to satisfy Eq. (B4) in all regions of interest.
For the I ℓω , ρ in/up ℓω coefficients we use data from Ref. [36], which is evaluated using a Mathematica im-plementation of the Mano-Suzuki-Takasugi method. An extensive review of that method can be found in Ref. [50].

L Ψ dd -terms
By setting d = d ′ in Eq. (A8) and manipulating the integration range from −∞ < ω < ∞ to 0 < ω < ∞, we obtain the following expression for the L Ψ dd -terms: (d) Relative difference in the Boulware state integral from (A9), after rewriting the integral so that it ranges from ω = 0 to ω = ωcut for several different ℓ = 0, 10, 20, 30. To evaluate the relative difference, we sum up to ℓcut and then to ℓcut + 1, ranging form ℓcut = 1 up to ℓcut = 100. For ℓcut > 45 all significant digits exactly cancel out when evaluating the relative difference.
Putting everything together we conclude that which leads to the conclusion that, whenever T a ̸ = T b , the whole term decays superexponentially as e − 1 8 (Tb−Ta) 2 ω 2 . On the other hand, if T a = T b ≡ T , which is the case we study in this work, (B13) simplifies to e − (µ 2 a +ν 2 b )T 2 4 , ω → ∞.
Then both, the real and imaginary parts of this quantity decay linearly with ω, since the absolute value of the denominator grows linearly with the frequency. Given that, by considering the extra ω −1 in the integrand in Eq. (A9), we conclude that it decays as ω −2 . Hence, the integrals defining the M Ψ -terms are convergent.
We remark that such leading asymptotic behavior in Eq. (B14) is, in general, not achieved by M ω = 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. 10a. We can further confirm that convergence by looking at the integral in Eq. (A9) as a function of ω cut : as can be seen in Fig. 10b, up to ℓ = 20 we have convergence up to four significant digits, while for larger values of ℓ, the convergence becomes worse, resulting in no significant digits at all. Yet, since the resulting integral decays superexponentially with ℓ, as presented in Fig. 10c, one can still obtain up to four significant digits, as can be seen in Fig. 10d.