# Quantum correlation of Hawking radiations for mass-superposed BTZ black
holes

Hongbin Zhang  Baocheng Zhang<EMAIL_ADDRESS>School of Mathematics
and Physics, China University of Geosciences, Wuhan 430074, China

###### Abstract

Recent studies using the Unruh-DeWitt detector model have revealed the quantum
signatures of black hole mass superposition. In this paper, we use the equal-
time density-density correlation function to reinvestigate the mass
superposition of BTZ black holes. Interestingly, the quantum signatures of BTZ
black hole mass superposition can be revealed by the equal-time density-
density correlation peaks. Additionally, different mass superpositions for the
same mass ratio can be distinguished, a phenomenon not previously documented.
Specifically, the correlation function method can pinpoint the location where
Hawking radiations are generated, even for mass-superposed BTZ black holes.
Our findings indicate that this location is approximately at two and a half
times the radius of the horizon, where the correlation function reaches its
maximum. This supports the concept of a quantum atmosphere as the origin locus
of Hawking radiation, even in quantum-superposed spacetimes.

Key words: quantum correlations, black holes, mass superposition, Wightman
function.

###### Contents

  1. I Introduction
  2. II Wightman functions
  3. III Density-density correlation functions
  4. IV Quantum signatures
    1. IV.1 Obtaining quantum signals through density-density correlation functions
    2. IV.2 Distinguishing quantum signals with same mass ratio
  5. V Creation and evolution of Hawking radiations
    1. V.1 Correlation across the horizon
    2. V.2 Time evolution of correlation
  6. VI Conclusions

##  I Introduction

Einstein’s general relativity predicts that a sufficiently compact mass can
deform spacetime to form a black hole wald1999gravitational . Black holes
provide an arena where interesting effects associated with gravity, both
classical and quantum, can manifest. The boundary beyond which nothing can
escape is called the event horizon. Quantum field theory in curved spacetime
birrell1984quantum  predicts that black holes are not truly black, as they
emit thermal radiation known as Hawking radiation hawking1974black . Pairs of
correlated particles are created near the black hole event horizon; one
particle escapes to infinity as Hawking radiation, while the other is
swallowed by the singularity inside the black hole horizon, known as the
partner particle. Hawking’s discovery is considered a milestone of modern
theoretical physics hawking1975particle .

The exact location where Hawking radiation originates is an interesting
question currently being researched. Unruh first discussed the origin of
Hawking radiation from gravitational black holes unruh1977origin . Subsequent
studies analyzed analog black holes schutzhold2008origin  and used two-point
correlation functions to reanalyze gravitational black holes
schutzhold2010quantum . Later, Giddings introduced the concept of a quantum
atmosphere, displaced from the horizon, as the origin of Hawking radiation
giddings2016hawking . Understanding the origin of Hawking quanta is crucial,
as it can address the information loss problem hawking1976breakdown ;
almheiri2021entropy ; zhang2009hidden ; zhang2011entropy . Consequently,
research on the origin of Hawking radiation has been extended to various
scenarios hod2016hawking ; dey2017black ; buss2018quantum ; dey2019black ;
eune2019test ; kaczmarek2024signatures ; ong2020quantum .

Recently, the phenomenon of black hole mass superposition has been studied in
the context of spacetime superpositions zych2019bell ; henderson2020quantum .
Examples include quantum superpositions of Minkowski spacetime foo2023quantum
and the signatures of rotating black holes in quantum superpositions
suryaatmadja2023signatures . Researchers are also attempting to apply these
concepts to analog black holes barcelo2022analogue ; barcelo2022chronology ,
building on earlier work unruh1981experimental ; barcelo2011analogue ;
shi2023quantum ; zhang2023shielding , and testing quantum effects in gravity
through tabletop experiments christodoulou2019possibility ;
marletto2017gravitationally ; bose2017spin . Black hole mass superpositions
provide a novel approach to studying quantum black holes and advancing the
theory of quantum gravity.

Black holes in mass superpositions also exhibit Hawking radiation, measurable
by the Unruh-DeWitt detector crispino2008unruh . This detector reveals
signatures of quantum-gravitational effects when black hole mass
superpositions are formed with specific mass ratios foo2022quantum . The
quantum signatures of black hole mass superpositions support Bekenstein’s
conjecture bekenstein1973black ; bekenstein2020quantum  that a complete theory
of quantum gravity should treat black holes as quantum objects. According to
Bekenstein, the black hole horizon area and mass are adiabatic invariants with
an associated discrete quantization hodPhysRevLett.81.4293 . Bekenstein and
Mukhanov proposed that mass quantization should alter the nature of Hawking
radiation, resulting in a discrete emission spectrum with evenly spaced energy
levels gr-qc/9505012 . This idea has drawn significant attention to the study
of quantum black holes 1 ; 2 ; 3 ; 4 ; 5 ; 6 . The quantization of black hole
mass implies that black holes can also exist in superpositions. The
possibility of black hole mass superpositions has been explored within the
framework of loop quantum gravity 7 ; 8 ; 9 ; 10 ; 11 .

Since mass-superposed black holes can emit Hawking radiation, investigating
the origin of this radiation is crucial. Correlation functions are a valuable
tool for studying the origin of Hawking radiation. They have been used to
analyze the occurrence position and timing of Hawking radiation in analog
black holes balbinot2008nonlocal ; carusotto2008numerical ; fabbri2021ramp ;
liao2019proposal ; michel2016phonon , despite the very low Hawking radiation
temperature in such systems. Density-density correlation has become a key
piece of evidence for measuring Hawking radiation from analog black holes in
laboratory settings steinhauer2016observation ; munoz2019observation . This
method has also been applied to determine the origin of Hawking radiation from
gravitational black holes parentani2010vacuum ; balbinot2022quantum ;
anderson2022horizons ; balbinot2023quantum ; fontana2023stress ;
balbinot2023hawking . In this paper, we will study the density-density
correlation functions for mass-superposed BTZ black holes to gain more insight
into mass superposition and the original location of Hawking radiation.

This paper is organized as follows. In the second section, we introduce the
BTZ black hole metrics and analyze the Wightman functions for both classical
BTZ black holes and mass-superposed BTZ black holes. In the third section,
based on the Wightman functions, we construct the density-density correlation
function for mass-superposed BTZ black holes. Using the density-density
correlation function outside the event horizon of mass-superposed BTZ black
holes, we demonstrate the quantum signatures of black hole mass superpositions
in the fourth section. In the fifth section, we investigate the origin of
Hawking radiation using the density-density correlation across the horizon of
both classical and mass-superposed BTZ black holes. Finally, we present our
conclusions in the sixth section.

##  II Wightman functions

In (2+1)-dimensions spacetime, Bañados, Teitelboim, and Zanelli found a
solution of the Einstein gravitational field equation called BTZ black hole
banados1992black . Geometrically, (2+1)-dimensions AdS space can be obtained
from flat space $\mathbb{R}^{2,2}$. AdS3 space is a submanifold carlip19952+

$\displaystyle-l^{2}=-T_{1}^{2}-T_{2}^{2}+X_{1}^{2}+X_{2}^{2},$ (1)

in flat space $\mathbb{R}^{2,2}$

$\displaystyle\mathrm{d}s^{2}=-\mathrm{d}T_{1}^{2}-\mathrm{d}T_{2}^{2}+\mathrm{d}X_{1}^{2}+\mathrm{d}X_{2}^{2}.$
(2)

The above embedding coordinates $(X_{1},X_{2},T_{1},T_{2})$ can be chosen as
smith2014looking

$\displaystyle\begin{split}T_{1}&=\tilde{r}\cosh\tilde{\phi},\\\
X_{1}&=\tilde{r}\sinh\tilde{\phi},\\\
T_{2}&=\sqrt{\tilde{r}^{2}-l^{2}}\sinh\frac{\tilde{t}}{l},\\\
X_{2}&=\sqrt{\tilde{r}^{2}-l^{2}}\cosh\frac{\tilde{t}}{l},\end{split}$ (3)

and then we get the AdS-Rindler metric

$\displaystyle\mathrm{d}s^{2}$
$\displaystyle=-(\frac{\tilde{r}^{2}}{l^{2}}-1)\mathrm{d}\tilde{t}^{2}+\frac{1}{(\frac{\tilde{r}^{2}}{l^{2}}-1)}\mathrm{d}\tilde{r}^{2}+\tilde{r}^{2}\mathrm{d}\tilde{\phi}^{2}.$
(4)

The BTZ metric can be derived by the new coordinates $r=\tilde{r}\sqrt{M}$,
$t=\tilde{t}/\sqrt{M}$ and $\phi=\tilde{\phi}/\sqrt{M}$ foo2022quantum ,

$\displaystyle\mathrm{d}s^{2}$
$\displaystyle=-(\frac{r^{2}}{l^{2}}-M)\mathrm{d}t^{2}+\frac{1}{(\frac{r^{2}}{l^{2}}-M)}\mathrm{d}r^{2}+r^{2}\mathrm{d}\phi^{2},$
(5)

where $M$ is the mass of BTZ black hole, $l$ is the AdS length scale in an
asymptotically AdS vacuum, and $\Lambda=-{1/l^{2}}$ is the cosmological
constant. In particular, the event horizon is located at $r_{h}=\sqrt{M}l$.
Here, we emphasize that in this paper, the discussed BTZ black holes are
considered without angular momentum for the sake of convenience in discussing
mass-superposed black holes. In these scenarios, the horizon undergoes a
superposition characterized by $r_{s}=\sqrt{M_{1}}l$ and
$r_{b}=\sqrt{M_{2}}l$, defined by the two mass values $M_{1},M_{2}$
($M_{1}<M_{2}$) of the black hole mass superposition.

To obtain the density-density correlation function we have to know the
Wightman function of BTZ black hole spacetime. This Wightman function can be
constructed from the well-known vacuum $AdS_{3}$ Wightman function
lifschytz1994scalar ; carlip19952+

$\displaystyle W_{AdS}(x,x^{\prime})=\frac{1}{4\pi
l\sqrt{2}}\left(\frac{1}{\sqrt{\sigma(x,x^{\prime})}}-\frac{\zeta}{\sqrt{\sigma(x,x^{\prime})+2}}\right),$
(6)

where the parameter $\zeta$ is the boundary condition at infinity in $AdS$
spacetime, taking the values $\\{0,1,-1\\}$ corresponding to transparent,
Dirichlet, and Neumann boundary conditions, respectively. The function
$\sigma(x,x^{\prime})$ is the squared geodesic distance between $x$ and
$x^{\prime}$ in the flat $\mathbb{R}^{2,2}$ space given by

$\displaystyle\sigma(x,x^{\prime})=\frac{1}{2l^{2}}[-\left(T_{1}-T_{1}^{\prime}\right)^{2}-\left(T_{2}-T_{2}^{\prime}\right)^{2}+\left(X_{1}-X_{1}^{\prime}\right)^{2}+\left(X_{2}-X_{2}^{\prime}\right)^{2}].$
(7)

The Wightman functions for the BTZ black hole spacetime can be derived from
AdS vacuum by

$\displaystyle
W_{BTZ}\left(x,x^{\prime}\right)=\frac{1}{\sum_{p}\eta^{2p}}\sum_{p}\sum_{q}\eta^{p}\eta^{q}W_{\mathrm{AdS}}\left(\Gamma^{p}x,\Gamma^{q}x^{\prime}\right),$
(8)

where $\Gamma^{p}$ is the periodic identification of the angle coordinate
$\tilde{\phi}\rightarrow\tilde{\phi}+2\pi p\sqrt{M}$ and the parameter
$\eta=+1,-1$. Inserting the $AdS_{3}$ embedding coordinates
$(X_{1},X_{2},T_{1},T_{2})$ together with the periodic identification into Eq.
$(\ref{eq.7})$, we gain the squared geodesic distance in the BTZ black hole
spacetime as

$\displaystyle\sigma\left({\Gamma}^{p}x,{\Gamma}^{q}x^{\prime}\right)=$
$\displaystyle\frac{\tilde{r_{1}}\tilde{r_{2}}}{l^{2}}\cosh\left[2\pi(p-q)\sqrt{M}\right]-1$
$\displaystyle-\sqrt{\frac{\tilde{r_{1}}^{2}}{l^{2}}-1}\sqrt{\frac{\tilde{r_{2}}^{2}}{l^{2}}-1}\cosh\frac{\tilde{t_{1}}-\tilde{t_{2}}}{l},$
(9)

and the cross term of squared geodesic distance as

$\displaystyle\sigma\left({\Gamma_{1}}^{p}x,{\Gamma_{2}}^{q}x^{\prime}\right)=$
$\displaystyle\frac{\tilde{r_{1}}\tilde{r_{2}}}{l^{2}}\cosh\left[2\pi(p\sqrt{M_{1}}-q\sqrt{M_{2}})\right]-1$
$\displaystyle-\sqrt{\frac{\tilde{r_{1}}^{2}}{l^{2}}-1}\sqrt{\frac{\tilde{r_{2}}^{2}}{l^{2}}-1}\cosh\frac{\tilde{t_{1}}-\tilde{t_{2}}}{l},$
(10)

where ${\Gamma_{1}}^{p}$ and ${\Gamma_{2}}^{q}$ are two different periodic
isometries with two different mass $M_{1}$ and $M_{2}$ foo2022quantum . For
simplicity, we study the Wightman functions with transparent boundary
condition $\zeta=0$. In the end, we can get the Wightman function for the BTZ
black hole spacetime as

$\displaystyle
W_{BTZ}\left(x,x^{\prime}\right)=\frac{1}{\sum_{p}\eta^{2p}}\frac{1}{4\pi
l\sqrt{2}}\sum_{p}\sum_{q}\frac{1}{\sqrt{\sigma\left({\Gamma}^{p}x,{\Gamma}^{q}x^{\prime}\right)}},$
(11)

and the Wightman function of the mass-superposed BTZ black hole spacetime as

$\displaystyle
W_{BTZ}^{(1,2)}\left(x,x^{\prime}\right)=\frac{1}{\sum_{p}\eta^{2p}}\frac{1}{4\pi
l\sqrt{2}}\sum_{p}\sum_{q}\frac{1}{\sqrt{\sigma\left({\Gamma_{1}}^{p}x,{\Gamma_{2}}^{q}x^{\prime}\right)}}.$
(12)

Now we consider a massless scalar field $\hat{\psi}$ in (2+1)-dimension AdS
spacetime. The Wightman function of the BTZ black hole with mass $M$ in Eq.
$(\ref{eq.11})$ can be rewritten as henderson2018harvesting

$\displaystyle
W_{BTZ}\left(x,x^{\prime}\right)=\frac{1}{\mathcal{N}}\sum_{p,q}\eta^{p}\eta^{q}\langle
0|\hat{\psi}\left(\Gamma^{p}x\right)\hat{\psi}\left(\Gamma^{q}x^{\prime}\right)|0\rangle,$
(13)

and the Wightman function of the mass-superposed BTZ black hole with
superposed mass $M_{1}$ and $M_{2}$ in Eq. $(\ref{eq.12})$ can be rewritten as

$\displaystyle
W_{BTZ}^{(1,2)}\left(x,x^{\prime}\right)=\frac{1}{\mathcal{N}}\sum_{p,q}\eta^{p}\eta^{q}\langle
0|\hat{\psi}\left(\Gamma_{1}^{p}x\right)\hat{\psi}\left(\Gamma_{2}^{q}x^{\prime}\right)|0\rangle,$
(14)

where the parameter $\eta=+1,-1$ corresponding to the twisted and untwisted
fields respectively, $\mathcal{N}=\sum_{p}\eta^{2p}$ is normalization factor,
and $|0\rangle$ is the AdS vacuum. The establishment of the connection between
the AdS vacuum state and the Whiteman functions will facilitate our
construction of the density-density correlation function.

##  III Density-density correlation functions

In this section, we will utilize Wightman functions to build the density-
density correlation functions. One objective of constructing these correlation
functions is to gather additional insights into BTZ black holes in mass
superpositions. Another aim is to pinpoint where the signal indicating the
presence of Hawking radiation emerges in both the classical BTZ black hole and
the mass-superposed BTZ black hole.

We choose free-falling observers that are locally inertial observers in the
BTZ black hole spacetime. The free-falling observer can measure the energy
density balbinot2022quantum

$\hat{\rho}=\hat{T}_{ab}u^{a}u^{b},$ (15)

where

$\hat{T}_{ab}(\hat{\phi})=\partial_{a}\hat{\phi}\partial_{b}\hat{\phi}-\frac{g_{ab}}{2}g^{cd}\partial_{c}\hat{\phi}\partial_{d}\hat{\phi}$
(16)

is the energy-momentum operator of the massless scalar field $\hat{\phi}$ that
describes an outgoing flux of particles at the Hawking temperature and
$u^{a}=(1,V)$ is the four-velocity of the free-falling observer. For the
convenience to make the physical description from these observers, we
transform the BTZ metric in Eq. $(\ref{eq.5})$ into $Painlev\acute{e}$
coordinates painleve1921cr

$\displaystyle\mathrm{d}s^{2}$
$\displaystyle=-f\mathrm{d}\tau^{2}-2V\mathrm{d}\tau\mathrm{d}r+\mathrm{d}r^{2}+r^{2}\mathrm{d}\phi^{2},$
(17)

where $f=\frac{r^{2}}{l^{2}}-M$ and $V=-\sqrt{1-(\frac{r^{2}}{l^{2}}-M)}$. So
we gain the vanishing of the trace of $\hat{T}_{ab}$

$g^{ab}\hat{T}_{ab}=-\hat{T}_{\tau\tau}-2V\hat{T}_{\tau r}+f\hat{T}_{rr}=0,$
(18)

and the energy density

$\hat{\rho}=\hat{T}_{ab}u^{a}u^{b}=\hat{T}_{\tau\tau}+2V\hat{T}_{\tau
r}+V^{2}\hat{T}_{rr},$ (19)

when we consider $\phi$ as a constant. Then we combine Eq. $(\ref{eq.18})$
with Eq. $(\ref{eq.19})$ to give the energy density as

$\hat{\rho}=(f+V^{2})\hat{T}_{rr}=\hat{T}_{rr},$ (20)

which is obtained using the $Painlev\acute{e}$ time

$\tau=t+\int\frac{\sqrt{1-(r^{2}/l^{2}-M)}}{r^{2}/l^{2}-M}\mathrm{d}r=t+r^{\ast}.$
(21)

According to massless scalar field in Eq. $(\ref{eq.16})$, the energy density
can be rewritten as

$\hat{\rho}=(1-\frac{f}{2})\partial_{r}\hat{\phi}\partial_{r}\hat{\phi}+\frac{1}{2}\partial_{\tau}\hat{\phi}\partial_{\tau}\hat{\phi}+V\partial_{\tau}\hat{\phi}\partial_{r}\hat{\phi}.$
(22)

Finally, we gain the density-density correlation function from Eq.
$(\ref{eq.22})$ as

$\displaystyle G(x,x^{\prime})=\langle
0|\hat{\rho}(x)\hat{\rho}(x^{\prime})|0\rangle$ $\displaystyle=\langle
0|[(1-\frac{f_{1}}{2})\partial_{r_{1}}\hat{\phi}\partial_{r_{1}}\hat{\phi}+\frac{1}{2}\partial_{\tau_{1}}\hat{\phi}\partial_{\tau_{1}}\hat{\phi}+V_{1}\partial_{\tau_{1}}\hat{\phi}\partial_{r_{1}}\hat{\phi}]$
$\displaystyle\ast[(1-\frac{f_{2}}{2})\partial_{r_{2}}\hat{\phi}\partial_{r_{2}}\hat{\phi}+\frac{1}{2}\partial_{\tau_{2}}\hat{\phi}\partial_{\tau_{2}}\hat{\phi}+V_{2}\partial_{\tau_{2}}\hat{\phi}\partial_{r_{2}}\hat{\phi}]|0\rangle,$
(23)

where $f_{1}=f(r_{1})$, $f_{1}=f(r_{2})$, and

$\hat{\phi}=\frac{1}{\sqrt{\mathcal{N}}}\sum_{p}\eta^{p}\hat{\psi}\left(\Gamma^{p}x\right)$
(24)

is an automorphic field of BTZ black hole spacetime foo2022quantum . In order
to make our paper clear and easy to read, we define four functions that are
related to density-density correlation functions by

$\displaystyle\begin{split}&G_{1}=\partial_{r_{1}}\partial_{r_{2}}\langle
0|\hat{\phi}(x)\hat{\phi}(x^{\prime})|0\rangle=\partial_{r_{1}}\partial_{r_{2}}W_{BTZ},\\\
&G_{2}=\partial_{r_{1}}\partial_{\tau_{2}}\langle
0|\hat{\phi}(x)\hat{\phi}(x^{\prime})|0\rangle=\partial_{r_{1}}\partial_{\tau_{2}}W_{BTZ},\\\
&G_{3}=\partial_{\tau_{1}}\partial_{r_{2}}\langle
0|\hat{\phi}(x)\hat{\phi}(x^{\prime})|0\rangle=\partial_{\tau_{1}}\partial_{r_{2}}W_{BTZ},\\\
&G_{4}=\partial_{\tau_{1}}\partial_{\tau_{2}}\langle
0|\hat{\phi}(x)\hat{\phi}(x^{\prime})|0\rangle=\partial_{\tau_{1}}\partial_{\tau_{2}}W_{BTZ},\end{split}$
(25)

where the Wightman function is just that in Eq. $(\ref{eq.13})$ for the
classical BTZ black hole or that in Eq. $(\ref{eq.14})$ for the mass-
superposed BTZ black hole. Finally, we obtain the form of the density-density
correlation function as

$\displaystyle G(x,x^{\prime})=$
$\displaystyle(1-\frac{f_{1}}{2})(1-\frac{f_{2}}{2}){G_{1}}^{2}+\frac{1}{2}(1-\frac{f_{1}}{2}){G_{2}}^{2}+V_{2}(1-\frac{f_{1}}{2})G_{1}G_{2}+\frac{1}{2}(1-\frac{f_{2}}{2}){G_{3}}^{2}$
$\displaystyle+\frac{1}{4}{G_{4}}^{2}+\frac{V_{2}}{2}G_{3}G_{4}+V_{1}(1-\frac{f_{2}}{2})G_{1}G_{3}+\frac{V_{1}}{2}G_{2}G_{4}+V_{1}V_{2}G_{1}G_{4}.$
(26)

##  IV Quantum signatures

###  IV.1 Obtaining quantum signals through density-density correlation
functions

We will use the correlation functions obtained in the last section to give the
quantum signatures of the black hole mass superposition. Firstly, we express
the squared geodesic distance with the radial coordinate $r=\tilde{r}\sqrt{M}$
and $Painlev\acute{e}$ time $\tau=t+r^{\ast}$ as

$\displaystyle\sigma\left({\Gamma}^{p}x,{\Gamma}^{q}x^{\prime}\right)=\frac{r_{1}r_{2}}{Ml^{2}}\cosh\left[2\pi(p-q)\sqrt{M}\right]-1$
$\displaystyle-\sqrt{\frac{r_{1}^{2}}{Ml^{2}}-1}\sqrt{\frac{r_{2}^{2}}{Ml^{2}}-1}\cosh\frac{\tau_{1}-\tau_{2}+r_{2}^{\ast}-r_{1}^{\ast}}{l},$
(27)

and the cross term of squared geodesic distance as

$\displaystyle\sigma\left({\Gamma_{1}}^{p}x,{\Gamma_{2}}^{q}x^{\prime}\right)=\frac{r_{1}r_{2}}{\sqrt{M_{1}}\sqrt{M_{2}}l^{2}}\cosh\left[2\pi(p\sqrt{M_{1}}-q\sqrt{M_{2}})\right]-1$
$\displaystyle-\sqrt{\frac{r_{1}^{2}}{M_{1}l^{2}}-1}\sqrt{\frac{r_{2}^{2}}{M_{2}l^{2}}-1}\cosh\frac{\tau_{1}-\tau_{2}+r_{2}^{\ast}-r_{1}^{\ast}}{l}.$
(28)

Then, we substitute Eq. $(\ref{eq.27})$ and Eq. $(\ref{eq.28})$ into Eq.
$(\ref{eq.26})$ to gain the correlation functions of the classical BTZ black
hole and the mass-superposed BTZ black hole. In particular, correlations can
be calculated from the black hole vacuum outside the event horizon
anderson2022horizons ; smith2014looking ; gallock2021harvesting . It is noted
that when the superposed mass takes the same value ($M_{1}=M_{2}=M$) in Eq.
$(\ref{eq.28})$, the density-density correlation reduces to the same form as
in Eq. $(\ref{eq.27})$.

Black hole mass superposition induces novel effects that are absent in
classical spacetime. As shown in Ref. foo2022quantum , the detector exhibits
resonances at specific mass ratios when interacting with the Hawking radiation
from mass-superposed BTZ black holes. These resonances can be understood in
light of Bekenstein’s conjecture regarding the quantization of black hole mass
bekenstein1973black . Quantum signatures can also be observed through the
density-density correlation function, as depicted in Fig. 1, which
corroborates findings in Ref. foo2022quantum  for the mass ratio
$M_{1}/M_{2}=1/n^{2}$. Furthermore, it is evident that the value of the equal-
time density-density correlation peak gradually decreases with increasing $n$
in the mass ratio. This indicates that larger differences between the two mass
values in the black hole mass superposition lead to reduced density-density
correlation peaks.

Figure 1:  The absolute value of equal-time $G(x,x^{\prime})$ outside the
horizon of mass-superposed black holes for four different mass pairs:
$(M_{1}=1,M_{2}=9)$, $(M_{1}=1,M_{2}=16)$, $(M_{1}=1,M_{2}=25)$, and
$(M_{1}=1,M_{2}=36)$, corresponding to mass ratios, $M_{1}/M_{2}=1/n^{2}$ with
$n=3,4,5,6$, respectively. Other parameters are taken as $r_{1}=r_{2}=10$,
$l=1$, $p=q=20$ and $\tau_{1}=\tau_{2}$. Figure 2:  The absolute value of
qual-time $G(x,x^{\prime})$ outside the horizon of mass-superposed black holes
for four different mass pairs: $(M_{1}=1,M_{2}=9/4)$, $(M_{1}=1,M_{2}=16/9)$,
$(M_{1}=1,M_{2}=25/16)$, and $(M_{1}=1,M_{2}=36/25)$, corresponding to mass
ratios $M_{1}/M_{2}=(n-1)^{2}/n^{2}$ with $n=3,4,5,6$, respectively. Other
parameters are kept consistent with those in Fig. $(\ref{fig1})$.

Regarding the resonance phenomenon observed at the mass ratio
$M_{1}/M_{2}=(n-1)^{2}/n^{2}$ as reported in Ref. foo2022quantum , we
replicate these findings through calculations of the equal-time density-
density correlation, illustrated in Fig. 2. In contrast to the scenario
depicted in Fig. 1 at the mass ratio $M_{1}/M_{2}=1/n^{2}$, we observe that
the value of the equal-time density-density correlation peak increases
gradually with higher values of $n$ in the mass ratio.

Comparing Fig. 1 with Fig. 2, it becomes apparent that the equal-time density-
density correlation value for the mass ratio $M_{1}/M_{2}=(n-1)^{2}/n^{2}$ is
greater than that for $M_{1}/M_{2}=1/n^{2}$ when other parameters remain
constant. This difference arises because the disparity between the two mass
values in the former mass ratio is smaller than in the latter. Specifically,
for each mass ratio, the correlation peak values for black hole mass
superpositions increase as the two mass values become closer to each other.

The equal-time density-density correlation also exhibits peaks at certain mass
ratios that lie between the special values $1/n^{2}$ or $(n-1)^{2}/n^{2}$, but
these peaks are smaller than those observed at these specific mass ratios.
Furthermore, the values of $p$ and $q$ in Eqs. (27) and (28) affect the
occurrence of these correlation peaks. Higher values of $p$ and $q$ result in
more peaks in the density-density correlation for black hole mass
superpositions. In particular, when the radial coordinates $r_{1}=r_{2}$, the
equal-time density-density correlation peaks diminish as the radial coordinate
moves further away from the horizon $r_{b}=\sqrt{M}l$, defined by the larger
mass of the two black hole mass superpositions.

###  IV.2 Distinguishing quantum signals with same mass ratio

It has been observed that resonant peaks occur for $M_{1}<M_{2}$ or
$M_{1}>M_{2}$ at specific mass ratios, but distinguishing between different
mass pairs with the same mass ratio, such as $(M_{1}=0.5,M_{2}=2)$,
$(M_{1}=1,M_{2}=4)$ corresponding to $M_{1}/M_{2}=1/n^{2}$ with $n=2$, has not
been clear using earlier methods foo2022quantum . However, the density-density
correlation method can achieve this distinction.

In Table 1, we compute and present four absolute values of $G(x,x^{\prime})$
for superposed mass values: $(M_{1}=4,M_{2}=16)$, $(M_{1}=1,M_{2}=4)$,
$(M_{1}=4,M_{2}=1)$, and $(M_{1}=0.5,M_{2}=2)$. From this table, it is evident
that different absolute values are obtained for the density-density
correlation even at the same mass ratios in the first, second, and fourth
rows. Furthermore, these absolute values increase as the two mass values of
the black hole mass superpositions increase.

Table 1: The absolute value of $G(x,x^{\prime})$ outside the horizon of mass-superposed black holes for four different mass pairs: $(M_{1}=0.5,M_{2}=2)$, $(M_{1}=1,M_{2}=4)$, $(M_{1}=4,M_{2}=16)$, and $(M_{1}=4,M_{2}=1)$. Other parameters are kept consistent with those in Fig. $(\ref{fig1})$.  Mass ratios $(M_{1},M_{2})$ | $4,16$ | $1,4$ | $4,1$ | $0.5,2$
---|---|---|---|---
$G(x,x^{\prime})(10^{-6})$ | $6.048$ | $1.163$ | $1.163$ | $0.480$

(a)

(b)

Figure 3:  The absolute value of equal-time $G(x,x^{\prime})$ outside the
horizon of mass-superposed black holes for two different mass pairs:
$(M_{1}=1,M_{2}=4)$ (a), $(M_{1}=4,M_{2}=1)$ (b). Other parameters are kept
consistent with those in Fig. $(\ref{fig1})$.

Moreover, it is noted from the table that the absolute value of
$G(x,x^{\prime})$ for the mass values $(M_{1}=1,M_{2}=4)$ is identical to that
for $(M_{1}=4,M_{2}=1)$. This indicates that these configurations are
indistinguishable, demonstrating an exchange symmetry in the density-density
correlation functions for the mass-superposed parameters $M_{1}$ and $M_{2}$,
as illustrated in Fig. 3. Specifically, it is seen from Eq. (28) that this
symmetry arises when the parameters $p$ and $q$ are set to the same value in
our calculations.

##  V Creation and evolution of Hawking radiations

###  V.1 Correlation across the horizon

The above analysis indicates that the mass ratio affects the density-density
correlation outside the horizon of mass-superposed black holes. In this
section, we will study the origin of Hawking radiation through the density-
density correlation, as that made for the origin of Hawking radiation from
quantum atmosphere parentani2010vacuum ; balbinot2022quantum ;
anderson2022horizons ; balbinot2023quantum ; fontana2023stress .

(a)

(b)

Figure 4:  Three-dimensional plot of equal-time $G(x,x^{\prime})$ across the
horizon of the classical BTZ black hole (a) and the mass-superposed BTZ black
hole (b). The parameters are taken as $M=4$ for the classical BTZ black hole
(a) and $M_{1}=1,M_{2}=4$ for the mass-superposed BTZ black hole (b). Other
parameters are taken as $r_{1}=r_{2}=10$, $l=1$, $p=q=10$ and
$\tau_{1}=\tau_{2}=1$.

To determine the origin of Hawking particles from the quantum atmosphere in
both classical BTZ black holes and mass-superposed BTZ black holes, it is
essential to examine the density-density correlation function across the
horizon balbinot2008nonlocal ; carusotto2008numerical ; fabbri2021ramp ;
liao2019proposal ; michel2016phonon . This requires rescaling the time
coordinate to $\tilde{t}=t\sqrt{M}$ and deriving the new squared geodesic
distance as

$\displaystyle\sigma\left({\Gamma}^{p}x,{\Gamma}^{q}x^{\prime}\right)=\frac{r_{1}r_{2}}{Ml^{2}}\cosh\left[2\pi(p-q)\sqrt{M}\right]$
$\displaystyle-1-\sqrt{\frac{r_{1}^{2}}{Ml^{2}}-1}\sqrt{\frac{r_{2}^{2}}{Ml^{2}}-1}\cosh\sqrt{M}\left(\frac{\tau_{1}-\tau_{2}+r_{2}^{\ast}-r_{1}^{\ast}}{l}\right)$
(29)

and the new cross term of squared geodesic distance as

$\displaystyle\sigma\left({\Gamma_{1}}^{p}x,{\Gamma_{2}}^{q}x^{\prime}\right)=\frac{r_{1}r_{2}}{\sqrt{M_{1}}\sqrt{M_{2}}l^{2}}\cosh\left[2\pi(p\sqrt{M_{1}}-q\sqrt{M_{2}})\right]$
$\displaystyle-1-\sqrt{\frac{r_{1}^{2}}{M_{1}l^{2}}-1}\sqrt{\frac{r_{2}^{2}}{M_{2}l^{2}}-1}\cosh\frac{\sqrt{M_{1}}\tau_{1}-\sqrt{M_{2}}\tau_{2}+\sqrt{M_{2}}r_{2}^{\ast}-\sqrt{M_{1}}r_{1}^{\ast}}{l}.$
(30)

Then, we substitute Eq. $(\ref{eq.29})$ and Eq. $(\ref{eq.30})$ into Eq.
$(\ref{eq.26})$ to get the new density-density correlation function
$G_{1}(x,x^{\prime})$ of classical BTZ black holes and mass-superposed BTZ
black holes.

In Fig. 4, we depict the density-density correlation functions
$G(x,x^{\prime})$ across the event horizon of both the classical BTZ black
hole and the mass-superposed BTZ black hole. It is evident that a distinct
peak appears approximately at a position 2 to 3 times the distance of the
event horizon $r_{h}=\sqrt{M}l=2$ for the classical BTZ black hole, and at
about 2.5 times the distance of the horizon $r_{b}=\sqrt{M}l=2$, defined by
the larger mass of the two mass values in the black hole mass superposition.

(a)

(b)

(c)

(d)

Figure 5:  The absolute value of $G_{1}(x,x^{\prime})$ across the horizon of
mass-superposed BTZ black holes at four increasing times: $\tau_{1}=5.2$ (a),
$\tau_{2}=5.9$ (b), $\tau_{3}=6.4$ (c), and $\tau_{4}=9$ (d), with the
superposed mass set as $(M_{1}=1,M_{2}=4)$. Other parameters are kept
consistent with those in Fig. $(\ref{fig4})$.

###  V.2 Time evolution of correlation

For classical gravitational black holes, including the BTZ black hole, the
equal-time density-density correlation cannot depict time evolution since the
time term vanishes at equal times $\tau_{1}=\tau_{2}$, as shown in Eq. (27).
However, the time term remains present in the density-density correlation for
mass-superposed BTZ black holes, as illustrated in Eq. (30). Therefore, we can
analyze the evolution of Hawking radiation through the correlation function in
mass-superposed BTZ black holes.

In Fig. 5, we compute and present the absolute values of the density-density
correlation at four increasing times $\tau_{1}=5.2$, $\tau_{2}=5.9$,
$\tau_{3}=6.4$, and $\tau_{4}=9$. The peak at $r_{s}=r_{1}=1$ represents the
horizon location defined by the smaller mass of the two mass values in the
black hole mass superposition. For correlations where one point is located
less than $r_{s}$ and the other is located beyond the horizon $r_{b}$ defined
by the larger mass of the two mass values, the density-density correlation
peak initially increases to a maximum and then decreases. Moreover, when one
point is between $r_{s}$ and $r_{b}$, and the other is beyond $r_{b}$, the
correlation peaks increase with time. These changes in peak values manifest
the evolution of particle creation in the spacetime of mass-superposed BTZ
black holes.

Particularly noteworthy is that the density-density correlation peak moves
closer to the horizon $r_{s}=\sqrt{M}l=1$ as time progresses. This observation
indicates that partner particles of Hawking radiation in mass-superposed BTZ
black holes do not fall into the singularity but accumulate near the horizon
$r_{s}$ over time. This contrasts with findings for gravitational
Schwarzschild black holes in Ref. balbinot2022quantum .

##  VI Conclusions

In this paper, we study the equal-time density-density correlation function in
the context of classical BTZ black holes and mass-superposed BTZ black holes.
We find the quantum signatures of BTZ black hole mass superposition at
specific mass ratios, $\sqrt{M_{1}/M_{2}}=1/n$ and
$\sqrt{M_{1}/M_{2}}=(n-1)/n$, by observing the peaks in the equal-time
density-density correlation function outside the event horizon. These peaks
diminish when the mass values in the black hole mass superpositions are
significantly different. However, the peak values increase as the mass values
of the BTZ black hole mass superposition increase. Notably, mass-superposed
black holes with the same mass ratio but different mass values can be
distinguished by the values of the density-density correlation peaks, which
has not been discussed before.

We also investigate the location and evolution of Hawking radiation for mass-
superposed BTZ black holes using the equal-time density-density correlation
function. Since the location of the peaks in the correlation function
represents the creation of Hawking particle pairs, we find that the
characteristic peak appears at approximately 2 to 2.5 times the distance of
the event horizon $r_{h}=\sqrt{M}l=2$ for classical BTZ black holes, and at
about 2.5 times the horizon distance $r_{b}=\sqrt{M}l=2$ for mass-superposed
BTZ black holes. As time progresses, the peaks in the density-density
correlation of the mass-superposed BTZ black holes gradually increase to a
maximum and then decrease. Simultaneously, the peaks move toward the horizon
$r_{s}$ over time. This movement indicates the evolution of the creation of
Hawking particle pairs to some extent.

A clearer description requires further understanding of the superposed horizon
in mass-superposed black holes, but it is possible to extend the present
research to superpositions of gravitationally collapsed BTZ black holes
saini2018gravitational , exotic BTZ black holes townsend2013thermodynamics ,
or the analog gravity simulation of spacetime superpositions
barcelo2022analogue ; barcelo2022chronology .

## Acknowledgments

This work is supported by National Natural Science Foundation of China (NSFC)
with Grant No. 12375057, and the Fundamental Research Funds for the Central
Universities, China University of Geosciences (Wuhan).

## References

  * (1) R.M. Wald,  Gravitational collapse and cosmic censorship,  Black holes, gravitational radiation and the universe: essays in honor of CV Vishveshwara, Springer Netherlands (1999) pages 69–86.
  * (2) N.D. Birrell and P. Davies,  Quantum fields in curved space,  Cambridge University Press, Cambridge, England (1984).
  * (3) S.W. Hawking,  Black hole explosions?, _Nature_ 248 (1974) 30.
  * (4) S.W. Hawking,  Particle creation by black holes, _Commun. Math. Phys._ 43 (1975) 199.
  * (5) W. Unruh,  Origin of the particles in black-hole evaporation, _Phys. Rev. D_ 15 (1977) 365 .
  * (6) R. Schützhold and W. Unruh,  Origin of the particles in black hole evaporation, _Phys. Rev. D_ 78 (2008) 041504 [gr-qc/0804.1686].
  * (7) R. Schützhold and W. Unruh,  Quantum correlations across the black hole horizon, _Phys. Rev. D_ 81 (2010) 124033 [gr-qc/1002.1844].
  * (8) S.B. Giddings,  Hawking radiation, the Stefan–Boltzmann law, and unitarization, _Phys. Lett. B_ 754 (2016) 39 [hep-th/1511.08221].
  * (9) S.W. Hawking,  Breakdown of predictability in gravitational collapse, _Phys. Rev. D_ 14 (1976) 2460.
  * (10) A. Almheiri, T. Hartman, J. Maldacena, E. Shaghoulian, and A. Tajdini,  The entropy of Hawking radiation, _Rev. Mod. Phys._ 93 (2021) 035002 [hep-th/2006.06872].
  * (11) B. Zhang, Q.Y. Cai, L.You, and M.S. Zhan,  Hidden messenger revealed in Hawking radiation: a resolution to the paradox of black hole information loss, _Phys. Lett. B_ 675 (2009) 98 [hep-th/0903.08932].
  * (12) B. Zhang, Q.Y. Cai, M.S. Zhan, and L.You,  Entropy is conserved in Hawking radiation as tunneling: a revisit of the black hole information loss paradox, _Ann. Phys._ 326 (2011) 350 [hep-th/0906.5033].
  * (13) S. Hod,  Hawking radiation and the Stefan–Boltzmann law: The effective radius of the black-hole quantum atmosphere, _Phys. Lett. B_ 757 (2016) 121 [gr-qc/1607.02510].
  * (14) R. Dey, S. Liberati, and D. Pranzetti,  The black hole quantum atmosphere, _Phys. Lett. B_ 774 (2017) 308 [gr-qc/1701.06161].
  * (15) C. Buss and M. Casals,  Quantum correlator outside a schwarzschild black hole, _Phys. Lett. B_ 776 (2018) 168 [hep-th/1709.05990].
  * (16) R. Dey, S. Liberati, Z. Mirzaiyan, and D. Pranzetti,  Black hole quantum atmosphere for freely falling observers, _Phys. Lett. B_ 797 (2019) 134828 [gr-qc/1906.02958].
  * (17) M, Eune and W, Kim,  Test of quantum atmosphere in the dimensionally reduced Schwarzschild black hole, _Phys. Lett. B_ 798 (2019) 135020 [gr-qc/1908.03374].
  * (18) A.Z. Kaczmarek and D. Szcześniak,  Signatures of the black hole quantum atmosphere in nonlocal correlations, _Phys. Lett. B_ 848 (2024) 138364 [gr-qc/2306.09941].
  * (19) Y.C. Ong and M.R.R. Good,  Quantum atmosphere of Reissner-Nordström black holes, _Phys. Rev. Res._ 2 (2020) 033322 [gr-qc/2003.10429].
  * (20) M. Zych, F. Costa, I. Pikovski, and Č. Brukner,  Bell’s theorem for temporal order, _Nat. Commun._ 10 (2019) 3772 [quant-ph/1708.00248].
  * (21) L.J. Henderson, A. Belenchia, E. Castro-Ruiz, C. Budroni, M. Zych, Č. Brukner, and R.B. Mann,  Quantum temporal superposition: the case of quantum field theory, _Phys. Rev. Lett._ 125 (2020) 131602 [quant-ph/2002.06208].
  * (22) J. Foo, C.S. Arabaci, M. Zych, and R.B. Mann,  Quantum superpositions of Minkowski spacetime, _Phys. Rev. D_ 107 (2023) 045014 [gr-qc/2208.12083].
  * (23) C. Suryaatmadja, C.S. Arabaci, M.P.G. Robbins, J. Foo, M. Zych, and R.B. Mann,  Signatures of rotating black holes in quantum superposition, [gr-qc/2310.10864].
  * (24) C. Barceló, L.J. Garay, and G. García-Moreno,  Analogue gravity simulation of superpositions of spacetimes, _Eur. Phys. J. C_ 82 (2022) 1 [gr-qc/2104.15078].
  * (25) C. Barceló, J. E.Sánchez, G. García-Moreno, and G. Jannes,  Chronology protection implementation in analogue gravity, _Eur. Phys. J. C_ 82 (2022) 299 [gr-qc/2201.11072].
  * (26) W. Unruh,  Experimental black-hole evaporation?, _Phys. Rev. Lett._ 46 (1981) 1351.
  * (27) C. Barcelo, S. Liberati, and M. Visser,  Analogue gravity, _Living Rev. Relativity_ 14 (2011) 1 [gr-qc/0505065].
  * (28) Y.H. Shi, R.Q. Yang, Z. Xiang, Z.Y. Ge, H. Li, Y.Y. Wang, K. Huang, Y. Tian, X. Song, D. Zheng et al,  Quantum simulation of hawking radiation and curved spacetime with asuperconducting on-chip black hole, _Nat. Commun._ 14 (2023) 3263 [quant-ph/2111.11092].
  * (29) H. Zhang and B. Zhang,  Shielding of penrose superradiance in optical black holes, _Eur. Phys. J. C_ 83 (2023) 1 [gr-qc/2303.14418].
  * (30) M. Christodoulou and C. Rovelli,  On the possibility of laboratory evidence for quantum superposition of geometries, _Phys. Lett. B_ 792 (2019) 64 [gr-qc/1808.05842].
  * (31) C. Marletto and V. Vedral,  Gravitationally induced entanglement between two massive particles is sufficient evidence of quantum effects in gravity, _Phys. Rev. Lett._ 119 (2017) 240402 [quant-ph/1707.06036].
  * (32) S. Bose, A. Mazumdar, G.W. Morley, H. Ulbricht, M. Toroš, M. Paternostro, A.A. Geraci, P.F. Barker, M.S.Kim, and G. Milburn,  Spin entanglement witness for quantum gravity, _Phys. Rev. Lett._ 119 (2017) 240401 [quant-ph/1707.06050].
  * (33) L.C.B. Crispino, A. Higuchi, and G.E.A. Matsas,  The Unruh effect and its applications, _Rev. Mod. Phys._ 80 (2008) 787 [gr-qc/0710.5373].
  * (34) J. Foo, C.S. Arabaci, M. Zych, and R.B. Mann,  Quantum signatures of black hole mass superpositions, _Phys. Rev. Lett._ 129 (2022) 181301 [gr-qc/2111.13315].
  * (35) J.D. Bekenstein,  Black holes and entropy, _Phys. Rev. D_ 7 (1973) 2333.
  * (36) J.D. Bekenstein,  The quantum mass spectrum of the kerr black hole,  in JACOB BEKENSTEIN: The Conservative Revolutionary World Scientific, Singapore, (2020), pp.331–334.
  * (37) S. Hod,  Bohr’s Correspondence Principle and the Area Spectrum of Quantum Black Holes, _Phys. Rev. Lett._ 81 (1998) 4293 [gr-qc/9812002].
  * (38) J. D. Bekenstein and V. Mukhanov,  Spectroscopy of the quantum black hole, _Phys. Lett. B_ 360 (1995) 7 [gr-qc/9505012].
  * (39) J. D. Bekenstein,  Quantum black holes as atoms, in Proceedings of the 8th Marcel Grossmann Meeting on Recent Developments in Theoretical and Experimental General Relativity, Gravitation and Relativistic Field Theories (MG 8) (World Scientific Publishing Company, Singapore, 1997), pp. 92–111, [gr-qc/9710076].
  * (40) E. C. Vagenas,  Area spectrum of rotating black holes via the new interpretation of quasinormal modes, _J. High Energy Phys._ 11 (2008) 073 [gr-qc/0804.3264].
  * (41) R. Banerjee, B. R. Majhi, and E. C. Vagenas,  Quantum tunneling and black hole spectroscopy, _Phys. Lett. B_ 686 (2010) 279 [hep-th/0907.4271].
  * (42) B. Pourhassan, M. Faizal, Z. Zaz, and A. Bhat,  Quantum fluctuations of a BTZ black hole in massive gravity, _Phys. Lett. B_ 773 (2017) 325 [gr-qc/1709.09573].
  * (43) K. Lochan and S. Chakraborty,  Discrete quantum spectrum of black holes, _Phys. Lett. B_ 755 (2016) 37 [gr-qc/1509.09010].
  * (44) Y. Kwon and S. Nam,  Area spectra of the rotating BTZ black hole from quasinormal modes, _Classical Quantum Gravity_ 27 (2010) 125007 [hep-th/1001.5106].
  * (45) C. Rovelli,  Quantum Gravity (Cambridge University Press, Cambridge, England, 2004).
  * (46) H. Kastrup and T. Thiemann,  Spherically symmetric gravity as a completely integrable system, _Nucl. Phys. B_ 425 (1994) 665 [gr-qc/9401032].
  * (47) M. Campiglia, R. Gambini, and J. Pullin,  Loop quantization of spherically symmetric midi-superspaces, _Classical Quantum Gravity_ 24 (2007) 3649 [gr-qc/0712.0817].
  * (48) R. Gambini, J. Olmedo, and J. Pullin,  Quantum black holes in loop quantum gravity, _Classical Quantum Gravity_ 31 (2014) 095009 [gr-qc/1310.5996].
  * (49) R. Gambini and J. Pullin,  Loop Quantization of the Schwarzschild Black Hole, _Phys. Rev. Lett._ 110 (2013) 211301 [gr-qc/1302.5265].
  * (50) R. Balbinot, A. Fabbri, S. Fagnocchi, A. Recati, and I. Carusotto,  Nonlocal density correlations as a signature of Hawking radiation from acoustic black holes, _Phys. Rev. A_ 78 (2008) 021603 [cond-mat/0711.4520].
  * (51) I. Carusotto, S. Fagnocchi, A. Recati, R. Balbinot, and A. Fabbri,  Numerical observation of Hawking radiation from acoustic black holes in atomic Bose–Einstein condensates, _New J. Phys._ 10 (2008) 103001 [cond-mat/0803.0507].
  * (52) A. Fabbri and R. Balbinot,  Ramp-up of Hawking radiation in Bose-Einstein-condensate analog black holes, _Phys. Rev. Lett._ 126 (2021) 111301 [gr-qc/2005.04443].
  * (53) L. Liao, E.C.I. Van Der Wurff, D. Van Oosten, and H.T.C. Stoof,  Proposal for an analog Schwarzschild black hole in condensates of light, _Phys. Rev. A_ 99 (2019) 023850 [gr-qc/1806.00023].
  * (54) F. Michel, J-F. Coupechoux, and R. Parentani,  Phonon spectrum and correlations in a transonic flow of an atomic Bose gas, _Phys. Rev. D_ 94 (2016) 084027 [gr-qc/1605.09752].
  * (55) J. Steinhauer,  Observation of quantum Hawking radiation and its entanglement in an analogue black hole, _Nat. Phys._ 12 (2016) 959 [gr-qc/1510.00621].
  * (56) J.R.M. de Nova, K. Golubkov, V.I. Kolobov, and J. Steinhauer,  Observation of thermal Hawking radiation and its temperature in an analogue black hole, _Nature_ 569 (2019) 688 [gr-qc/1809.00913].
  * (57) R. Parentani,  From vacuum fluctuations across an event horizon to long distance correlations, _Phys. Rev. D_ 82 (2010) 025008 [gr-qc/1003.3625].
  * (58) R. Balbinot and A. Fabbri,  Quantum correlations across the horizon in acoustic and gravitational black holes, _Phys. Rev. D_ 105 (2022) 045010 [gr-qc/2107.00702].
  * (59) P.R. Anderson and J. Traschen,  Horizons and correlation functions in 2d Schwarzschild-de Sitter spacetime, _J. High Energy Phys._ 2022 (2022) 192 [hep-th/2012.08494].
  * (60) R. Balbinot and A. Fabbri,  Quantum energy momentum tensor and equal time correlations in a Reissner-Nordström black hole, _Phys. Rev. D_ 108 (2023) 045004 [gr-qc/2303.11039].
  * (61) M. Fontana and M. Rinaldi,  Stress-energy tensor correlations across regular black holes horizons, _Phys. Rev. D_ 108 (2023) 125003 [gr-qc/2302.08804].
  * (62) R. Balbinot and A. Fabbri,  The Hawking effect in the particles–partners correlations, _Physics_ 5 (2023) 968 [gr-qc/2306.16952].
  * (63) M. Banados, C. Teitelboim, and J. Zanelli,  Black hole in three-dimensional spacetime, _Phys. Rev. Lett._ 69 (1992) 1849 [hep-th/9204099].
  * (64) S. Carlip,  The (2+ 1)-dimensional black hole, _Classical Quantum Gravity_ 12 (1995) 2853 [gr-qc/9506079].
  * (65) A.R.H. Smith and R.B. Mann,  Looking inside a black hole, _Classical Quantum Gravity_ 31 (2014) 082001 [gr-qc/1309.4125].
  * (66) G. Lifschytz and M. Ortiz,  Scalar field quantization on the (2+ 1)-dimensional black hole background, _Phys. Rev. D_ 49 (1994) 1929 [gr-qc/9310008].
  * (67) L.J. Henderson, R.A. Hennigar, R.B. Mann, A.R.H. Smith, and J. Zhang,  Harvesting entanglement from the black hole vacuum, _Classical Quantum Gravity_ 35 (2018) 21LT02 [quant-ph/1712.10018].
  * (68) P. Painlevé,  C.R.Acad.sci. (paris) 173 (1921) 677.
  * (69) K. Gallock-Yoshimura, E. Tjoa, and R.B. Mann,  Harvesting entanglement with detectors freely falling into a black hole, _Phys. Rev. D_ 104 (2021) 025001 [quant-ph/2102.09573].
  * (70) A, Saini and D, Stojkovic,  Gravitational collapse and Hawking-like radiation of a shell in AdS spacetime, _Phys. Rev. D_ 97 (2018) 025020 [hep-th/1711.08182].
  * (71) P.K. Townsend and B. Zhang,  Thermodynamics of “exotic” Banados-Teitelboim-Zanelli black holes, _Phys. Rev. Lett._ 110 (2013) 241302 [hep-th/1302.3874].


# Communication: Non-adiabatic derivative coupling elements for the coupled
cluster singles and doubles model

Eirik F. Kjønstad<EMAIL_ADDRESS>Department of Chemistry and The
PULSE Institute, Stanford University, Stanford, California 94305, USA
Department of Chemistry, Norwegian University of Science and Technology, 7491
Trondheim, Norway  Henrik Koch  Department of Chemistry, Norwegian University
of Science and Technology, 7491 Trondheim, Norway  Scuola Normale Superiore,
Piazza dei Cavaleri 7, 56126 Pisa, Italy

###### Abstract

We present an efficient implementation of analytical non-adiabatic derivative
coupling elements for the coupled cluster singles and doubles model. The
derivative coupling elements are evaluated in a biorthonormal formulation in
which the nuclear derivative acts on the right electronic state, where this
state is biorthonormal with respect to the set of left states. This stands in
contrast to earlier implementations based on normalized states and a gradient
formula for the derivative coupling. As an illustration of the implementation,
we determine a minimum energy conical intersection between the $n\pi^{\ast}$
and $\pi\pi^{\ast}$ states in the nucleobase thymine.

††preprint: APS/123-QED

##  I Introduction

The nuclear dynamics that follows photoexcitation typically involves non-
adiabatic population transfer between several electronic states. For example,
in the nucleobase thymine, photoexcitation to the bright $\pi\pi^{\ast}$ state
is followed by rapid ($60$ fs) non-adiabatic population transfer to the dark
$n\pi^{\ast}$ state.[1] As is well known, the approximate description of the
electronic structure can have a dramatic qualitative impact on the simulated
nuclear dynamics, often complicating the task of correctly identifying the
actual physics behind the processes observed in pump-probe experiments.[2, 3]
A recent example is the ongoing debate about the dynamics that follows
excitation to the bright $B_{3u}$ state in pyrazine.[4, 5, 6, 7, 8] The
ambiguities involved in interpreting time-resolved spectra illustrate the need
for highly accurate description of the electronic structure.

A number of electronic structure methods has a long history of being applied
to treat non-adiabatic effects, including complete active space[9] (CAS)
methods, density functional theory[10] (DFT), and algebraic diagrammatic
construction[11] (ADC). These methods are often complementary, where some are
able to describe static correlation in the ground state and ground state
intersections (CAS) while others better capture dynamical correlation but are
unable to treat static correlation in the ground state as well as actual
crossings with the ground state (DFT, ADC). In the latter category, there is
still a need for a method that has systematically improvable accuracy that
extends beyond a perturbative description of double excitations.

Coupled cluster theory is now well-established as the method of choice
whenever this level of accuracy is required and the ground state is accurately
described by a single determinant. However, initial progress towards its use
in nonadiabatic dynamics simulations was slowed down with the realization[12,
13] that the method produces non-physical results at electronic degeneracies
when the states that cross span the same symmetry. Later work by the present
authors and collaborators showed that these artifacts were caused by the loss
of electronic state orthogonality (matrix defects)[14] and that they could be
fully removed by enforcing orthogonality relations between the electronic
states.[15, 16] Our current understanding is that coupled cluster methods are
able to describe conical intersections when the states span different
symmetries but corrections[13, 15, 16] are required when the states span the
same symmetry. However, these conclusions are based on studies of the
potential energy surfaces and not from considerations of the predicted
physics. It still remains an open question to what extent the artifacts at
same-symmetry intersections negatively affect the predicted dynamics in
trajectory-based simulation methods like surface hopping[17] and _ab initio_
multiple spawning.[18]

Already in 1999 Christiansen[19] derived expressions for the derivative
coupling elements in coupled cluster theory, but the first implementation was
given later by Tajti and Szalay[20] at the singles and doubles level (CCSD).
These authors did not, however, implement the expressions in Ref. 19. Instead,
the coupling was evaluated from the gradient of the two states as well as the
gradient of a fictitious summed state; this summed-state approach was also
used in a more recent implementation of the CCSD coupling elements.[21] In
addition, they proposed modifications to account for the fact that the coupled
cluster states are not normalized, building on earlier work by Gauss and
coworkers[22] who had found that normalization is important when evaluating
the diagonal Born-Oppenheimer correction to the energy. The need for
normalization in dynamics, which is not trivial to achieve, was later
questioned by Shamasundar.[23] In a recent publication, we confirmed this by
showing that a biorthonormal formalism exists in which there is no dependence
on the norm of the electronic states.[24]

In the present work, we provide a derivation (which is equivalent to Ref. 19)
and implementation, at the CCSD level of theory, of the derivative coupling
between ground and excited states as well as between excited states. The
derivation follows the Lagrangian approach for the derivative coupling
proposed by Hohenstein in the context of CAS configuration interaction
(CASCI),[25] while the present implementation builds on an efficient
implementation of analytical gradients, exploiting Cholesky decomposed
electronic repulsion integrals, recently published by the authors and
collaborators.[26]

##  II Theory

###  II.1 Lagrangian

The derivative coupling between states $i$ and $j$ is[19, 24]

$\displaystyle\bm{F}_{ij}=\langle\psi_{i}^{L}\hskip 1.00006pt|\hskip
1.00006pt\nabla\psi_{j}^{R}\rangle,\quad i,j=0,1,2,\ldots,$ (1)

where $L$ and $R$ signify that these are the left and right electronic states,
and the gradient $\nabla$ is taken with respect to the coordinates of the
atomic nuclei.

Analytical expressions for $\bm{F}_{ij}$ may be derived by using the
Lagrangian technique. Here, we use the Lagrangian proposed by Hohenstein.[25]
For the coupled cluster case, this Lagrangian can be expressed as[24]

$\displaystyle\mathscr{L}_{ij}=\mathscr{O}_{ij}+\text{conditions}$ (2)

where

$\displaystyle\mathscr{O}_{ij}=\langle\psi_{i}^{L}(\bm{x}_{0})\hskip
1.00006pt|\hskip 1.00006pt\psi_{j}^{R}(\bm{x})\rangle.$ (3)

Here we have made the dependence on the nuclear geometry explicit:
$\bm{x}_{0}$ is the geometry where the derivative is to be evaluated, while
$\bm{x}$ is allowed to vary. Upon differentiating $\mathscr{L}_{ij}$, the
derivative operation $\nabla$ only acts on the ket vector. As a result, the
derivative of $\mathscr{L}_{ij}$ at $\bm{x}_{0}$ is identical to $\bm{F}_{ij}$
at $\bm{x}_{0}$.[25, 24]

The conditions in $\mathscr{L}_{ij}$ are those that are required to specify
the right state $\psi_{j}^{R}$ for all values of $\bm{x}$. These are: the
Hartree-Fock equations, for specifying the orbitals; the amplitude equations,
for specifying the ground state cluster amplitudes; and the excited state
eigenvalue equations, for specifying the excited state amplitudes. Written out
in detail, the Lagrangian reads

$\displaystyle\begin{split}\mathscr{L}_{ij}&=\mathscr{O}_{ij}+\sum_{\mu}\bar{\zeta}_{\mu}\langle\mu\hskip
1.00006pt|\hskip 1.00006pt\bar{H}\hskip 1.00006pt|\hskip
1.00006pt\mathrm{HF}\rangle\\\
&+\sum_{\mu}\bar{\gamma}_{\mu}\bigl{(}\langle\mu\hskip 1.00006pt|\hskip
1.00006pt[\bar{H},R_{j}]\hskip 1.00006pt|\hskip
1.00006pt\mathrm{HF}\rangle-\omega_{j}R_{\mu}^{j}\bigr{)}\\\
&+\bar{\xi}(1-\langle L_{j}\hskip 1.00006pt|\hskip
1.00006ptR_{j}\rangle)+\sum_{ai}\bar{\kappa}_{ai}\mathscr{F}_{ai},\end{split}$
(4)

where we have suppressed the dependence on $\bm{x}_{0}$ for notational
convenience.

This expression for $\mathscr{L}_{ij}$ introduces various quantities. The
coupled cluster conditions are expressed in terms of the similarity-
transformed Hamiltonian

$\displaystyle\bar{H}=\exp(-T)\exp(\kappa)H\exp(-\kappa)\exp(T),$ (5)

where we have introduced the orbital rotation operator

$\displaystyle\kappa=\sum_{ai}\kappa_{ai}E_{ai}^{-},\quad
E_{ai}^{-}=E_{ai}-E_{ia},$ (6)

as well as the cluster operator

$\displaystyle T=\sum_{\mu}t_{\mu}\tau_{\mu}$ (7)

The scalars $t_{\mu}$ are known as cluster amplitudes, and the $\tau_{\mu}$
denote excitation operators. The $E_{ai}$ are singlet one-electron excitation
operators and $E_{ia}$ are corresponding deexcitation operators. Here,
$\kappa(\bm{x}_{0})=0$ by assumption.

The electronic states are conveniently expressed as

$\displaystyle|\psi_{k}^{R}\rangle$
$\displaystyle=\mathcal{R}_{k}\exp(T)|\mathrm{HF}\rangle$ (8)
$\displaystyle\langle\psi_{k}^{L}|$
$\displaystyle=\langle\mathrm{HF}|\mathcal{L}_{k}\exp(-T)$ (9)

where

$\displaystyle\mathcal{R}_{k}$
$\displaystyle=R_{0}^{k}+R_{k}=R_{0}^{k}+\sum_{\mu}R_{\mu}^{k}\tau_{\mu}$ (10)
$\displaystyle\mathcal{L}_{k}$
$\displaystyle=L_{0}^{k}+L_{k}=L_{0}^{k}+\sum_{\mu}L_{\mu}^{k}\tau_{\mu}^{\dagger}.$
(11)

We will also find it useful to write

$\displaystyle|\mathcal{R}_{k}\rangle$
$\displaystyle=R^{k}_{0}|\mathrm{HF}\rangle+|R_{k}\rangle$ (12)
$\displaystyle\langle\mathcal{L}_{k}|$
$\displaystyle=\langle\mathrm{HF}|L^{k}_{0}+\langle L_{k}|.$ (13)

Furthermore, we have let

$\displaystyle\omega_{k}=\langle L_{k}\hskip 1.00006pt|\hskip
1.00006pt[\bar{H},R_{k}]\hskip 1.00006pt|\hskip 1.00006pt\mathrm{HF}\rangle$
(14)

and defined the Fock matrix as

$\displaystyle\mathscr{F}_{pq}=h_{pq}+\sum_{k}(2g_{pkkq}-g_{pqkk}).$ (15)

Here, $h_{pq}$ and $g_{pqrs}$ are the one- and two-electron integrals of the
Hamiltonian. Following the conventional notation, we let $p,q,r,$ and $s$
denote generic orbitals; $i,j,k,$ and $l$ denote occupied orbitals; $a,b,c,$
and $d$ denote virtual orbitals. Lagrangian multipliers are denoted with a bar
($\bar{\zeta}_{\mu},\bar{\kappa}_{ai}$, $\bar{\xi}$).

The left-state quantities in $\mathscr{L}_{ij}$, that is, $\psi_{i}^{L}$ and
$L_{j}$, are _constants_ that define $\mathscr{L}_{ij}$. They are evaluated at
$\bm{x}_{0}$. Thus, the Lagrangian’s dependencies are understood as

$\displaystyle\mathscr{L}_{ij}=\mathscr{L}_{ij}(\bm{x},\bm{t},\bm{R}_{j},\bm{\kappa},\bar{\bm{\zeta}},\bar{\xi},\bar{\bm{\kappa}};\bm{x}_{0}),$
(16)

where the semicolon denotes that $\mathscr{L}_{ij}$ depends only
parametrically on $\bm{x}_{0}$.

###  II.2 Lagrangian stationarity conditions

The derivative coupling becomes the partial derivative of $\mathscr{L}_{ij}$
when the Lagrangian is stationary with respect to all variables and
multipliers that depend implicitly on $\bm{x}$. We begin by considering
stationarity for $\bm{R}_{j}$:

$\displaystyle\begin{split}\frac{\partial\mathcal{L}_{ij}}{\partial
R_{\sigma}^{j}}&=L^{i}_{\sigma}+\sum_{\mu}\bar{\gamma}_{\mu}A_{\mu\sigma}\\\
&-\omega_{j}\bar{\gamma}_{\sigma}-\sum_{\nu}L^{j}_{\nu}A_{\nu\sigma}\sum_{\mu}\bar{\gamma}_{\mu}R_{\mu}^{j}-\bar{\xi}L^{j}_{\sigma}=0,\end{split}$
(17)

where

$\displaystyle A_{\mu\nu}=\langle\mu\hskip 1.00006pt|\hskip
1.00006pt[\bar{H},\tau_{\nu}]\hskip 1.00006pt|\hskip
1.00006pt\mathrm{HF}\rangle.$ (18)

Using vector notation, this condition reads

$\displaystyle\bm{0}=\bm{L}_{i}^{T}+\bar{\bm{\gamma}}^{T}(\bm{A}-\omega_{j})-(\omega_{j}\bar{\bm{\gamma}}^{T}\bm{R}_{j}+\bar{\xi})\bm{L}_{j}^{T}$
(19)

Clearly, with $\bar{\xi}=-\omega_{j}\bar{\bm{\gamma}}^{T}\bm{R}_{j}$, the last
term in the equation vanishes, and we obtain stationarity provided

$\displaystyle\bar{\bm{\gamma}}^{T}=\frac{1}{\omega_{j}-\omega_{i}}\bm{L}_{i}^{T}.$
(20)

We thus see that the excited state multipliers ($\bar{\xi},\bar{\bm{\gamma}}$)
can be expressed in terms of the excited states ($\bm{L}_{i}$, $\bm{R}_{j}$)
and the associated excitation energies ($\omega_{i},\omega_{j}$).

Stationarity with respect to $\bm{t}$ yields

$\displaystyle\bm{0}={{}^{t}}\bm{\eta}^{T}+\bar{\bm{\zeta}}^{T}\bm{A}$ (21)

where

$\displaystyle{{}^{t}}\eta_{\sigma}=\langle\mathcal{L}_{i}\hskip
1.00006pt|\hskip 1.00006pt\tau_{\sigma}\hskip 1.00006pt|\hskip
1.00006pt\mathcal{R}_{j}\rangle+(\bm{\mathcal{F}}(\bar{\bm{\gamma}})\bm{R}_{j})_{\sigma},$
(22)

with the well-known[27] $F$-matrix defined as

$\displaystyle\mathcal{F}(\bar{\bm{\gamma}})_{\mu\nu}=\langle\bar{\gamma}\hskip
1.00006pt|\hskip 1.00006pt[[\bar{H},\tau_{\mu}],\tau_{\nu}]\hskip
1.00006pt|\hskip
1.00006pt\mathrm{HF}\rangle,\quad\langle\bar{\gamma}|=\langle\sigma|\bar{\gamma}_{\sigma}.$
(23)

Similarly, stationarity with respect to $\bm{\kappa}$ yields

$\displaystyle\bm{0}={{}^{\kappa}}\bm{\eta}^{T}+\bar{\bm{\kappa}}^{T}\bm{A}^{\mathrm{HF}},$
(24)

where

$\displaystyle\begin{split}{{}^{\kappa}}\eta_{ai}&=\langle\mathcal{L}_{i}\hskip
1.00006pt|\hskip 1.00006ptE_{ai}^{-}\hskip 1.00006pt|\hskip
1.00006pt\mathcal{R}_{j}\rangle+\langle\bar{\zeta}\hskip 1.00006pt|\hskip
1.00006pt[E_{ai}^{-},\bar{H}]\hskip 1.00006pt|\hskip
1.00006pt\mathrm{HF}\rangle\\\ &+\langle\bar{\gamma}\hskip 1.00006pt|\hskip
1.00006pt[[E_{ai}^{-},\bar{H}],R_{j}]\hskip 1.00006pt|\hskip
1.00006pt\mathrm{HF}\rangle,\end{split}$ (25)

and where $\bm{A}^{\mathrm{HF}}$ is the Hartree-Fock Hessian. The amplitude
and orbital conditions, given by Eqs. (21) and (24), are solved numerically
for $\bar{\bm{\zeta}}$ and $\bar{\bm{\kappa}}$.

###  II.3 Derivative coupling elements

Once $\bar{\bm{\zeta}}$ and $\bar{\bm{\kappa}}$ are known, we can evaluate the
coupling by taking the partial derivative of $\mathscr{L}_{ij}$ with respect
to the nuclear components $\\{q\\}$. This yields[19, 24]

$\displaystyle F_{ij}^{q}=\frac{\langle L_{i}\hskip 1.00006pt|\hskip
1.00006pt[\bar{H}^{q},R_{j}]\hskip 1.00006pt|\hskip
1.00006pt\mathrm{HF}\rangle}{\omega_{j}-\omega_{i}}+\langle\bar{\zeta}\hskip
1.00006pt|\hskip 1.00006pt\bar{H}^{q}\hskip 1.00006pt|\hskip
1.00006pt\mathrm{HF}\rangle+\bar{\kappa}_{ai}\mathscr{F}_{ai}^{q},$ (26)

where

$\displaystyle\bar{H}^{q}=\exp(-T)H^{q}\exp(T).$ (27)

Here $H^{q}$ denotes the partial derivative of $H$ with respect to the $q$th
nuclear coordinate, $x_{q}$. By expanding the commutator in Eq. (26), we
obtain the equivalent expression

$\displaystyle F_{ij}^{q}=\frac{\langle L_{i}\hskip 1.00006pt|\hskip
1.00006pt\bar{H}^{q}\hskip 1.00006pt|\hskip
1.00006ptR_{j}\rangle}{\omega_{j}-\omega_{i}}+\langle\tilde{\zeta}\hskip
1.00006pt|\hskip 1.00006pt\bar{H}^{q}\hskip 1.00006pt|\hskip
1.00006pt\mathrm{HF}\rangle+\bar{\kappa}_{ai}\mathscr{F}_{ai}^{q},$ (28)

where

$\displaystyle\tilde{\bm{\zeta}}=\bar{\bm{\zeta}}-\bm{J},\quad
J_{\mu}=\frac{\langle L_{i}\hskip 1.00006pt|\hskip 1.00006ptR_{j}\hskip
1.00006pt|\hskip
1.00006pt\mu\rangle}{\omega_{j}-\omega_{i}}=\frac{j_{\mu}}{\omega_{j}-\omega_{i}}.$
(29)

Clearly, $\bm{F}_{ij}$ is the sum of an excited state gradient and a ground
state gradient, plus an orbital relaxation term. The expression in Eq. (28) is
convenient when invoking an existing molecular gradient code.

So far we have assumed that the right state ($\psi_{j}^{R}$) is an excited
state. This raises the question of how to evaluate the coupling when
$\psi_{j}^{R}$ is the ground state ($j=0$). When this is the case, the excited
state condition in $\mathscr{L}_{ij}$ can be removed. As a result, the
$\bm{t}$ stationarity simplifies to

$\displaystyle\bm{0}=\bm{L}_{i}^{T}+\bar{\bm{\zeta}}^{T}\bm{A},$ (30)

so that

$\displaystyle\bar{\bm{\zeta}}^{T}=-\frac{1}{\omega_{i}}\bm{L}_{i}^{T}=\frac{1}{E_{0}-E_{i}}\bm{L}_{i}^{T},$
(31)

where $E_{k}$ denotes the electronic energy of the $k$th state. The orbital
multiplier equation is also simplified by the removal of
$\langle\bar{\gamma}|$, but this equation must still be solved numerically.
Once $\bar{\bm{\kappa}}$ is known, we can evaluate $\bm{F}_{i0}$ as

$\displaystyle F_{i0}^{q}=\frac{\langle L_{i}\hskip 1.00006pt|\hskip
1.00006pt\bar{H}^{q}\hskip 1.00006pt|\hskip
1.00006pt\mathrm{HF}\rangle}{E_{0}-E_{i}}+\bar{\kappa}_{ai}\mathscr{F}_{ai}^{q}.$
(32)

###  II.4 Significance of orbital connections

Hamiltonian derivatives are treated in the same way as for molecular energy
gradients. That is, we take $H$ to be expressed, for all $\bm{x}$, in a non-
unique orthonormal MO (OMO) basis which is defined by an orbital
connection.[28] Any orbital connection can be used, but the choice may
actually affect the expression for the derivative coupling. In fact, as we
will explain below, the formula in Eq. (26) is only correct when we use the
natural connection.[28, 19, 24] For other connections, such as the widely-used
symmetric connection, the partial derivative of $\mathscr{O}_{ij}$ is non-zero
and must be added to the expression for $\bm{F}_{ij}$.[25]

To show this, we express the derivative of $\mathscr{O}_{ij}$ in terms of the
orbital connection. Given a connection matrix $\bm{T}$, we define the OMOs as

$\displaystyle\psi_{p}=\sum_{q}T_{pq}\varphi_{q},$ (33)

where the unmodified MOs (UMOs) are given as

$\displaystyle\varphi_{q}=\sum_{\alpha}C_{\alpha
q}(\bm{x}_{0})\chi_{\alpha}(\bm{x}).$ (34)

Here, $\\{C_{\alpha q}\\}$ denotes MO coefficients, and $\\{\chi_{\alpha}\\}$
denotes atomic orbitals. The UMOs are generally only orthonormal at
$\bm{x}_{0}$, that is,

$\displaystyle S_{rs}=\langle\varphi_{r}\hskip 1.00006pt|\hskip
1.00006pt\varphi_{s}\rangle\neq\delta_{rs},\quad\bm{x}\neq\bm{x}_{0}.$ (35)

This is, of course, why an orbital connection is required in the first place;
consistently evaluating the derivative is most easily done in a Fock space
defined by an orbital basis that is orthonormal for all values of $\bm{x}$.

Figure 1: LiH/cc-pVQZ derivative coupling calculated for CCSD and FCI. For
CCSD, we present couplings both for the direct evaluation of the coupling
(this work) and for the summed-gradient values in Ref. 20.

Now, the derivative of $\mathscr{O}_{ij}$ can be written[28]

$\displaystyle\begin{split}\mathscr{O}_{ij}^{q}=\frac{\partial\mathscr{O}_{ij}}{\partial
x_{q}}\Big{|}_{0}=\sum_{rs}D_{rs}^{ij}Y_{rs}^{q},\end{split}$ (36)

where $\bm{D}^{ij}$ is the transition state density at $\bm{x}_{0}$, and

$\displaystyle Y_{rs}^{q}=\Big{\langle}\psi_{r}\hskip 1.00006pt\Big{|}\hskip
1.00006pt\frac{\partial\psi_{s}}{\partial x_{q}}\Big{\rangle}\Big{|}_{0}.$
(37)

For the natural connection, we have, by construction,[28]

$\displaystyle Y_{rs}^{q}=0,$ (38)

and so we can conclude that[28, 19]

$\displaystyle\mathscr{O}_{ij}^{q}=0.$ (39)

Next, let us consider the symmetric connection. In this case,
$\bm{T}=\bm{S}^{-1/2}$, which implies that

$\displaystyle\frac{\partial T_{rs}}{\partial
x_{q}}=-\frac{1}{2}\frac{\partial S_{rs}}{\partial
x_{q}}\Big{|}_{0}=-\frac{1}{2}(W_{rs}^{q}+W_{sr}^{q}),$ (40)

where

$\displaystyle W_{rs}^{q}=\Big{\langle}\varphi_{r}\hskip
1.00006pt\Big{|}\hskip 1.00006pt\frac{\partial\varphi_{s}}{\partial
x_{q}}\Big{\rangle}\Big{|}_{0}.$ (41)

Consequently,

$\displaystyle
Y_{rs}^{q}=W_{rs}^{q}-\frac{1}{2}(W_{rs}^{q}+W_{sr}^{q})=\frac{1}{2}(W_{rs}^{q}-W_{sr}^{q}),$
(42)

and so

$\displaystyle\begin{split}\mathscr{O}_{ij}^{q}&=\sum_{rs}D_{rs}^{ij}\Bigl{(}\frac{1}{2}(W_{rs}^{q}-W_{sr}^{q})\Bigr{)}\\\
&=\sum_{rs}\Bigl{(}\frac{1}{2}(D_{rs}^{ij}-D_{sr}^{ij})\Bigr{)}W_{rs}^{q}.\end{split}$
(43)

For the symmetric connection, therefore, the derivative of $\mathscr{O}_{ij}$
is equal to the anti-symmetrized density matrix contracted with a ket-
derivative of an overlap matrix.[25] This overlap derivative is evaluated as

$\displaystyle W_{rs}^{q}=\sum_{\alpha\beta}C_{\alpha r}C_{\beta
s}\Big{\langle}\chi_{\alpha}\hskip 1.00006pt\Big{|}\hskip
1.00006pt\frac{\partial\chi_{\beta}}{\partial x_{q}}\Big{|}_{0}\Big{\rangle}.$
(44)

For the natural connection, $\bm{W}^{q}$ is of course not needed for
$\mathscr{O}_{ij}^{q}$ (which is zero). However, $\bm{W}^{q}$ _is_ required
for the reorthonormalization terms associated with the Hamiltonian. For the
natural connection, the ket-derivative $\bm{W}^{q}$ plays the same role that
the braket-derivative $\bm{S}^{q}$ does for the symmetric connection.[28]
These reorthonormalization terms are the same for derivative couplings and
molecular energy gradients, so we refer the reader to the literature for more
details.[26]

###  II.5 Relation to previous implementations

In the literature, the derivative coupling has been implemented through a
summed-state formula[20, 21] which is closely related to the one presented in
this work. However, we have not been able to show that the two formulations
are equivalent, except in the FCI limit. As we will see, our values for the
coupling deviates to some extent from the values presented by Tajti and Szalay
for the LiH molecule.[20]

##  III Implementation

###  III.1 Evaluation of the derivative coupling

The derivative coupling has been implemented in a development version of the
$e^{T}$ program.[29] The implementation builds on the recent implementation by
Schnack-Petersen _et al_.[26] for ground and excited state molecular
gradients. Our implementation uses existing routines for molecular gradients
and two-electron densities,[26] as well as several other quantities already
implemented in the $e^{T}$ program,[29] such as the $F$-matrix
($\bm{\mathcal{F}}(\bar{\bm{\gamma}})$), the Hartree-Fock Hessian
($\bm{A}_{\mathrm{HF}}$), and the second and third terms of
${{}^{\kappa}}\bm{\eta}$. We apply central differences to obtain $\bm{W}^{q}$
numerically, exploiting Libint 2[30] to evaluate the AO overlap integrals.

We have implemented the first term in the ${{}^{t}}\bm{\eta}$ vector and in
the ${{}^{\kappa}}\bm{\eta}$ vector, that is, the terms that arise when
differentiating $\mathscr{O}_{ij}$ with respect to $\bm{t}$ and $\bm{\kappa}$.
In the case of CCSD, ${{}^{t}}\bm{\eta}$ can be expressed as

$\displaystyle\begin{split}{{}^{t}\eta}_{ai}^{1}&=L_{ai}^{i}R_{0}^{j}+j_{ai}\\\
&=L_{ai}^{i}R_{0}^{j}+\sum_{bj}L_{bj}^{i}R_{bjai}^{j}=D^{ij}_{ai}\end{split}$
(45) $\displaystyle{{}^{t}\eta}_{aibj}^{1}$
$\displaystyle=L_{aibj}^{i}R_{0}^{j},$ (46)

where $\bm{D}^{ij}$ is the one-electron transition density. Finally:

$\displaystyle{{}^{\kappa}\eta}_{ai}^{1}=D^{ij}_{ai}-D^{ij}_{ia}.$ (47)

We use an existing implementation to obtain the transition density
$\bm{D}^{ij}$.[29]

Finally, we have implemented the normalization factor $N_{j}^{L}$, since this
allows us to validate our implementation by comparison to the exact limit.
Programmable expressions for this quantity can be found elsewhere.[20]

Figure 2: Branching plane for CCSD/aug-cc-pVDZ conical intersection in H2S
($1{{}^{1}}A_{2}/1{{}^{1}}B_{1}$). We depict the relative electronic energies
(left) and the norm of the coupling vector (right). Figure 3: Branching plane
for CCSD/aug-cc-pVDZ conical intersections in HOF
($1{{}^{1}}A^{\prime\prime}/2{{}^{1}}A^{\prime\prime}$). We depict the real
part of relative electronic energies (left) and the norm of the coupling
vector (right).

###  III.2 Optimization of minimum energy conical intersections

As numerical illustrations of the new implementation, we have applied Bearpark
_et al_.’s algorithm for determining minimum energy conical intersections
(MECIs), where a gradient is constructed so that it is zero when two
conditions are fulfilled: the energy difference vanishes and the energy
gradient along the seam is zero.[31] In particular, we minimize the gradient

$\displaystyle\bm{G}=\mathcal{P}\nabla
E_{2}+2(E_{2}-E_{1})\frac{\bm{g}}{||\bm{g}||},$ (48)

where

$\displaystyle\bm{g}=\nabla(E_{2}-E_{1})$ (49)

and where $\mathcal{P}$ is the projection onto the complement of the
$\bm{g}$-$\bm{h}$ plane. The $\bm{h}$ vector is

$\displaystyle\bm{h}=(E_{2}-E_{1})\bm{F}_{12}.$ (50)

The gradient $\bm{G}$ is used in combination with a Broyden-Fletcher-Goldfarb-
Shanno (BFGS) solver already implemented in $e^{T}$ for geometry
optimizations.[26]

Figure 4: Branching plane for CCSD/cc-pVDZ Cs minimum energy conical
intersection in thymine ($n\pi^{\ast}$/$\pi\pi^{\ast}$). We depict the
relative electronic energies (left) and the norm of the coupling vector
(right).

##  IV Numerical examples

###  IV.1 Comparison to earlier implementation: LiH

In Figure 1, we show the derivative coupling element for the LiH system as a
function of the Li-H bond distance, computed with three methods: CCSD using
the direct formula (present work), CCSD using summed-state formula (numbers
taken from Tajti and Szalay[20]), and the exact FCI derivative couplings
(obtained with OpenMolcas[32]). All calculations are performed with the
Dunning[33] basis cc-pVQZ.

All three methods agree closely for all bond distances. However, there is a
slight deviation between our results and that given in Ref. 20, see Figure 1
(right). This may be caused by both insufficient numerical convergence (as
indicated by the uneven deviation from FCI) as well as differences in the
analytical derivative couplings, as noted in Section II.5.

In order to ensure a consistent comparison to FCI, where states are normalized
by default, we approximate the coupling from normalized coupled cluster
states, averaging over the left and right coupling elements:

$\displaystyle\begin{split}\bar{\bm{F}}_{ij}^{\mathrm{norm}}&=\frac{\bm{F}_{ij}^{\mathrm{norm}}-\bm{F}_{ji}^{\mathrm{norm}}}{2}\\\
&=\frac{\langle N_{i}^{L}\psi_{i}^{L}\hskip 1.00006pt|\hskip 1.00006pt\nabla
N_{j}^{R}\psi_{j}^{R}\rangle-\langle N_{j}^{L}\psi_{j}^{L}\hskip
1.00006pt|\hskip 1.00006pt\nabla N_{i}^{R}\psi_{i}^{R}\rangle}{2}\\\
&=\frac{N_{i}^{L}N_{j}^{R}\bm{F}_{12}-N_{j}^{L}N_{i}^{R}\bm{F}_{21}}{2}\\\
&\approx\frac{N_{i}^{L}(N_{j}^{L})^{-1}\bm{F}_{12}-N_{j}^{L}(N_{i}^{L})^{-1}\bm{F}_{21}}{2}\end{split}$
(51)

Recall that this normalization procedure is only required when we compare to
methods with normalized states.

###  IV.2 Branching planes in three-atomic systems: SH2, HOF

To provide some indication as to the behavior of the coupling in the vicinity
of conical intersections, we have calculated branching planes for points of
intersection in SH2 ($1{{}^{1}}A_{2}/1{{}^{1}}B_{1}$) and HOF
($1{{}^{1}}A^{\prime\prime}/2{{}^{1}}A^{\prime\prime}$); see Figures 2 and 3,
respectively. As expected, we find a divergence at the point of intersection
in SH2 and no visible artifacts. This is consistent with the fact that this is
an intersection between states spanning different symmetries.[14] The HOF
intersection, on the other hand, is defective because the states have the same
symmetry. Note that the coupling still diverges as one approaches the defect.

###  IV.3 Minimum energy conical intersection: thymine

Finally, we have applied the optimization algorithm described in Section III.2
to locate the $n\pi^{\ast}$/$\pi\pi^{\ast}$ minimum energy conical
intersection in thymine, restricted to nuclear geometries with $C_{s}$
symmetry; see Figure 4. In this calculation, we have used the cc-pVDZ basis.
As for SH2, this is a different-symmetry intersection and there is no sign of
non-physical artifacts.

##  V Summary and outlook

In this work we have presented an efficient implementation of derivative
coupling elements that will enable us to perform large-scale simulations of
nonadiabatic dynamics at the CCSD level of theory. Chemical systems of
interest are now within the reach of CCSD dynamics using _e.g._ the multiple
spawning framework;[18] for example, a single-point calculation on thymine
with a cc-pVDZ basis, including gradients of the $n\pi^{\ast}$ and
$\pi\pi^{\ast}$ states, as well as the coupling between them, can be performed
in a matter of minutes on a modern CPU node (see Schnack-Petersen _et al._[26]
for representative timings).

We emphasize that for systems where the intersecting states span the same
symmetry, the wavepacket may end up in regions that encompasses a defective
intersection. We then expect that corrections must be applied to the standard
CC methods in order to extract meaningful results, though this will depend on
the size of the defective intersection seam, which, in turn, depends on the
truncation level. Work on extending the present implementation to the
similarity constrained coupled cluster method (SCCSD), where such defects are
completely eliminated,[15, 16] is in progress. Note that the Lagrangian
approach makes such an extension straight-forward; we simply need to add the
orthogonality condition to the Lagrangian and solve the resulting response
equations.

The case that _can_ be treated with standard coupled cluster theory is that of
intersections where the states span different symmetries (_e.g._ the
$n\pi^{\ast}$ and $\pi\pi^{\ast}$ states in thymine). We may expect that such
systems can be accurately described in dynamics simulations where coupled
cluster theory provides the underlying electronic structure. This is the
subject of a forthcoming article.

###### Acknowledgements.

We thank David M. G. Williams for enlightening discussions. This work has
received funding from the European Research Council (ERC) under the European
Union’s Horizon 2020 Research and Innovation Programme (grant agreement No.
101020016). E.F.K. and H.K. both acknowledge funding from the Research Council
of Norway through FRINATEK project 275506. We acknowledge computing resources
through UNINETT Sigma2 – the National Infrastructure for High Performance
Computing and Data Storage in Norway, through project number NN2962k.

## References

  * Wolf  _et al._ [2017] T. Wolf, R. H. Myhre, J. Cryan, S. Coriani, R. Squibb, A. Battistoni, N. Berrah, C. Bostedt, P. Bucksbaum, G. Coslovich, _et al._ , “Probing ultrafast $\pi$$\pi$*/n$\pi$* internal conversion in organic chromophores via k-edge resonant absorption,” Nat. Commun. 8, 1–7 (2017).
  * Domcke, Yarkony, and Köppel [2011] W. Domcke, D. R. Yarkony, and H. Köppel,  _Conical intersections: theory, computation and experiment_ , Vol. 17 (World Scientific, 2011).
  * Curchod and Martinez [2018] B. F. Curchod and T. J. Martinez, “Ab initio nonadiabatic quantum molecular dynamics,” Chem. Rev. 118, 3305–3336 (2018).
  * Kanno  _et al._ [2015] M. Kanno, Y. Ito, N. Shimakura, S. Koseki, H. Kono, and Y. Fujimura, “Ab initio quantum dynamical analysis of ultrafast nonradiative transitions via conical intersections in pyrazine,” Phys. Chem. Chem. Phys. 17, 2012–2024 (2015).
  * Horio  _et al._ [2016] T. Horio, R. Spesyvtsev, K. Nagashima, R. A. Ingle, Y.-i. Suzuki, and T. Suzuki, “Full observation of ultrafast cascaded radiationless transitions from S2($\pi\pi^{\ast}$) state of pyrazine using vacuum ultraviolet photoelectron imaging,” J. Chem. Phys. 145, 044306 (2016).
  * Mignolet  _et al._ [2018] B. Mignolet, M. Kanno, N. Shimakura, S. Koseki, F. Remacle, H. Kono, and Y. Fujimura, “Ultrafast nonradiative transition pathways in photo-excited pyrazine: Ab initio analysis of time-resolved vacuum ultraviolet photoelectron spectrum,” Chem. Phys. 515, 704–709 (2018).
  * Sun  _et al._ [2020] K. Sun, W. Xie, L. Chen, W. Domcke, and M. F. Gelin, “Multi-faceted spectroscopic mapping of ultrafast nonadiabatic dynamics near conical intersections: A computational study,” J. Chem. Phys. 153, 174111 (2020).
  * Scutelnic  _et al._ [2021] V. Scutelnic, S. Tsuru, M. Pápai, Z. Yang, M. Epshtein, T. Xue, E. Haugen, Y. Kobayashi, A. I. Krylov, K. B. Møller, _et al._ , “X-ray transient absorption reveals the 1Au (n$\pi$*) state of pyrazine in electronic relaxation,” Nat. Commun. 12, 1–8 (2021).
  * Roos, Taylor, and Sigbahn [1980] B. O. Roos, P. R. Taylor, and P. E. Sigbahn, “A complete active space scf method (casscf) using a density matrix formulated super-ci approach,” Chem. Phys. 48, 157–173 (1980).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, “Self-consistent equations including exchange and correlation effects,” Phys. rev. 140, A1133 (1965).
  * Schirmer [1982] J. Schirmer, “Beyond the random-phase approximation: A new approximation scheme for the polarization propagator,” Phys. Rev. A 26, 2395 (1982).
  * Hättig [2005] C. Hättig, “Structure optimizations for excited states with correlated second-order methods: CC2 and ADC(2),” Advances in quantum chemistry 50, 37–60 (2005).
  * Köhn and Tajti [2007] A. Köhn and A. Tajti, “Can coupled-cluster theory treat conical intersections?” J. Chem. Phys. 127, 044105 (2007).
  * Kjønstad  _et al._ [2017] E. F. Kjønstad, R. H. Myhre, T. J. Martinez, and H. Koch, “Crossing conditions in coupled cluster theory,” J. Chem. Phys. 147, 164105 (2017).
  * Kjønstad and Koch [2017] E. F. Kjønstad and H. Koch, “Resolving the notorious case of conical intersections for coupled cluster dynamics,” The Journal of Phys. Chem. Lett. 8, 4801–4807 (2017).
  * Kjønstad and Koch [2019] E. F. Kjønstad and H. Koch, “An orbital invariant similarity constrained coupled cluster model,” J. Chem. Theory Comput. 15, 5386–5397 (2019).
  * Tully [1990] J. C. Tully, “Molecular dynamics with electronic transitions,” J. Chem. Phys. 93, 1061–1071 (1990).
  * Ben-Nun, Quenneville, and Martinez [2000] M. Ben-Nun, J. Quenneville, and T. J. Martinez, “Ab initio multiple spawning: Photochemistry from first principles quantum molecular dynamics,” J. Phys. Chem. A 104, 5161–5175 (2000).
  * Christiansen [1999] O. Christiansen, “First-order nonadiabatic coupling matrix elements using coupled cluster methods. i. theory,” J. Chem. Phys. 110, 711–723 (1999).
  * Tajti and Szalay [2009] A. Tajti and P. G. Szalay, “Analytic evaluation of the nonadiabatic coupling vector between excited states using equation-of-motion coupled-cluster theory,” J. Chem. Phys. 131, 124104 (2009).
  * Faraji, Matsika, and Krylov [2018] S. Faraji, S. Matsika, and A. I. Krylov, “Calculations of non-adiabatic couplings within equation-of-motion coupled-cluster framework: Theory, implementation, and validation against multi-reference methods,” J. Chem. Phys. 148, 044103 (2018).
  * Gauss  _et al._ [2006] J. Gauss, A. Tajti, M. Kállay, J. F. Stanton, and P. G. Szalay, “Analytic calculation of the diagonal Born-Oppenheimer correction within configuration-interaction and coupled-cluster theory,” J. Chem. Phys. 125, 144111 (2006).
  * Shamasundar [2018] K. Shamasundar, “Diagonal born–oppenheimer correction for coupled-cluster wave-functions,” Mol. Phys. 116, 1483–1495 (2018).
  * Kjønstad and Koch [2021] E. F. Kjønstad and H. Koch, “Biorthonormal formalism for nonadiabatic coupled cluster dynamics,” J. Chem. Theory Comput. 17, 127–138 (2021).
  * Hohenstein [2016] E. G. Hohenstein, “Analytic formulation of derivative coupling vectors for complete active space configuration interaction wavefunctions with floating occupation molecular orbitals,” J. Chem. Phys. 145, 174110 (2016).
  * Schnack-Petersen  _et al._ [2022] A. K. Schnack-Petersen, H. Koch, S. Coriani, and E. F. Kjønstad, “Efficient implementation of molecular CCSD gradients with Cholesky-decomposed electron repulsion integrals,” J. Chem. Phys. 156, 244111 (2022).
  * Koch and Jørgensen [1990] H. Koch and P. Jørgensen, “Coupled cluster response functions,” J. Chem. Phys. 93, 3333–3344 (1990).
  * Olsen  _et al._ [1995] J. Olsen, K. L. Bak, K. Ruud, T. Helgaker, and P. Jørgensen, “Orbital connections for perturbation-dependent basis sets,” Theor. Chem. Acc. 90, 421–439 (1995).
  * Folkestad  _et al._ [2020] S. D. Folkestad, E. F. Kjønstad, R. H. Myhre, J. H. Andersen, A. Balbi, S. Coriani, T. Giovannini, L. Goletto, T. S. Haugland, A. Hutcheson, I.-M. Høyvik, T. Moitra, A. C. Paul, M. Scavino, A. S. Skeidsvoll, Å. H. Tveten, and H. Koch, “eT 1.0: An open source electronic structure program with emphasis on coupled cluster and multilevel methods,” J. Chem. Phys. 152, 184103 (2020).
  * Valeev [2020] E. Valeev, “Libint: A library for the evaluation of molecular integrals of many-body operators over gaussian functions, version 2.7. 0-beta. 6,” (2020).
  * Bearpark, Robb, and Schlegel [1994] M. J. Bearpark, M. A. Robb, and H. B. Schlegel, “A direct method for the location of the lowest energy point on a potential surface crossing,” Chem. Phys. Lett. 223, 269–274 (1994).
  * Aquilante  _et al._ [2020] F. Aquilante, J. Autschbach, A. Baiardi, S. Battaglia, V. A. Borin, L. F. Chibotaru, I. Conti, L. De Vico, M. Delcey, I. Fdez. Galván, _et al._ , “Modern quantum chemistry with [open] molcas,” J. Chem. Phys. 152, 214117 (2020).
  * Dunning [1989] T. H. Dunning, “Gaussian basis sets for use in correlated molecular calculations. I. the atoms boron through neon and hydrogen,” J. Chem. Phys. 90, 1007–1023 (1989).