Harvesting correlations in Schwarzschild and collapsing shell spacetimes

We study the harvesting of correlations by two Unruh-DeWitt static detectors from the vacuum state of a massless scalar field in a background Vaidya spacetime consisting of a collapsing null shell that forms a Schwarzschild black hole (hereafter Vaidya spacetime for brevity), and we compare the results with those associated with the three preferred vacua (Boulware, Unruh, Hartle-Hawking-Israel vacua) of the eternal Schwarzschild black hole spacetime. To do this we make use of the explicit Wightman functions for a massless scalar field available in (1+1)-dimensional models of the collapsing spacetime and Schwarzschild spacetimes, and the detectors couple to the proper time derivative of the field. First we find that, with respect to the harvesting protocol, the Unruh vacuum agrees very well with the Vaidya vacuum near the horizon even for finite-time interactions. Second, all four vacua have different capacities for creating correlations between the detectors, with the Vaidya vacuum interpolating between the Unruh vacuum near the horizon and the Boulware vacuum far from the horizon. Third, we show that the black hole horizon inhibits any correlations, not just entanglement. Finally, we show that the efficiency of the harvesting protocol depend strongly on the signalling ability of the detectors, which is highly non-trivial in presence of curvature. We provide an asymptotic analysis of the Vaidya vacuum to clarify the relationship between the Boulware/Unruh interpolation and the near/far from horizon and early/late-time limits. We demonstrate a straightforward implementation of numerical contour integration to perform all the calculations.


Introduction
In recent years fruitful progress in our understanding of quantum field theory (QFT) and fundamental physics has been made by applying insights from quantum information theory. In particular, a great deal of attention has been focused on entanglement in quantum field theory. For example, different regions of QFT vacua contain both classical correlations and entanglement even if the regions are causally disconnected [1,2]. Entanglement plays a crucial role in many fundamental phenomena such as the black hole information problem [3][4][5][6]. Viewed from a quantum information perspective, entanglement is a resource that can be used to perform information-processing tasks [7][8][9].

JHEP08(2020)155
The pioneering work of Valentini [10] and later by Reznik et al. [11,12] showed that one can indeed extract entanglement from the QFT vacuum using a pair of initially uncorrelated quantum systems: in recent years this protocol became known as entanglement harvesting. Entanglement harvesting has been shown to be sensitive to accelerations [13], time dependence of the interaction and the number of spacetime dimensions [14], spacetime curvature [15][16][17][18], spacetime topology [19], boundary conditions [20,21], as well as indefinite causal ordering [22][23][24]. Entanglement harvesting has also been investigated in more practical, experimental settings [25][26][27][28][29][30][31]. In such studies, it is common to employ a form of particle-detector model such as the Unruh-DeWitt (UDW) model [32,33], which is a simplified scalar model derived from realistic light-matter interaction. This model captures the essential physics if no angular momentum is exchanged [34]. Furthermore, the model can be easily modified into more complicated variants ('UDW'-like models), such as non-linear couplings [35], coupling to a vector field ('dipole coupling') [34] or to the conjugate momentum of the scalar field ('derivative coupling') [36].
Virtually all studies of entanglement harvesting -and indeed of detector responseare for spacetimes that are static. Rather little is known beyond this case, with two notable exceptions being a study of detector response for a rotating black hole in (2+1) dimensions [37], and a new study of entanglement harvesting in the presence of a gravitational wave [38]. Very recently, the transition rate of a UDW detector coupled to a test scalar field in a (1 + 1)-dimensional Vaidya spacetime background describing collapse of a null shell into a black hole was investigated [39]. The objective was to see whether the traditional expectation that the vacuum of a dynamical collapsing star spacetime is well-described by the so-called Unruh vacuum associated with a test field on a static Schwarzschild background [32]. The conclusion was that in the causal future of the shell this is indeed the case at appropriate limits, and in particular in the long time ('late time') limit. Furthermore, it was also found that the Unruh vacuum gives an upper bound for radiation flux at future null infinity, and a stationary observer carrying a detector at fixed orbit in the exterior black hole region measures a Planckian spectrum in the late time limit [39].
We consider here for the first time the correlation properties of the quantum vacuum in a spacetime containing gravitational collapse. Specifically, we investigate the harvesting protocol in (1+1)-dimensional Vaidya spacetime consisting of a collapsing null shell that forms a black hole (hereafter simply called Vaidya spacetime for brevity), and compare it with the eternal Schwarzschild black hole spacetime. Our study is motivated from two considerations. First, it has been known for quite some time [15] that a single inertial detector cannot distinguish a Minkowski thermal bath from the conformal vacuum of de Sitter expanding universe, but in some regimes two detectors can do so via their entanglement dynamics. Further studies have indicated that the entanglement harvesting properties of a pair of detectors yields novel behaviour that cannot be discerned from that of a single detector [17][18][19][20][21]. It is therefore natural to ask if the Unruh vacuum is still a good approximation of the vacuum describing gravitational collapse (henceforth called the Vaidya vacuum) when two detectors are involved. Second, the detector transition rate was found to be well-approximated by the Unruh vacuum near the horizon and also 'at late times' [39]. While this calculation demonstrated clearly the Planckian behaviour seen in the Unruh vac-JHEP08(2020)155 uum, the detector is turned on both sharply and for infinitely long times after the black hole has formed, thus the notion of 'late time' also subsumes 'long time' limit. We would like to explore this for both single and double detector scenarios for more realistic, strictly finitetime interactions to see how the timescales of the interaction affect the detector dynamics.
We present four main results regarding the dynamics of two detectors interacting with a quantum massless scalar field in a (1+1) Vaidya spacetime in comparison with corresponding results from the three preferred vacua -Boulware, Unruh, Hartle-Hawking-Israelof a Schwarzschild spacetime. First we find that, with respect to harvesting protocol, the Unruh vacuum agrees very well with the Vaidya vacuum even for finite-time interactions near the horizon, complementing recent results for long-interaction times [36]. Second, all four vacua have different capacities for creating correlations between the detectors: the Vaidya vacuum interpolates between the Unruh vacuum near the horizon and the Boulware (and hence Minkowski vacuum) far from the horizon. Third, by investigating the mutual information of the detectors, we show that the black hole horizon inhibits any correlations, not just entanglement, complementing results found in [17]. Finally, we also show that the efficiency of the harvesting protocol depends strongly on the ability of two detectors to signal, which is highly non-trivial in presence of curvature. Our study also includes an asymptotic analysis of the Vaidya vacuum, which clarifies how its respective approximations by Unruh and Boulware/Minkowski vacua are related to the early/late time and near/far from horizon limits of the detector-field interaction.
We shall employ derivative coupling for the detector-field interaction in order to remove the problematic infrared (IR) divergence associated with massless scalar fields in two-dimensional Schwarzschild spacetimes [40]. In addition to reproducing some aspects of the short-distance Hadamard property of the Wightman distribution in (3 + 1) dimensions, this also allows us to make a fairer comparison with Vaidya spacetime (which has no IR divergence). It turns out that we shall require a rather careful numerical treatment after we obtain the general expression for the joint density matrix elements of the two detectors for three reasons. First, the problem lacks a lot of symmetry exploited in other investigations for performing simplifications using Fourier transforms or stationarity; thus numerical treatment becomes essential (especially for the Vaidya vacuum). Second, the Wightman distribution for the Vaidya spacetime has very a complicated pole structure, making the standard i procedure numerically unstable. 1 Finally, transition rate calculations [36,39] involve one-dimensional integrals that allow for better control of the numerics. 2 We take this opportunity to present a straightforward way of circumventing these issues in appendix A without resorting to bottom-up numerical schemes, which makes use of the smoothness and strong support property of Gaussian switching function and basic complex analysis. This method should be applicable to other studies in quantum field theory and beyond where integrals over Green's functions, propagators, kernels, etc. are involved.

JHEP08(2020)155
Our paper is organized as follows. We first review the geometrical aspects of Klein-Gordon QFT in Schwarzschild spacetime in section 2, and then in Vaidya spacetime in section 3. In section 4 we outline the setup for Unruh-DeWitt detector model coupled to (1 + 1)-dimensional background spacetime. In section 5 we present our main results.
In this paper we employ the natural units c = = 1 throughout. We take the metric g to be such that g(V, We also write x ≡ x µ to denote the spacetime coordinates.

Klein-Gordon field in Schwarzschild spacetime
In the next two sections we first review the geometrical and quantum field-theoretic aspects of a quantum massless scalar field on Schwarzschild and Vaidya background spacetimes, 3 as done in [39].

Schwarzschild spacetime: geometry
We will start from the maximal extension of Schwarzschild spacetime, also known as Kruskal-Szekeres extension (M K , g K ), where M K = R 2 ×S 2 . In terms of Kruskal-Szekeres coordinates (U, V, θ, φ), the metric reads , and φ ∈ [0, 2π). For convenience we define M := GM, where M > 0 is the ADM mass of the black hole, and r > 0 can be written in terms of U, V , i.e.
where W(z) is the Lambert W-function (W(z)e W(z) = z) [41]. This spacetime is static, spherically symmetric, asymptotically flat and globally hyperbolic. The point r = 0 is a curvature singularity. The Schwarzschild spacetime (M K , g K ) admits four Killing vector fields. Three of these vector fields, which we denote by ζ 1 , ζ 2 , ζ 3 , are globally spacelike and generate spherical symmetry: while the fourth Killing field, denoted ξ, is given by (2.4) 3 We believe [39] is one of the most concise yet clearest and most illuminating exposition about QFT in a black hole background suitable for our purposes; thus we find it appropriate that we review their descriptions and notation as closely as we can. This Killing field ξ is timelike for r > 2M , spacelike for 0 < r < 2M and null at the hypersurface r = 2M . The null hypersurface r = 2M thus defines a bifurcate Killing horizon. This bifurcate Killing horizon separates M K into four regions, conventionally labelled Region I, II, III, and IV as shown in figure 1.

JHEP08(2020)155
Regions I and II are defined by V > 0 and part of the Killing horizon H + that separates the two regions. In this region, one can use ingoing Eddington-Finkelstein coordinates (U, v, θ, φ) defined by V = e v/(4M ) , with v ∈ R. In this coordinate system, we can view regions I and II as an asymptotically flat and globally hyperbolic spacetime in itself, denoted (M E , g E ), where M E is a submanifold of M K and g E is the induced metric obtained from inclusion map i : M E → M K by pullback, i.e. g E = i * g K . In this coordinate system, the metric reads The Killing vectors for M E are ζ 1 , ζ 2 , ζ 3 , and also ξ under the restriction V > 0, which can now be written as ξ = − 1 4M U ∂ U + ∂ v . Finally, Region I of Schwarzschild spacetime describing the exterior of a static spherically symmetric star or eternal black hole is defined by U < 0 and V > 0. Thus in addition to the ingoing Eddington-Finkelstein coordinates v, one can also introduce the outgoing Eddington-Finkelstein coordinates u defined by U = −e −u/(4M ) . Now using the so-called tortoise radial coordinate r * we can construct two null coordinates We can thus regard Region I as a standalone asymptotically flat and globally hyperbolic spacetime, denoted (M S , g S ), where M S is a submanifold of M K and g S is the induced metric obtained from the inclusion map i : M S → M K by pullback, i.e. g S = i * g K . In this JHEP08(2020)155 coordinate system, the metric reads The Killing vectors for M S are ζ 1 , ζ 2 , ζ 3 , and also ξ under the restriction U < 0 and V > 0, which can now be written as ξ = ∂ t .

Schwarzschild spacetime: Klein-Gordon field
A real massless Klein-Gordon field φ : M → R conformally coupled to gravity in (n + 1)dimensional spacetime M satisfies the Klein-Gordon equation where R is the Ricci scalar curvature. In order to construct an appropriate vacuum state of the theory, we will solve for the classical mode solutions. The general solution takes the form The mode (eigen)functions {u k (x)} satisfy the orthogonality conditions where (f, g) is the Klein-Gordon inner product of f, g given by with respect to the Cauchy surface Σ. The definition of a vacuum state of the field depends on the choice of timelike Killing vector field: given a timelike Killing vector ξ, the mode function u k is said to be positive frequency with respect to ξ if u k (x) solves the eigenvalue equation where L ξ is the Lie derivative with respect to ξ and In our problem, there are three distinguished vacuum states that are invariant under the Killing vector ξ: (a) Boulware vacuum |0 B : this state is defined in Region I and has modes that are positive and negative frequency with respect to Schwarzschild timelike Killing field ∂ t (restriction of ξ to Region I). It is considered unphysical as it is not regular on both future and past horizons H ± . However, this state will be useful for the discussion of the vacuum in the Vaidya background later. (c) Hartle-Hawking-Israel (HHI) vacuum |0 H : this state is defined on the full Kruskal-Szekeres extension and has modes that are positive frequency with respect to both past and future horizon generators ∂ U and ∂ V . This is a state representing a black hole in thermal equilibrium with a radiation bath, such that the restriction of the state to Region I is KMS at the Hawking temperature T H = (8πM ) −1 . Note that T H is the temperature measured by an observer at infinity (see section 5.5).
We restrict our attention to the dimensionally reduced (1+1)-dimensional Schwarzschild spacetime by removing the angular part of the metric in (3+1) dimensions. This allows us to obtain the closed-form expression of the vacuum states in terms of the Wightman two-point distributions by invoking conformal invariance of the Klein-Gordon equation in (1+1) dimensions. The positive frequency modes associated with each vacuum read [42] Hartle-Hawking-Israel : e −iωU , e −iωV , We are interested in constructing the Wightman two-point distribution for each vacuum state (denoted |0 α ), defined by where Λ > 0 is an IR cutoff inherent in (1+1) massless scalar field theory.

Comment on IR ambiguity and derivative coupling
We pause here to comment briefly on the IR divergence in two-dimensional massless field theory that will motivate the choice of detector-field interaction in this paper.

JHEP08(2020)155
It is well-known that a two-dimensional massless scalar field in Minkowski space exhibits an infrared (IR) ambiguity. More specifically, from eq. (2.8) one can show that in (1+1)-dimensional Minkowski space, a massless scalar field quantized in Minkowski coordinates (t, x) corresponding to inertial laboratory frame has Fourier mode decomposition given byφ This decomposition allows us to define Minkowski vacuum |0 M : the Wightman distribution associated to Minkowski vacuum W M (x, x ) can then be shown to have a logarithmic divergence: where > 0 is a ultraviolet (UV) regulator and Λ > 0 is an infrared (IR) regulator. This IR divergence can also be seen from the Fourier mode decomposition, where the integral in eq. (2.16) is divergent for k = 0. We can choose the principal branch of the logarithm so that The second term is formally divergent 5 as Λ → 0. This IR divergence will also appear for the Schwarzschild spacetime as all two-dimensional spacetimes are conformally flat and the Klein-Gordon equation (2.8) is conformally invariant for n = 1; hence the same IR divergence appears in eqs. (2.15a)-(2.15c).
A priori, this IR divergence is problematic for detector dynamics as the density matrix for the detector(s) would depend on the IR cut-off chosen (see e.g. [14,36,43]). Typically, one either chooses Λ based on a characteristic length scale of the system under consideration, or removes it by hand via other arguments. For instance the additive constant that appears in eq. (2.18) can be dropped using the argument that entanglement measures such as concurrence and negativity are by definition infrared-safe [21]: they involve subtraction of two matrix elements that contain the same IR-divergent additive constant. Therefore the formally infinite additive constant drops out of the entanglement calculation. This is analogous to how entanglement entropy in QFT in general contains state-dependent divergences, but a quantity such as relative entropy is finite due to cancellation of divergences (see e.g. [44]).
In this paper, we shall follow [36] instead and consider in section 4 a particular model of detector-field interaction known as derivative coupling. 6 It is a modification of the UDW model where the pullback of the vacuum Wightman distribution along the detector's trajectory is given by the proper time derivative of the field,

JHEP08(2020)155
The proper time derivatives will remove the IR ambiguity from the (1+1)-dimensional Wightman function and the detectors' joint density matrix, and it leads to the same shortdistance behaviour of the Wightman distribution in (3+1) dimensions [36,39]. Furthermore, the derivative Wightman distribution for each vacuum state is indeed invariant under time translation generated by the timelike Killing field ξ, in contrast to the usual Wightman function where the Unruh vacuum is strictly speaking not invariant due to the IR cut-off [39]. Last but not least, analysis of entanglement harvesting for a derivative-coupled UDW-like model where no IR ambiguity occurs has been shown to give similar qualitative results in flat space and (1+1)-dimensional spacetimes with moving mirror [21,34].

Klein-Gordon field in Vaidya spacetime
The geometry of Vaidya spacetime is given by the Lorentzian manifold (M V , g V ) with topology M V = R 2 × S 2 , with the metric written in terms of the ingoing Eddington-Finkelstein-type coordinates (v, r, θ, φ): The Vaidya metric allows for a general class of mass function M (v), and a particularly simple model for null collapse is prescribed by the mass function where Θ(v) is the Heaviside step function and M ≥ 0 is a mass parameter corresponding to ADM mass of the black hole when it is formed by the null shell. The spacetime is isometric to Minkowski space for v < 0 and isometric to Schwarzschild spacetime for v > 0. The conformal diagram is shown in figure 2. For convenience, we shall simply call this particular class of Vaidya metric with mass function (3.2) as Vaidya spacetime. Similar to the Schwarzschild case, the (1+1)-dimensional model for null collapse is obtained by removing the angular coordinates in eq. (3.1). This will enable us to find the vacuum Wightman distribution with respect to the vacuum state of the theory, which we will call the Vaidya vacuum state |0 V , analytically. Solving for the Klein-Gordon equation subjected to Dirichlet boundary conditions at r = 0, the Wightman function for Vaidya spacetime is given by [39] W where u is related to the Kruskal (dimensionless) null coordinate U by with W(z) the Lambert W-function. The function u(U ) can be obtained by matching modes along the null shockwave v = 0 [39], making use in particular the expression for r(U, V ) in eq. (2.2) at the junction.

Unruh-DeWitt model and entanglement harvesting
In this section we review the basics of the Unruh-DeWitt detector model for the description of entanglement harvesting protocol. Although not the original Unruh-DeWitt model (hence sometimes said to be 'UDW-like'), we will call the derivative-coupling version an Unruh-DeWitt model as well for convenience.

Derivative coupling Unruh-DeWitt model
Let M be a (n + 1)-dimensional spacetime manifold and consider two observers Alice and Bob each carrying a pointlike Unruh-DeWitt detector along their respective timelike trajectories 7 x j : R → M , where j ∈ {A, B}. We consider both detectors to be non-inertial, static detectors at fixed radii R j > 2M outside the black hole horizon with R A ≤ R B . Consequently, the detectors experience different gravitational redshifts at their respective locations. Each detector is a two-level quantum system with local interaction HamiltonianĤ j of the detector-field system given by derivative-coupling Hamiltonian where t is a time coordinate for the spacetime. The superscript t is to make clear that the Hamiltonian generates time translation with respect to t. The local interaction between JHEP08(2020)155 each detector and the fieldĤ t j (t) is simpler when it is written as Hamiltonian that generate time translation with respect to the proper time τ , i.e.
Here u µ is the 4-velocity of the detector parametrized by proper time τ , and the two Hamiltonians that generate time translations with respect to t and τ are related by time- λ j denotes the coupling strength andμ j is the monopole moment, given in terms of the detector proper time:μ with Ω j the detector gap,σ ± j the ladder operators of su(2) algebra, and χ j (τ ) is the switching function that controls the duration of interaction. For simplicity we will consider both detectors to be identical, i.e. λ j = λ, Ω j = Ω, with the same Gaussian switching functions where σ prescribes the duration of interaction and τ 0 defines the origin of the proper time. Due to gravitational redshift, the time evolution is best described using a common time provided by the background coordinate system, i.e. usingĤ t I (t). Since we are interested in detector dynamics in Region I, we can use the Schwarzschild time t as a common time coordinate. The time evolution operator is then given bŷ where we have used eq. (4.3) and τ A and τ B are proper times parametrizing different timelike trajectories x A and x B respectively. 8 We also fix the proper times of each detector τ A , τ B such that τ A = τ B = 0 when the Schwarzschild time t = 0, which is possible because the spacetime admits a Cauchy surface given by constant-t slices. We note that due to pointlike nature of the detectors, the derivative coupling particle detector model adopted here is fully covariant [45,46]. For weak coupling, we can perform a Dyson series expansion whose first two terms areÛ

JHEP08(2020)155
and whereÛ (k) is of order λ k . Note that the second order correctionÛ (2) is time-ordered with respect to coordinate time t. Our interest is in vacuum entanglement harvesting, so the initial state is taken to be the uncorrelated state where |g j is the ground state of detector j that satisfiesσ + j |g j = |e j ,σ − j |e j = |g j ; |e j is the excited state of detector j, and |0 α is a vacuum state of the field described in section 2 and 3. The time evolved density matrix is given by ρ =Û ρ 0Û † , and using the Dyson series expansion (4.7) we obtain where ρ (k) is of order λ k : The choice of initial state in eq. (4.9) implies that ρ (1) = 0 since the one-point function 0|φ(x)|0 = 0 for all x. Therefore, the leading order contribution in perturbation theory is ρ (2) . In order to compute the entanglement between the two detectors, we find the joint reduced density matrix of the detectors by tracing out the field's degrees of freedom: Using the ordered basis {|g A |g B , |g A |e B , |e A |g B , |e A |e B }, the matrix representation of ρ AB to leading order reads where the matrix elements are given by where A α is given by (2.19). The local 'noise' terms L ii correspond to the transition probability of detector j, so sometimes we will write this as Pr j (Ω, σ) := L jj . The nonlocal term M depends on the trajectories of both detectors. In the expression for M, we JHEP08(2020)155 have defined In particular we have γ ba = γ −1 ab . For convenience we choose the convention that r B ≥ r A (detector B is at larger radial coordinate than detector A). The constant γ AB in the upper limit of M appears because the time-ordering inÛ (2) needs to account for the redshift factor τ j (t) = f (r j )t. More explicitly, if t = t(τ A ) and t = t (τ B ), it follows from the time-ordering ofÛ (2) ρ 0 (and also ρ 0Û (2) † ) that hence the upper limit in the expression for M in eq. (4.15). Finally, in order to measure the amount of entanglement between the two qubit detectors, there are several faithful entanglement measures we can use. For simplicity, we will use concurrence C[ρ AB ] [47]. For the time-evolved density matrix in our scenario, this has the form [17,19] C to leading order in the coupling. One could also consider entanglement negativity [48] (see e.g. [14,49] for its use in the harvesting setup), but we choose concurrence because it cleanly separates the effect of the non-local term M and local noise L ii on bipartite entanglement. Furthermore, in the special case when the bipartite qubit state is pure, which is the case in this paper, concurrence is a faithful and monotone entanglement measure [47]. We will also be interested in another type of correlation called mutual information, which quantifies the total amount of classical and quantum correlations between the two detectors. The mutual information between the two detectors is defined by where S[ρ] = − Tr ρ log ρ is the von Neumann entropy, and ρ i := Tr j ρ ij is the reduced state for detector i after tracing out the detector j's internal degree of freedom. For the joint density matrix in our scenario [49] to leading order in perturbation theory, where This quantity is useful because if there is no entanglement between the two detectors, then we know that any correlation between them must be either classical correlation or non-entanglement quantum correlations known as quantum discord [50,51]. JHEP08(2020)155

The two-point Wightman distributions for derivative coupling
The remaining task is to calculate the derivative-coupling Wightman function for each of the four vacua |0 α , α = B, U, H, V . Let us use the shorthand . Taking a proper-time derivative of eqs. (2.15a)-(2.15c) and eq. (3.3), we obtain for Schwarzschild vacua Note that for simplicity we have written A U and A H in terms of U instead of U . For the Vaidya vacuum the two-point function has two additional terms (4.23) due to the boundary condition at r = 0.

Comments on switching time and computation of joint density matrix
Now we have all the ingredients to study the correlations between two detectors after interacting with the field. We pause here to make several comments on the procedure of computing the time-evolved density matrix ρ AB to leading order in perturbation theory.
First, note that in our construction the collapsing null shell occurs at v = 0 (this could be generalized to arbitrary v = v 0 but we do not do this here). In terms of the Eddington-Finkelstein coordinates, this means that t + r * = 0. Due to the matching condition at v = 0, it is imperative that for detectors in Region I, the switching time τ = τ 0 is chosen such that it respects v = 0. In particular, if Alice's detector is located at r = kr H for k > 1 and r H = 2M , then inverting the null coordinate v we get the constraint Accounting for redshift, this constraint can be written in terms of detector's proper time: Therefore, if we demand that the Gaussian strong support to be bσ (b > 0), the requirement that this support is contained entirely in Region I imposes the constraint that

JHEP08(2020)155
In this paper we consider 5σ (analogous to "five-sigma standard deviation" in particle physics) to be appropriate and useful for ensuring (4.24), so we set b = 5, though this standard is mathematically somewhat arbitrary. 9 Second, we know that the three standard Schwarzschild vacua have time-translation invariant Wightman functions with respect to the Killing time ξ. Therefore, the excitation probability Pr j (Ω, σ) is invariant under a constant shift of the switching time τ 0 . However, this is not the case for the non-local terms, as the two detectors at two different radii experience different gravitational redshift. Therefore the pullback of the Wightman functions to each detector's trajectory W(x A (τ ), x B (τ )) will not be stationary, i.e. it is not a function of τ − τ .
Third, to our knowledge most UDW model literature to date involves sufficiently simple settings in which numerical integration can be performed relatively straightforwardly, and in some nice cases closed-form expressions can be obtained (see e.g. remarkable calculations in [14,34] for harvesting scenario, or [36,39] for transition rate calculations). In these cases, often the symmetry of the problem allow exact expressions, and in the case of Unruh effect calculations, transition rate is simpler because it is a one-dimensional integral obtained using stationarity of the Wightman distributions. In other contexts such as [17,18,20], the nice properties of AdS 3 spacetime allow analytic computation of both the Wightman functions and reduction of numerical integrals to one-dimensional integrals. The most formidable calculations of the two-point functions of this kind are done e.g. in [53,54], though the objectives are different.
Here we are working with (1) a derivative coupling Wightman distribution, and also (2) a time-dependent collapsing spacetime, which renders the density matrix elements intractable analytically. Therefore a numerical approach is required to make progress. However, it is not hard to check by direct computation that the usual i prescription easily leads to numerical instabilities, and for Vaidya spacetime where the Wightman function has a very complicated pole structure, this is practically impossible without very careful and deliberate control of the integration schemes around the poles. In certain cases, such as the flat space Minkowski vacuum, it may be possible to deal with this by a suitable rewriting of the response function (see e.g. [36,55,56]) in such a way that the i prescription is completely eliminated. However this is an exception to the rule; for example, a spacetime with a static mirror at the origin cannot be dealt with this way as the mirror introduces new poles [21].
In view of the above difficulties, we will compute the joint detector density matrix elements ρ AB using numerical contour integration. Formally, this is equivalent to the i prescription but instead of 'shifting the poles' and taking → 0 (which is numerically unstable in general), we will perform numerical integration that involves a contour in the complex plane. By making a suitable choice of contour that takes into account the exponential suppression of the Gaussian switching functions, we will be able to simplify the JHEP08(2020)155 numerical integration considerably so that no complicated scheme is required. Furthermore, this also serves as a simple demonstration of how contour integration can be useful in a multi-dimensional integral settings that is relatively straightforward to implement as compared to bottom-up numerical schemes. 10 We describe this procedure in appendix A.

Results
In this section we will calculate the amount of correlations that can be extracted by the two qubit detectors and compare the differences between the three preferred states for Schwarzschild background, namely Boulware vacuum |0 B , Unruh vacuum |0 U and Hartle-Hawking-Israel (HHI) state |0 H . We will then compare this to the case where the two detectors are in the black hole exterior Region I of Vaidya spacetime, corresponding to detectors interacting after the black hole collapse has occurred. We will consider both concurrence and mutual correlations as measures of classical and quantum correlations between the two detectors. We close by commenting on the KMS property and detailed balance condition associated with these four vacua.
For numerical computations, we need to choose the nearest distance to the horizon for illustrating the physics very close to the horizon. Let d ij := d(r i , r j ) be the proper distance between two radial coordinates r i , r j . We will impose the condition that the nearest Alice's detector can get to the horizon is given by the proper distance where r H = 2M is the Schwarzschild radius and σ is the switching timescale. We can therefore effectively think of d A = 0.1σ as having Alice's detector to be just above the horizon. We will measure distances in units of σ except in section 5.5 where we need to vary σ. Since by convention we take Bob's detector to be farther from the horizon than Alice's detector, we always have r B ≥ r A . In principle, we could go nearer to, say, d A 0.01σ (since we cannot numerically evaluate the density matrix at r = r H ); however this would take much more optimization and computational time to work with whilst not providing new insights. Our choice is simply a matter of (practical) convenience and simplicity that still includes the relevant physics.
As a side note, we remark that for derivative coupling UDW model, the coupling constant λ has units of [Length] n−1 2 while for amplitude coupling λ has unit [Length] n−3 2 . Consequently, our results will be in terms of dimensionless coupling constantλ := λσ 1−n 2 , where n denotes the number of spatial dimensions. It happens that for derivative coupling in (1+1) dimensions, we haveλ = λ and we will writeλ throughout to remind ourselves that in general coupling constant of UDW model has dimension-dependent units.

Harvesting entanglement
In figure 3 we show the concurrence as a function of proper distance of Alice's detector from the horizon d A . Both Alice and Bob are static, non-inertial observers at fixed Schwarzschild JHEP08(2020)155  radii r A , r B respectively, separated by a fixed proper distance d(r A , r B ) = 2σ in 3(a,b) and d(r A , r B ) = 3σ in figure 3(c,d). We first compare how the four vacua can entangle the two qubits after finite-time interaction.
First, observe that from figure 3(a) that there is an inhibition of entanglement extraction close to the horizon for all states, a result conjectured to hold in general [17] based on a study of this scenario for (2+1) BTZ black holes. Our results support the claim that this is a generic feature of a black hole background, since our choice of detector-field coupling and the choice of states are vastly different, and our example includes the Vaidya vacuum, which is not time-translation invariant in both the state and the black hole background. We also note that the region where concurrence is zero is slightly smaller for the Boulware vacuum and slightly larger for HHI vacuum, which is an indication that all black hole vacua do not have equal 'entangling power' (to use the phrase in [15]). Furthermore, we see that the Unruh vacuum approximates the Vaidya vacuum very well near the horizon even for finite interactions: this result therefore extends the utility of the Unruh vacuum in modelling the vacuum state for collapsing spacetime. For completeness, we note as well that with larger proper separation between the two detectors, entanglement harvesting is JHEP08(2020)155  diminished and since each vacuum entangles differently, it is possible for some vacua to not be able entangle at some distance but other vacua could, as shown in figure 3(d). This result is the black hole equivalent of that found both for accelerating detectors and for comoving detectors in an expanding universe [13,15].
Secondly, for finite time interactions the Unruh vacuum no longer approximates well the Vaidya vacuum as the detectors move far away from the horizon. In figure 3(b), we see that all four vacua distinguish themselves and all have different entangling power, and in particular we note that Vaidya vacuum is an interpolation of Unruh and Boulware vacuum: that is, the Vaidya vacuum is well-approximated by Boulware vacuum as measured by faraway observers while it is well-approximated by Unruh vacuum near the horizon. In figure 4 we separate the local noise contribution due to detector A's excitation and non-local contribution M. In this particular example, the entangling power of the Vaidya vacuum is larger than the Unruh vacuum further from the horizon because the Boulware vacuum has a larger non-local term and smaller local noise. Again we observe that both local noise and non-local terms associated with the Vaidya vacuum interpolate between the Unruh and Boulware vacua; this suggests that the excellent approximation of the Vaidya vacuum by JHEP08(2020)155 = λ is dimensionless coupling constant and we set Ωσ = 2, with the two detectors separated by proper distance d(r A , r B ) = 2σ. The detectors are turned on at τ 0 = 12σ. (a) the entanglement death zone is saturated for large masses, at d A ∼ 4σ for M/σ ≥ 10, and at large masses the concurrence tends to the zerotemperature Boulware limit. (b) shrinking of the difference in entangling power between the four vacua as mass increases. either Unruh or Boulware vacuum is generic and not unique to the entanglement dynamics of the two detectors. We will see that this is true also for mutual information.
In figure 5 we depict how concurrence varies with black hole mass. From figure (5)(a), we can make two observations. First, we see that as black hole mass increases, the entanglement death zone increases until at some point it is saturated for large enough mass. In our example, all vacua for M/σ ≥ 10 have the same death zone, given by d A ∼ 4σ. Second, for large masses, the differences between the different vacua shrinks very quickly: in our example, for M/σ ≥ 2 the four vacua (marked by different line style 11 ) are practically indistinguishable from one another in the plot. Figure 5 shows how a small increase in mass (in units of σ) already shrinks the difference considerably, hence increasing the mass reduces the difference in entangling power of the four vacua. Note that for large masses, the curves for the four vacua overlap and approach the Boulware limit at large distances. This behaviour has a natural interpretation: it can be understood from the fact that as mass increases, the Hawking temperature decreases and hence in the limit of very large mass, the concurrence approaches that of zero-temperature vacuum in the sense of KMS condition [57], i.e. the Boulware vacuum (see section 5.5 for brief discussion on thermality). Such small-mass distinctions are also present for the BTZ black hole [17].

Harvesting mutual information
So far we have looked at entanglement dynamics of the two detectors in a black hole background for various state. We now consider another correlation measure: mutual information. In figures 6 and 7 we plot mutual information as a function of distance from the horizon for the same parameter choices as figure 3. There are four observations that can be made here. First, from the right-hand diagrams in figure 6, again we see that the Vaidya vacuum interpolates between the Unruh and Boulware vacua insofar as mutual information is concerned. Second, we observe that far from the horizon, how much mutual information can be extracted from the Vaidya vacuum relative to the Unruh vacuum depends on detector separation. Comparing the right-hand diagrams in figure 6, larger detector separation tends to make the Vaidya vacuum better at 'imparting' mutual information to the detectors than smaller separation, unlike concurrence where the Vaidya vacuum persistently 'outperforms' the Unruh vacuum at large distances.
Third, the dependence of mutual correlation harvesting on black hole mass is similar to the concurrence, as we show in figure 7 (compare this with figure 5). We again see that the differences in how much mutual information can be extracted from the four vacua shrinks quickly with increasing mass, and that for large masses the mutual information harvesting approaches the Boulware limit. Mutual information also seems to exhibit some saturation near the horizon: in our example, M/σ ≥ 10 seems to share similar mutual information for d A ∼ 7.5σ, and starts to show differences farther away. = λ is dimensionless coupling constant. We set Ωσ = 2 with detector separation d(r A , r B ) = 2σ, and the detectors are turned on at τ 0 = 12σ so that the Gaussian switching peak is very far from the shell. Note that for smaller mass black hole, the 'bifurcation' between Unruh and Vaidya vacuum occurs at larger radius.
The third and the most important observation of this section concerns the behaviour of mutual information dynamics near the horizon. From figures 6(a,c) we see that, independently of the mass, mutual information vanishes quickly as detectors approach the horizon. Furthermore, near the horizon (for every d A > 0) the amount of mutual information that can be harvested can be nonzero, even with zero entanglement (cf. figure 3). Thus very close to the horizon the detectors contain mostly non-entangling correlations. However, the mutual information quickly vanishes as the detectors approach the horizon regardless of detector separation. Figures 6(a,c) show that even though mutual information can be extracted outside the horizon, the black hole horizon effectively extinguishes all correlations -classical or quantum, at least for static detectors.
This result suggests an operational interpretation for detectors in static trajectory (constant r) that a black hole horizon is a null surface that breaks correlations: that is, two static detectors are increasingly unable to harvest any correlations from the vacuum as they get closer to the horizon, and this is independent of which vacuum one uses for the description of the field's ground state. Our results therefore strengthen what was found in [17], in that not only does the horizon inhibit entanglement harvesting, in fact it also suppresses harvesting of any correlations, regardless of the choice of black hole vacuum state. This is most likely related to the fact that detectors cannot maintain static trajectory at r = 2M , thus it is very interesting to see how different the outcome would be when the detectors approach the horizon in free-falling motion; we leave this for future investigations.

Communication between detectors: how timelike are the correlations?
Another natural question to ask in the harvesting protocol is how much of the extracted correlations come from mutual signalling between the two detectors. So far in the literature, JHEP08(2020)155 little attention has been paid to communication between the two detectors when it comes to harvesting correlations in curved spacetimes. The ability of the detectors to communicate important because two uncorrelated quantum systems can be correlated or entangled via signalling or their mutual interactions. 12 In the UDW model, although two detectors interact locally with a common quantum field, due to relativistic causality the two detectors' ability to signal depends on their spacetime separations. In flat space, the notion of spacelike and timelike separation is straightforward and can be given in terms of the coordinate separation Y µ := x µ A − x µ B . In other words, if Y µ Y µ ≤ 0 then the two points are causally connected (see e.g. [14]). When non-compact switching or smearing is involved (such as via Gaussian functions), then one can define two detectors to be spacelike separated whenever the strong support of the switching or smearing functions of one detector is within another detector's causal complement. Relativistic causality of the underlying quantum field in flat space then demands that if the two points are spacelike separated (Y µ Y µ > 0), then the any local observables constructed out of the field operators vanish: Consequently, communication between detectors can be measured in terms of 'signalling estimators' constructed out of the field commutators [58]. It is worth noting that the microcausality condition (5.2) is quite simple to compute even when the background spacetime is curved if the Fourier mode decomposition of the underlying quantum field is known. This is interesting because in curved spacetimes it is generically very difficult in practice to characterize spacelike separation even classically: in presence of curvature, one has to show that the two points cannot be connected by any causal curve [59]. This is especially prohibitive in practice for our detectors because we have to check how the entire Gaussian strong supports of the detectors are contained in each other's causal complement. Here we have a situation where quantum theory simplifies our task of quantifying communication between two detectors at two causally disjoint spacetime regions.
Let us construct a signalling estimator inspired by the construction in [58]: we define where α = B, U, H, V label different vacua and we take the imaginary part since E is purely imaginary. The factor 1/2 is arbitrarily chosen so that in figure 8 the magnitude of E is comparable to the concurrence and aids visualization. Note that we have used the proper time derivative ∂ τφ (x(τ )) instead of the field operatorφ(x(τ )) because the commutator can be easily computed from the derivative Wightman function: as a distribution, this is given by In order to generate entanglement, one would need nonlocal operations in general [7]. Since the field commutator is state-independent we can drop the label α, and also we assumed that the two detectors are identical. Therefore, the estimator can be simplified into Certainly one could construct other estimators using operators associated with the field, but due to eq. (5.2) all observables constructed out of the field operators will behave similarly for spacelike separated regions. The estimators will only give different behaviour for different field observables when the detectors are timelike separated. 13 We superimpose the concurrence for various detector separations with the signalling estimator E as shown in figure 8(a)-(c), using Unruh vacuum as a reference state for concurrence. We can make two important observations here. First, when the detectors become increasingly spacelike, the signalling estimator E is strongly confined to where the concurrence is nonzero. Therefore, at large distances the estimator vanishes for spacelike separation. This is especially manifest in figure 8 where d AB = 8σ would correspond to JHEP08(2020)155 spacelike separated detectors if the spacetime were flat. In figure 8(d) we plot the estimators together for different spacelike separation and we see that the dominant part gets more concentrated to smaller spatial regions.
The second observation is that the estimator suggests that as we bring the two detectors very close to the black hole, the detectors very quickly become spacelike. The estimators fall off very quickly near the horizon. This is very remarkable because it shows that curvature modifies communication in highly non-trivial way: the signalling estimator E is effectively zero near the horizon, followed by an intermediate region the signalling is very enhanced (large |E|), before eventually falling off again as one moves towards spatial infinity. 14 This observation demonstrates that in entanglement harvesting protocol, curvature and communication between detectors are very much related: curvature modifies causal relations between the two detectors as we move the them across different distances from the horizon.
Overall, while we do not explore the (extremely vast) parameter space of our setup, the signalling analysis highlights a mechanism through which entanglement harvesting between two detectors at the black hole exterior occurs: the efficiency of the protocol depends strongly on the ability of detectors to signal between them. The very specific issue of quantum communication in (3+1)D Schwarzschild black hole where angular variables of the metric have important role has been very recently investigated using state-of-the-art calculations in [54].

Vaidya vacuum: near/far from horizon and early/late time limits
Our results thus far suggest that the Vaidya vacuum is an interpolation of the Unruh and Boulware vacua as we move from the horizon towards infinity. In order to better understand these, let us study the late-time and large distance limit of the respective Wightman functions. We stress that our notion of 'late time' is not the same as [39]: latetime means the detectors are turned on for finite duration (σ < ∞) but the peak of the switching function is large (τ 0 → ∞).
First, note that by taking the limit r 2M at fixed coordinate time t (or proper time τ ), the pullback of the Wightman function for the Boulware state approaches the Minkowski value: where A M (x A (τ ), x B (τ )) is the derivative coupling Wightman function for Minkowski vacuum (i.e. derivative version of eq. (2.17) and [36]). Second, for the HHI vacuum the Wightman function would approach that of a thermal bath in Minkowski space with temperature 14 We have also checked the estimator when compact switching in [52] with approximately equal area as the Gaussian switching is used instead and the signalling estimator remains very similar. This provides an indication that despite the non-compact property of the Gaussian switching, the calculations done in this paper work as required.

JHEP08(2020)155
T H = (8πM ) −1 (cf. [60]): For the Unruh vacuum, the Wightman function would approach the "average" of Minkowski vacuum and thermal bath: This averaging makes sense because the Wightman function is constructed by removing "half" of the HHI vacuum's radiation -the ingoing flux (see section 5.5 for further discussion). Now, let us try to make sense of the early time and far from horizon limits. For the Vaidya vacuum, we recall from eq. (3.4) that the Wightman function involves variable u = −4M (1 + W(−U/e)), where W(z) is the Lambert-W function. For fixed coordinate time t (or proper time τ ) and large radial coordinate r, which corresponds to large −U , the asymptotic behaviour of the principal branch of the Lambert-W function is [41] and hence u ∼ t − r = u, where u is a null coordinate in Minkowski space. Therefore, we conclude that in the large −U limit the Vaidya vacuum is well-approximated by the Boulware vacuum (and hence also by Minkowski vacuum of flat space). This happens when either (1) detectors are very far from the horizon (r is very large), or (2) τ 0 is very small (hence t(τ ) along the strong support is small), i.e. detectors turned on very early but still within Schwarzschid exterior Region I shown in figure 2. In particular, this calculation shows that for fixed switching peak τ 0 , once the detectors are sufficiently far away, the detectors cannot tell whether a black hole will form or not because they have the same joint density matrix as if the vacuum were Minkowski, even if they lie within Region I of figure 2.
Let us now make sense of the late-time and near-horizon limit. When the detectors are switched on very late (very large τ 0 ), the behaviour increasingly approaches the Unruh limit. To see this, note that for any large but fixed radius r, one can always make U very small by taking t (or τ ) very large. This happens when we make the Gaussian switching peak τ 0 very large 15 (hence t(τ ) along the strong support is large). In this case, one looks for the other branch of the Lambert-W function and the asymptotic behavior for small −U is [41] W(−U/e) ∼ − log e U = −1 + log U ≈ −1 + U , (5.10) hence u ∼ −4M U = U . This is precisely the null coordinate used for the definition of the Unruh vacuum (see eq. (2.15b) and eq. (4.22b)) Therefore, we conclude that the Vaidya JHEP08(2020)155 = λ is dimensionless coupling constant. We set Ωσ = 2, M σ = 1 2 and d(r A , r B ) = 3σ. The 'bifurcation point' at which the two states begin to show differences in concurrence is now at about d A ≈ 3σ when τ 0 = 5.5σ, as compared to d A ≈ 8σ when τ 0 = 12σ, indicated by the vertical dashed lines.
vacuum is well-approximated by the Unruh vacuum when −U is very small, i.e. either (1) when the detectors are very close to the horizon, or (2) when the detectors are switched on at very late times (even for finite, short interaction timescale σ).
The early/late time limit affecting the Boulware/Unruh approximation of Vaidya vacuum can be visualized in figure 9. We find that the primary factor governing the point at which the approximation breaks down is when the detector is switched on, i.e. the switching peak τ 0 . The earlier the switching time is, the Unruh/Vaidya difference becomes manifest nearer to the horizon. In figure 9, we see that this 'bifurcation' point now begins when detector A is at proper distance of d A ≈ 3σ away from the horizon, 16 as compared to d A ≈ 8σ when τ 0 = 12σ. In other words, for finite-time interaction, how far away from the horizon the Unruh vacuum well-approximates the Vaidya vacuum depends on how early/late the detectors are turned on relative to the null collapse time. This is precisely what we obtained earlier from the asymptotic analysis of the Wightman functions.

Detailed balance condition and non-equilibrium states
In this section we briefly comment on the thermality (or lack thereof) of the black hole vacua we have considered, using the notion of detailed balance.
The Kubo-Martin-Schwinger (KMS) condition prescribes the necessary condition for a quantum state to be thermal with respect to some timelike Killing vector ξ with KMS temperature T = β −1 [61,62]. A particularly operational formulation is given by the detailed balance condition [57], which says that in the adiabatic limit (long, carefully switched on JHEP08(2020)155 interaction), we have lim σ→∞ F (Ω, σ) lim σ→∞ F (−Ω, σ) = e −β loc Ω , (5.11) where T loc = β −1 loc is the local temperature as measured an observer in curved spacetime. In our scenario, the local temperature will be given by the Tolman temperature [63,64], where T H = (8πM ) −1 is the Hawking temperature of the Schwarzschild black hole. The function F (Ω, σ) is usually called the response function (or transition rate [36,40]), defined as the excitation probability Pr(Ω, σ) per unit time: Recall that the expression for the excitation probability of a single detector is given precisely by the matrix element L jj , which takes the form (dropping the index j since we only consider one detector, cf. eq. (4.14)) where σ is the interaction timescale encoded in χ(t). Note that by symmetry, P (−Ω, σ) corresponds to the de-excitation probability of a detector from its excited state to its ground state. Formally, the detailed balance condition (5.11) is formulated in terms of the adiabatic limit of the response function F . We are interested in the physical limit where the interaction timescale σ is slowly increased to mimic longer and longer interaction time, and our choice of Gaussian switching guarantees that the interaction is smooth enough. Therefore, the ratio of the responses can instead be computed in terms of the excitation-to-deexcitation (EDR) ratio Recently the EDR ratio has been used to study interesting new phenomena, such as the Unruh effect without thermality [65], anti-Unruh effects [66,67], and a novel anti-Hawking effect [68].
In figure 10(a) we show how the EDR ratio R increases with interaction timescale σ near the horizon, where we consider a static detector at r = 2.2M . As expected, the HHI vacuum approaches the detailed balance prediction (horizontal dotted line) since it corresponds to black hole in thermal equilibrium with a thermal bath. Likewise, the Boulware vacuum approaches the expected value of zero since it corresponds to a zero-temperature state with β → ∞ and hence e −β loc Ω → 0. By construction, we know that Unruh and Vaidya vacua correspond to non-equilibrium states since they are supposed to simulate radiation flux to infinity; here we again see that the Unruh and Vaidya vacua agree very JHEP08(2020)155 well in terms of their EDR near the horizon. This provides further evidence that at the level of single detector responses, the Unruh vacuum approximates very well the vacuum of a collapsing star spacetime near the horizon. Inspection of figure 10(a) reveals that in the late time limit, the response functions of both the Unruh and Vaidya vacua approach a constant value that is below the detailed balance prediction (5.15) reflecting the non-thermal nature of these states. We can say more about the long time limit of EDR for these two states. From the Wightman functions (cf. section 2), one can see that the Unruh vacuum can be thought of as an average of Boulware and HHI states since the positive frequency modes are defined by u = v − 2r * and U . It then follows that in the long time limit, the EDR for the Unruh vacuum reads Observe that the long-time transition rate for Unruh vacuum is precisely the average of the response of the Boulware and HHI vacua, and for Ω > 0 Unruh vacuum has Planckian spectrum but with the density of states halved from the HHI vacuum. Substituting into eq. (5.16), we obtain the long time limit of the EDR ratio for the Unruh vacuum: Indeed, R U (Ω) is precisely the asymptotic value of the EDR ratio that both the Unruh and Vaidya vacua approach in figure 10(a). What this tells us is that although the Unruh vacuum is not an equilibrium state, one can obtain the effective local temperature associated with only the outgoing flux of this state (the right-moving modes) by computing In figure 10(b), however, we see that the EDR ratio for the Vaidya vacuum approaches that of the Boulware vacuum once it is far from the horizon (in this example we take r = 21.1r H ) for the same choice of switching peak τ 0 as figure 10(a). This can be anticipated from earlier results in the entanglement harvesting calculation, where we saw that the transition probability approaches the Boulware vacuum far from the horizon. This is interesting because a faraway observer in Region I of Vaidya spacetime computing the EDR ratio would conclude that the vacuum state is KMS with zero temperature, associated with a Boulware vacuum. Therefore, a single detector faraway interacting for finite times will not be able to infer the thermal behaviour of the radiation outflux to infinity for the Vaidya vacuum 17 once it is far enough from the black hole (or equivalently, switched on early enough; see the asymptotic analysis in section 5.4).
Finally, we note that for finite time interactions, the EDR ratio for the Vaidya vacuum approaches that of the Unruh vacuum as one increases τ 0 . More concretely, if at fixed radius and fixed σ, we pick τ 0 such that the EDR ratio agrees with the Boulware vacuum (as we did in figure 10(b)) increasing τ 0 will bring the EDR curve towards the Unruh one. This complements our earlier asymptotic analysis in section 5.4. From this perspective, the Unruh vacuum can be properly thought of as a physical, 'late-time' limit of the Vaidya vacuum (corresponding to τ 0 → ∞ limit), and this limit is approached very quickly. This

JHEP08(2020)155
is in contrast to the notion of the late time limit in [36], which has more to do with when the detector is switched off after it is turned on, thus more closely related to how long the interaction is switched on.

Conclusion
We have studied the phenomenon of harvesting of correlations by two detectors from the vacuum state of a massless scalar field in background Vaidya spacetime, and compare the results with those associated to the three preferred vacua (Boulware, Unruh, Hartle-Hawking-Israel vacua) in Schwarzschild spacetime. We use the derivative coupling particle detector model where the Wightman functions have similar short-distance behaviour as the Wightman functions in the (3+1)-dimensional counterpart, as well as to resolve the infrared ambiguities associated to massless scalar fields in (1 + 1)-dimensional spacetimes. We perform these studies using a straightforward implementation of numerical contour integration, outlined in appendix A.
Let us summarize our results. First, we showed that from the perspective of harvesting of correlations between two detectors, near the horizon the Unruh vacuum agrees very well with the Vaidya vacuum even for finite-time interactions, complementing the longinteraction result from [36]. Second, all four vacua have different capacities for creating correlations between the detectors, with the Vaidya vacuum's capacity interpolating between that of the Unruh vacuum near the horizon and the Boulware vacuum far from the horizon. Third, our examination of mutual information indicated that the black hole horizon inhibits any correlations, not just entanglement, complementing the results found in [17]. Finally, efficiency of the harvesting protocol depends strongly on the signalling ability of the two detectors, which is highly non-trivial in the presence of curvature. We have also studied the asymptotic behaviour of the Vaidya vacuum analytically to understand how it approximates the Boulware/Minkowski vacuum in the early time/large distance limit, and approximates the Unruh vacuum in the late time/near-horizon limit. Our asymptotic analysis clarifies the distinction between the late-time and long-time limits in transition rate calculations [36,39].
A natural extension to our results would be to analyse the correlations between two detectors, one of which free-falling through the horizon during the interactions. Since our results on the exterior region shows that the horizon inhibits all forms of correlations from being extracted from the vacuum, how free-falling detectors break correlations between them is an operational question related to information problem in black hole thermodynamics. Another related question concerns the effect of null shockwaves from the perspective of supertranslation [69][70][71]: it would be interesting to study the correlations between two detectors in presence of supertranslations. Finally, it would also be interesting to see how the four vacua considered in this paper affect communication efficiencies such as channel capacities between detectors [72]. We leave these questions for future work. where Ω > 0. We can write this as

JHEP08(2020)155
For this particular case, we are fortunate because the closed-form expression for eq. (A.5) is known, which we can use to check our calculations. This is given by [19] We remark that there is another closed form expression derived differently in [75] that only works correctly for Ω > 0, whereas eq. (A.6) is valid for all Ω ∈ R.
Let us now compare this with numerical computation. 20 We will refer to an integral J j as Method j , and we will call J 0 Method 0.

A.1 Method 1: direct i integration
We denote J 1 to be the integral (A.5) evaluated by brute force, picking a small enough during integration. We will evaluate this for Ωσ = 1 for concreteness. In this case Method 0 gives The values of J 1 can be computed using various settings and optimizations. In our case, reasonable results are obtained using MinRecursion → 3, MaxRecursion → 20, AccuracyGoal → ∞, PrecisionGoal → 10. We also cut the integral at strong support, i.e. t, t ∈ (−5σ, 5σ) for better convergence; one can check that the results are generally worse if one chooses to numerically integrate over R. The results are shown in table 1.
Observe that ∼ 10 −2 σ reasonably approximates (A.7), but the rest of the values do not work. To our knowledge, any other settings within this scheme do not help much, and we believe that while in principle there should be a way to make this method work, it would require a great deal of effort and understanding of the back-end numerical analysis to make this worthwhile in terms of both the computational time and numerical stability. We stress JHEP08(2020)155 that the sorts of computations done in [36] or [20] have one particular advantage: they can be recast into one-dimensional integrals that can be dealt with much better numerically. For example, the Method → "DoubleExponentialOscillatory" used in [17] is not available for higher-dimensional integrals.

A.2 Method 2: numerical contour integration
The idea is basically to perform the following integral: where C( ) is a contour deformed to the upper complex plane around the pole t = t. The contour is shown in figure 11 and shown contrasted to the i prescription. If we choose to instead perform the integral over t first, then the contour is deformed to the lower complex plane around the pole t = t . We let here to be the distance from the pole (in units of σ): that is, we integrate t from −∞ to t − , then from t − to t − + i, then to t + + i, followed by t + and finally from t → ∞. That is, we set to be the distance from the pole along the t axis. 21 Again we integrate over strong support (−5σ, 5σ) as in general integration over R is of lower quality. The results are shown in table 2. Notice that the results are a much better approximation to (A.7) than Method 1. Furthermore, to achieve the quality shown above, the only setting needed was MinRecursion → 3 and nothing more. This could be improved with more optimization. It is quite remarkable that this method works very well with minimal settings whereas the usual i approach of Method 1 fails terribly. Method 2 only starts to deviate very little when we get too close to the pole ( ∼ 10 −5 σ) due to numerical resolution.
We make four observations here. First, the fact that J 2 is numerically constant across a broad range of values of is a manifestation of a basic principle in complex analysis, namely the deformation theorem. The theorem states that within a holomorphic region we can deform the contour of an integral without changing the value of the integral, which JHEP08(2020)155  follows from Cauchy's integral theorem. Since the pole is along the t axis, any will give the same result since there is no other pole in the upper complex plane. Therefore Method 2 provides a very nice way of checking numerical stability: if the integral is no longer constant as we vary across a broad but reasonable range, (recall from table 2 that we would not want to be too small numerically), then perhaps one needs to check if something has gone wrong or the method itself no longer works stably. Second, because of the deformation theorem, in practice the contour shown in figure 11 is flexible: we chose this contour because we felt this to be the simplest for illustration. Third, notice that in computing J (Ω, σ, ), Jordan's lemma cannot be used due to the Gaussian switching function -hence we do not have the benefit of using the residue technique numerically. Finally, what we have performed here is effectively a two-dimensional contour integration, where the poles are continuous (one pole on the t axis for every t) on the (t, t ) plane. We pause to remark that actually there are two more methods that work well for Minkowski vacuum calculations, which are used in [36,40,55,56]. One of them in fact can be written in a form free of the UV regulator , thus it is either correct or incorrect. A brief investigation [73] indicated that for Minkowski vacua these two are competitive methods and behave very well. However they failed in the presence of a (possibly dynamical) Dirichlet boundary condition (such as a moving mirror [21]) because the mirror introduces new poles that reduce the utility of these methods. Nevertheless the resultant calculations remained valid because the transition rate is given by a one-dimensional integral in time. It was shown that numerical contour integration remains superior in the context of moving mirrors [73].
The results in this section and the observations above testify to the especial appeal of numerical contour integration in practical calculation.

A.3 Better contour for entanglement harvesting: Vaidya spacetime
For calculations carried out in this paper, the choice of contour in figure 11 is not good enough. The problem is that the contour we picked relies on finding the location of the poles, i.e. we are solving for t = f (t). Even for derivative Wightman functions for the Unruh and HHI vacua in eq. (4.22b) and (4.22c), f is in general not a linear function of t and depends on the black hole mass M . Even more importantly, for entanglement JHEP08(2020)155 harvesting the nonlocal term M depends on two different radial coordinates with different gravitational redshifts, so f is highly non-trivial.
An even bigger problem is caused by time ordering in M: one would have to constantly track whether within the strong support the poles are included or not when time ordering is applied. For derivative Wightman functions in Vaidya spacetime, this is worsened by the two additional equally complicated terms. It is not hard to check that the contour prescription we did earlier, where the deformation is somewhat close to the poles (say ∼ 10 −1 σ) did not quite work, let alone direct integration via i .
The deformation theorem and Gaussian suppression coming from the switching function provide us with a new contour that we can use. 22 The idea is that if the strong support of the Gaussian switching is (−bσ, bσ), then we set = bσ. In other words, we adopt the contour in figure 12. This contour has the advantage that we effectively do away with finding the poles: the poles are either within the within the strip or outside the strip, and a single contour covers all possible positions of the pole independent of the complexity of the function f that describes the location of the pole as a function of t. It also solves the issue of time ordering, by replacing the upper limit of t from 5σ to t.
As it turns out, with relatively minimal settings (such as MinRecursion ∼ 3 − 6), the numerical computation works very well even for Vaidya spacetime. We also chose for simplicity to deform the contour by one imaginary unit +i to the upper complex plane. In fact, this numerical calculation is stable enough for computation of long-time observables such as the EDR ratio for KMS detailed balance condition. We have also used infinite precision such as using fractions 1/2 instead of 0.5 whenever possible, and we set the global precision setting to 50 digits using PreRead command.
We note that with more complicated problems (and depending on the issue at hand), we expect that different optimizations and variations may be needed on top of what we have done here. The main point is that numerical contour integration provides a sufficiently robust and straightforward implementation without having to construct a separate numerical scheme from scratch, and the fact that optimization is possible at all (unlike the direct integration by i prescription). We also did not attempt to optimize computational time; this is perhaps left for future investigations.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.