# Quantify the Non-Markovian Process with Intervening Projections in a
Superconducting Processor

Liang Xiang  Zhejiang Province Key Laboratory of Quantum Technology and
Device, Department of Physics, Zhejiang University, Hangzhou, 310027, China
Zhiwen Zong  Zhejiang Province Key Laboratory of Quantum Technology and
Device, Department of Physics, Zhejiang University, Hangzhou, 310027, China
Ze Zhan  Zhejiang Province Key Laboratory of Quantum Technology and Device,
Department of Physics, Zhejiang University, Hangzhou, 310027, China  Ying Fei
Zhejiang Province Key Laboratory of Quantum Technology and Device, Department
of Physics, Zhejiang University, Hangzhou, 310027, China  Chongxin Run
Zhejiang Province Key Laboratory of Quantum Technology and Device, Department
of Physics, Zhejiang University, Hangzhou, 310027, China  Yaozu Wu  Zhejiang
Province Key Laboratory of Quantum Technology and Device, Department of
Physics, Zhejiang University, Hangzhou, 310027, China  Wenyan Jin  Zhejiang
Province Key Laboratory of Quantum Technology and Device, Department of
Physics, Zhejiang University, Hangzhou, 310027, China  Zhilong Jia  Key
Laboratory of Quantum Information, University of Science and Technology of
China, Hefei, 230026, China  Peng Duan  Key Laboratory of Quantum Information,
University of Science and Technology of China, Hefei, 230026, China  Jianlan
Wu  Zhejiang Province Key Laboratory of Quantum Technology and Device,
Department of Physics, Zhejiang University, Hangzhou, 310027, China  Yi Yin
<EMAIL_ADDRESS>Zhejiang Province Key Laboratory of Quantum Technology and
Device, Department of Physics, Zhejiang University, Hangzhou, 310027, China
Guoping Guo<EMAIL_ADDRESS>Key Laboratory of Quantum Information,
University of Science and Technology of China, Hefei, 230026, China  Origin
Quantum Computing, Hefei, 230026, China

###### Abstract

A Markov assumption considers a physical system memoryless to simplify its
dynamics. Whereas memory effect or the non-Markovian phenomenon is more
general in nature. In the quantum regime, it is challenging to define or
quantify the non-Markovianity because the measurement of a quantum system
often interferes with it. We simulate the open quantum dynamics in a
superconducting processor, then characterize and quantify the non-Markovian
process. With the complete set of intervening projections and the final
measurement of the qubit, a restricted process tensor can be determined to
account for the qubit-environment interaction. We apply the process tensor to
predict the quantum state with memory effect, yielding an average fidelity of
$99.86\%\pm 1.1\permil$. We further derive the Choi state of the rest process
conditioned on history operations and quantify the non-Markovianity with a
clear operational interpretation.

##  I Introduction

Schrödinger equation can formulate the evolution of a quantum system whose
state is isolated. A realistic quantum system is unavoidably coupled to
environment, which typically contains many degrees of freedom and hard to
characterize and control [1, 2, 3]. The dynamics of this open quantum system
sometimes can be simplified with a Markov assumption, such as in the Gorini-
Kossakowski-Sudarshan-Lindblad (GKSL) master equation [4, 5, 6]. It assumes
that the environment is stationary and the system-environment interaction is
weak and uncorrelated. Mathematically, a Markov evolution generated by the
GKSL master equation is composed of a series of completely positive trace-
preserving (CPTP) maps, which form a semi-group [5, 7].

If the evolution of a quantum system depends not only on its instantaneous
state but also on history, it is considered non-Markovian or has memory
effects. The information of the process history can be stored somewhere and
retrieved later to affect the system. The non-Markovian behavior in the
quantum regime is challenging to define and quantify, mostly because that
general measurements of a quantum system will unavoidably interfere with it
[8, 9]. With the advance of quantum device and control apparatus, the
revealing of the non-Markovian effects have been constantly reported [10, 11,
12, 13, 14, 15]. Previous studies of non-Markovianity mainly characterize the
open dynamics with the final tomographic measurement, i.e. the observer only
waits at the end of the stochastic process to measure the mixed quantum state,
and the non-Markovianity is revealed by examining the state changes over a
series of time steps. In general, there are two kinds of approaches to
characterizing the non-Markovianity [16, 17, 11, 12, 13, 18, 19, 20, 21]. One
is to detect the information flow between the system and its environment. The
other is to check the divisibility of the dynamical map.

Although previous experiments have shown evidence of the quantum non-Markovian
process, there is no consensus on the quantification of memory effect [9, 22].
The observations also lacks a clear causal structure [23, 24]. Here in this
work, we apply an alternative method to characterize and quantify the non-
Markovian process with intervening projective measurements, or positive
operator-valued measurements (POVMs) [1]. The projections are not only the
final measurement but also active interrogations on each intermediate step of
the process. This method is feasible because that projective measurements can
be integrated as local operators in the tomographic process tensor framework
[25, 26]. Experimentally, we design a three-time-step open quantum process in
a superconducting processor, where standard two-qubit gates are selected to
simulate the system-environment interaction. We implement arbitrary POVMs by
using the field-programmable-gate-array (FPGA)-based fast measurement and
control hardware and customized quantum instruction set architecture (ISA)
[27]. The restricted process tensor is determined using a complete set of
POVMs. The intervening projective measurements reveal the information of the
instantaneous state, and at the same time refresh the system
deterministically, from which we can compare the causal relation on different
time steps.

Note that a process tensor has been experimentally determined with unitary
gates on IBM’s cloud-based quantum processors [28]. Their work traces the
information flow between different stages of the process, and provides the
lower bounding of the memory effect. Here since projective measurements
destroy the system-environment entanglement and steer the system into a
definite pure state that is independent of its previous trajectories, we can
check if there is a conditional dependence of the future dynamics on the past
control operations [9, 26]. Consequently, we can directly determine whether a
process is Markovian or non-Markovian in a finite number of experiments. With
projective measurements providing complete information of the system during
the process, we finally quantify the non-Markovian process in a smaller subset
of time steps conditioned on earlier operations. Based on our experiment, we
illustrate the operational interpretation of the non-Markovianity.

##  II Open Quantum Dynamics Described by the Tensor Network

Figure 1:  Tensor diagrams of the three-time-step process, including initial
states of the system $\mathrm{\mathbf{S}}$ (green triangle) and the
environment $\mathrm{\mathbf{E}}$ (orange triangle), intervening local
operators (blue box) and $\mathrm{\mathbf{SE}}$ interactions (red rectangle).
(a) Closed dynamics of a single qubit $\mathrm{\mathbf{S}}$. (b) Process
tensor depiction of the open quantum dynamics. The process tensor
$\mathbf{\mathcal{T}}^{2:0}$ encodes all the environmental information in the
dashed frame. (c) Traditional depiction of the open quantum dynamics,
represented by a concatenation of CPTP maps. For each step, the dynamics of
$\mathrm{\mathbf{S}}$ is simulated by the reduced $\mathrm{\mathbf{SE}}$
unitary evolution (outlined in the dashed frame).

We here use the mathematical tool of tensor network [29, 30] to describe the
quantum process. In a simplified situation [Fig. 1(a)], a complex vector
describes the pure state of a qubit system $\mathrm{\mathbf{S}}$. It can be
equivalently denoted by a tensor $(\psi_{0})_{i}$ with one index
$\textit{i}=$0 or 1, encoding the ground or excited state. At each time step
$n$, the system evolution can be represented by a tensor
$(\mathbf{u}_{n:n-1})^{j}_{k}$, while an observer intervenes in the trajectory
with a local operator $(\mathbf{a}_{n-1})^{i}_{j}$. The dynamics of
$\mathrm{\mathbf{S}}$ can be simulated by summing over the corresponding
tensor indexes, called a ‘tensor contraction’ [31].

For an open quantum system, we design a process where a single qubit
$\mathrm{\mathbf{S}}$ interacts with its environment [Fig. 1(b)]. The state of
$\mathrm{\mathbf{S}}$ under stochastic quantum evolution is now represented by
a density matrix $\rho^{\mathrm{S}}_{s_{0},r_{0}}$. In our superconducting
processor, we choose a neighboring ancilla qubit $\mathrm{\mathbf{E}}$ to
simulate the environment of $\mathrm{\mathbf{S}}$. Initially,
$\mathrm{\mathbf{S}}$ and $\mathrm{\mathbf{E}}$ are prepared to a tensor
product of the ground states,
$\rho^{\mathrm{SE}}_{s_{0},r_{0};\gamma_{0},\epsilon_{0}}=\rho^{\mathrm{S}}_{s_{0},r_{0}}\otimes\rho^{\mathrm{E}}_{\gamma_{0},\epsilon_{0}}=|00\rangle$.
The intervened local operator on $\mathrm{\mathbf{S}}$ is represented by a
tensor $\mathcal{A}^{s_{0},r_{0}}_{s^{\prime}_{0},r^{\prime}_{0}}$, and the
$\mathrm{\mathbf{SE}}$ interaction is represented by a tensor
$\mathcal{U}^{s^{\prime}_{0},r^{\prime}_{0};\gamma^{\prime}_{0},\epsilon^{\prime}_{0}}_{s_{1},r_{1};\gamma_{1},\epsilon_{1}}$.
This $n$-step quantum process can be simulated by contracting
$\rho^{\mathrm{\mathbf{SE}}}$ with sequential $\mathcal{A}$s and
$\mathcal{U}$s as following,

$\displaystyle\rho_{n}{(\mathbf{A}_{n-1:0})}$
$\displaystyle=\operatorname{Tr}_{\mathrm{E}}\left[\mathcal{U}_{n:n-1}\mathcal{A}_{n-1}\cdots\mathcal{U}_{1:0}\mathcal{A}_{0}\left(\rho_{0}^{\mathrm{SE}}\right)\right]$
(1) $\displaystyle=\mathcal{T}^{n:0}\left[\mathbf{A}_{n-1:0}\right].$ (2)

After tracing out the environment indexes at the final step [Eq. (1)], we
derive the output state of $\mathrm{\mathbf{S}}$. Consequently, there is a map
from the sequence of local operators
$\\{\mathcal{A}_{0},\cdots,\mathcal{A}_{n-1}\\}$ ($\mathbf{A}_{n-1:0}$ for
short) to the output state $\rho_{n}$. This map is coined a ‘quantum comb’ in
the quantum circuit architecture [32], or a ‘process tensor’ in open quantum
dynamics [25].

The process tensor encodes the hidden environmental information in the dashed
frame of Fig. 1(b). In Eq. (2), $\mathbf{\mathcal{T}}^{n:0}$ is a multilinear
map on $\mathbf{A}_{n-1:0}$, i.e. its linearity holds independently to the
operator $\mathcal{A}$ at each step. Then $\mathcal{T}^{n:0}$ can be
experimentally determined by the quantum tomography technique [33]. At step
$n$, the local operator can be uniquely decomposed by a fixed set of linearly
independent basis operators $\\{\mathcal{B}_{l}\\}$,
$\mathcal{A}_{n-1}=\sum_{l}(\beta_{n-1})_{l}\mathcal{B}_{l}$. The sequence of
local operations can be further expanded in terms of their tensor products,

$\displaystyle\mathbf{A}_{n-1:0}=\sum_{\vec{l}}\beta_{\vec{l}}\bigotimes_{m=0}^{n-1}(\mathcal{B}_{m})_{l}$
(3)

where real numbers $\beta_{\vec{l}}$ are the coefficients of the tensor
products of $\mathcal{B}s$, and $\vec{l}:=\\{l_{0},l_{1},\cdots,l_{n-1}\\}$
denotes the combination of operators.

A process tensor has been recently determined using a set of unitary control
gates. The tensor can be applied to characterize a non-Markovian process [28].
In this work, we probe the system with intervening projective measurements or
POVM operators. Experimentally, our FPGA-based hardware and customized quantum
ISA [27] enable us to apply an arbitrary POVM operator on any time step of the
process. These POVMs effectively destroy the entanglement between the system
and the environment, and steer the system to a definite pure state independent
of its previous trajectories [9]. With a complete set of POVM operators
$\\{\mathcal{P}\\}$, a ‘restricted’ process tensor $\mathbf{\mathcal{T}}_{P}$
can be determined [34]. We then derive the process tensor in a subset of
steps, which is ‘contained’ in a larger one (the containment property) [25].
With the smaller process tensor, we quantify the system’s non-Markovianity by
calculating the von Neumann mutual information [35] of its Choi state [36,
37], and give a clear operational interpretation in the context of our
experiment.

##  III Results

In our experiment, the superconducting quantum processor is a chip composed of
six cross-shaped Transmon qubits [38]. The device was mounted in a dilution
refrigerator whose base temperature is around 10 mK. We choose two neighboring
qubits as the system $\mathrm{\mathbf{S}}$ and the environment
$\mathrm{\mathbf{E}}$. They are set at $\omega_{\mathrm{S}}=6.21$ GHz and
$\omega_{\mathrm{E}}=5.70$ GHz. The other device parameters are the same as
that in Ref. [39]. We specifically design two distinct processes consisting of
standard quantum gates in the quantum processor. One has the
$\mathrm{\mathbf{SE}}$ interaction brought by a CNOT gate after the first
operator and CZ gate the second. We swap these two gates for the other
process, where CZ is activated first and CNOT the second.

###  III.0.1 Implementation of POVMs

The sub-normalized state of the quantum system after a POVM operator
$\mathcal{P}$ is given by $\mathcal{P}\rho\mathcal{P}^{\dagger}$, where
$\mathcal{P}=|p\rangle\langle p|$ also denotes the state to which we project.
The probability of this projection is $\operatorname{Tr}[\mathcal{P}\rho]$. An
important consequence of the POVM operator on an open quantum system is the
destruction of its entanglement with the environment, as
$(\mathcal{P}^{\mathrm{\mathbf{S}}}\otimes\mathcal{I}^{\mathrm{\mathbf{E}}})\rho^{\mathrm{\mathbf{SE}}}(\mathcal{P}^{\mathrm{\mathbf{S}}}\otimes\mathcal{I}^{\mathrm{\mathbf{E}}})=\mathcal{P}^{\mathrm{\mathbf{S}}}\otimes\mathcal{I}^{\mathrm{\mathbf{E}}}\operatorname{Tr}_{\mathrm{\mathbf{S}}}[(\mathcal{P}^{\mathrm{\mathbf{S}}}\otimes\mathcal{I}^{\mathrm{\mathbf{E}}})\rho^{\mathrm{\mathbf{SE}}}]$,
where $\mathcal{I}^{\mathrm{\mathbf{E}}}$ is the identity operator on the
environment state.

The qubit is usually read out in a superconducting processor by projecting it
to the ground state ($\mathcal{P}_{z+}=|0\rangle\langle 0|$) or excited state
($\mathcal{P}_{z-}=|1\rangle\langle 1|$). Whereas in general cases,
$\mathcal{P}$ shall project the qubit to any state on the Bloch sphere, which
is parameterized as
$|p\rangle=\mathrm{cos}\frac{\theta}{2}|0\rangle+e^{i\phi}\mathrm{sin}\frac{\theta}{2}|1\rangle$.
Any $\mathcal{P}(\theta,\phi)$ can be realized by introducing the following
unitary transformation,

$\displaystyle\mathcal{P}\rho\mathcal{P}^{\dagger}$
$\displaystyle=(R\mathcal{P}_{z+}R^{\dagger})\rho(R\mathcal{P}_{z+}R^{\dagger})^{\dagger}$
$\displaystyle=R(\mathcal{P}_{z+}\operatorname{Tr}[\mathcal{P}_{z+}R^{\dagger}\rho
R])R^{\dagger}$ (4)

where $R(\theta,\phi)$ is a $\theta$-angle rotation around the vector
$\vec{n}$ in the $xy$-plane. The azimuth of $\vec{n}$ is $\phi+\pi/2$. As an
example, the experimental sequence to implement a POVM operator
$\mathcal{P}_{y-}$ is drawn in the dashed frame in Fig. 2(a). We first rotate
the qubit by $+\pi/2$ around the $x$-axis. Then we apply a fast dispersive
measurement [40] to project it to either the ground state $\mathcal{P}_{z+}$
(the upper ball) or excited state $\mathcal{P}_{z-}$ (the lower ball with
lighter color). Finally, we rotate the qubit back to the target state
$|p\rangle$ (ball on the rightmost Bloch sphere along the -$y$-axis). Note
that in Eq. 4 the probability of $\mathcal{P}$ is calculated by post-selecting
[41] the ground-state ($\mathcal{P}_{z+}$). Similarly, we can implement other
POVMs to project the qubit onto any axis. $\mathcal{P}_{x+}$ denotes the
projection on the positive $x$-axis; $\mathcal{P}_{yz+}$ or
$\mathcal{P}_{zy+}$ denotes the projection on the internal ($\theta=-\pi/4$,
$\phi=0$) or exterior ($\theta=\pi/4$, $\phi=0$) angle bisector of +$y$ and
+$z$ axis; $\mathcal{P}_{yz-}$ is the projection on the reversed direction of
$\mathcal{P}_{yz+}$, and so on.

Figure 2:  (a) QPT protocol to characterize the $\chi$-matrix of
$\mathrm{\mathbf{S}}$, with the implementation of $\mathcal{P}_{y-}$ as an
example. (b) Left panel: $\chi$-matrix of $\mathcal{P}_{y-}$. The magnitude
and complex phase of $\chi$ elements are shown as the height and color of
cubic columns. Magnitudes lower than 0.02 are drawn in gray for visibility.
Right panel: QPT fidelities of nine POVMs. (c) $\chi$-matrix of the two-qubit
CZ gate. (d) $\chi$-matrices of the reduced CZ-gate, with
$\mathrm{\mathbf{E}}$ being prepared to the ground state (left), excited state
(middle), and $(|0\rangle-i|1\rangle)/\sqrt{2}$ (right).

###  III.0.2 Traditional Picture of the Quantum Map

The CPTP map of each step of the process can be equivalently expressed in
terms of the Kraus decomposition [1, 35, 42],

$\displaystyle\Lambda(\rho)=\sum E_{m}\rho E_{n}^{\dagger}\chi_{mn}$ (5)

where the set of matrices $\\{E\\}$ are called the Kraus operators of
$\Lambda$. To benchmark the quantum gates, QPT protocol are normally used to
determine the matrix $\chi_{mn}$ [1, 43]. Here we choose Pauli operators
$\\{I,\sigma_{x},\sigma_{y},\sigma_{z}\\}$ as the Kraus operators. The gate
fidelity is calculated by $\operatorname{Tr}[\chi_{\mathrm{ideal}}\chi]$,
where $\chi_{\mathrm{ideal}}$ represents an ideal matrix.

We apply the QPT protocol to characterize the POVMs mentioned in Sec. III.0.1.
In our experiment, we prepare six qubit states by rotations on the ground
state $|0\rangle$. After the process [dashed frame of Fig. 2(a)], the output
states are determined using quantum state tomography (QST) [1, 44]. Through
out the paper, we use six symmetric POVMs
($\mathcal{P}\\{$$x\pm$,$y\pm$,$z\pm\\}$) to determine the quantum state,
though four informationally complete (IC) POVMs are sufficient [45]. Knowing
the input-output relation on a complete set of basis, it is sufficient to
determine $\chi$. We can see four bars standing symmetrically on both side of
the diagonal of the matrix [left panel of Fig. 2(b)], representing
coefficients of[1/4, -1/4, -1/4, 1/4] for Kraus operators [$(I,I)$,
$(I,\sigma_{y})$, $(\sigma_{y},I)$, $(\sigma_{y},\sigma_{y})$]. Not
surprisingly, the quantum map we determine coincides with $\mathcal{P}_{y-}$,

$\displaystyle\Lambda_{\mathcal{P}_{y-}}[\rho]$
$\displaystyle\approx\frac{1}{4}(I\rho I-\sigma_{y}\rho
I-I\rho\sigma_{y}+\sigma_{y}\rho\sigma_{y})$
$\displaystyle=\frac{I-\sigma_{y}}{2}\rho\frac{I-\sigma_{y}}{2}=\mathcal{P}_{y-}\rho\mathcal{P}_{y-}.$

Note that while the projective measurement is performed along a deterministic
axis, the two complementary POVM elements it contains is probabilistic and may
not necessarily be trace preserving (TP).

We also characterize other POVMs, which will be used in Sec. III.0.3 as a full
basis of operators to determine the process tensor. After normalization, the
fidelity of $\mathcal{P}_{y-}$ is 96.9% $\pm$0.12%, averaged over 20
independent QPTs. Fidelities of other POVMs are plotted in the right panel of
Fig. 2(b). We have confirmed that the POVM operators are achieved with high
fidelities before intervening them in a quantum process to determine the
process tensor.

Similarly, we characterize the $\chi$ matrices of the CZ gate [Fig. 2(c)] and
CNOT gate (not shown). Their fidelities are $96.3\%\pm 0.29\%$ and $93.1\%\pm
0.53\%$, respectively. We next characterize the quantum map of
$\mathrm{\mathbf{S}}$ when it interacts with the ancilla
$\mathrm{\mathbf{E}}$, termed the ‘reduced’ quantum map of
$\mathrm{\mathbf{S}}$. For CZ gate, the ‘reduced’ map of $\mathrm{\mathbf{S}}$
depends on the state of $\mathrm{\mathbf{E}}$, which can be written as,

$\displaystyle\Lambda_{\operatorname{Tr}_{\mathrm{\mathbf{E}}}[\mathrm{CZ}]}(\rho^{\mathrm{\mathbf{S}}})=\operatorname{Tr}_{\mathrm{\mathbf{E}}}[\mathcal{U}_{\mathrm{CZ}}\rho^{\mathrm{\mathbf{SE}}}\mathcal{U}^{\dagger}_{\mathrm{CZ}}].$
(6)

Figure 2(d) displays the $\chi$-matrices of three ‘reduced’ CZ gates with
$\mathrm{\mathbf{E}}$ in the ground state ($|0\rangle$), the excited state
($|1\rangle$), and the superposition state
($(|0\rangle-i|1\rangle)/\sqrt{2}$). The equivalent quantum maps are an
identity ($I$), a $180^{\circ}$ rotation around the $z$-axis ($Z$), and a
complete phase erasing (mixture of $I$ and $Z$ with equal probability),
respectively.

###  III.0.3 Tomography of the Restricted Process Tensor

Following we choose a complete set of basis POVM operators [34],

$\displaystyle\mathcal{F}=\mathcal{P}{\\{x+,x-,y+,z+,y-,z-,xy+,xz+,yz+\\}}.$
(7)

By performing the combinations $\\{\mathrm{\mathbf{B}}_{1:0}\\}$ =
$\bigotimes^{2}\\{\mathcal{P}\\}$ at the first two steps of the process and
measuring the outcome states of $\mathrm{\mathbf{S}}$ using QST, the
restricted process tensor can be determined as the following algorithm,

Algorithm Process Tomography of $\mathbf{\mathcal{T}}_{P}^{2:0}$

1:Stats = 3000 $\triangleright$ Repetitions of the sequence

2:$\mathrm{\mathbf{U}}_{2:0}$ = CNOT-CZ or CZ-CNOT $\triangleright$
$\mathrm{\mathbf{SE}}$ interactions

3:$\\{\mathrm{\mathbf{B}}_{1:0}\\}=\bigotimes^{2}\\{P_{x+},P_{x-},P_{y+},\ldots,P_{yz+}\\}_{9\times
9}$ $\triangleright$ three-step Basis

4:for $\mathcal{B}_{l_{0}},\mathcal{B}_{l_{1}}$ in
$\\{\mathrm{\mathbf{B}}_{1:0}\\}$ do

5: for $R_{m}$ in {I, X/2, Y/2, -X/2, -Y/2, X} do $\triangleright$ QST

6: for i in range(Stats) do

7: Apply $\mathcal{B}_{l_{0}}$ to $\mathrm{\mathbf{S}}$

8: Apply $\mathcal{U}_{0}$ to $\mathrm{\mathbf{SE}}$ $\triangleright$ 1st
$\mathrm{\mathbf{SE}}$ interaction

9: Apply $\mathcal{B}_{l_{1}}$ to $\mathrm{\mathbf{S}}$

10: Apply $\mathcal{U}_{1}$ to $\mathrm{\mathbf{SE}}$ $\triangleright$ 2nd
$\mathrm{\mathbf{SE}}$ interaction

11: Apply $\mathrm{R}_{m}$ to $\mathrm{\mathbf{S}}$

12: Record the system state of this sequence $S_{\vec{l},m,i}$

13: end for

14: Compute the state probability $P_{\vec{l},m}$ from $\\{S_{\vec{l},m,i}\\}$

15: end for

16: Compute the QST result $(\rho_{\vec{l}})_{2\times 2}$ from
$\\{P_{\vec{l},m}\\}$

17:end for

18:Find the solution $(\mathbf{\mathcal{T}}_{P}^{2:0})_{256\times 4}$ as a
linear map from $(\mathrm{\mathbf{B}}_{{\vec{l}}})_{9\cdot 9\times 4^{2}\cdot
4^{2}}$ to $(\rho_{{\vec{l}}})_{9\cdot 9\times 2\cdot 2}$

Note we fit the restricted process tensor using the least-square method
provided by scientific computing packages [46], and then obtain the best
approximation $\mathbf{\mathcal{T}}_{P}$ with positive semi-define constraint
and convex optimization (see Sec. III.0.4).

Once a process tensor $\mathbf{\mathcal{T}}_{P}^{2:0}$ is determined, we can
predict its output state when arbitrary projections
$\mathrm{\mathbf{A}}_{1:0}$ are inserted in the process. The fidelity of
predicted state is defined as [47],

$\displaystyle
F(\rho,\sigma)=(\operatorname{Tr}\sqrt{\sqrt{\rho}\sigma\sqrt{\rho}})^{2}$ (8)

where $\rho$ is the density matrix of $\mathrm{\mathbf{S}}$ measured in
experiment, and $\sigma$ is the predicted one. We use an over complete set of
$\mathrm{\mathbf{A}}_{1:0}$, consisting of POVMs in $\mathcal{P}\\{x+$, $x-$,
$y+$, $z+$, $y-$, $z-$, $xy+$, $xz+$, $yz+$, $xy-$, $xz-$, $yz-$, $yx+$,
$zx+$, $zy+$, $yx-$, $zx-$, $zy-\\}$, to benchmark the process tensor
prediction. For comparison, we also use the following traditional quantum map
method to predict the quantum process [Fig. 1(c)].

Method Quantum Map of the Markovian System

1:Assume that $\mathrm{\mathbf{S}}$ is initially at the ground state.

2:Calculate the state of $\mathrm{\mathbf{S}}$ after the CP map of
$\mathcal{A}_{0}$ using the $\chi$-matrix determined previously [Fig. 2(b)].

3:Calculate the state of $\mathrm{\mathbf{S}}$ after $\mathcal{U}_{1:0}$ using
the $\chi$-matrix of the reduced quantum map determined previously.

4:Repeat previous two steps for $\mathcal{A}_{1}$ and $\mathcal{U}_{2:1}$, to
obtain the final output state.

Note that in step 3, we choose the reduced map of $\mathcal{U}_{n:n-1}$
conditioned on $\mathrm{\mathbf{E}}$ being in $|0\rangle$. For example, the
$\chi$-matrix on the left panel of Fig. 2(d). We assume that
$\mathrm{\mathbf{E}}$ remains in $|0\rangle$ throughout the process, which is
the Markov assumption widely used to deal with the open quantum dynamics but
does not always hold.

In the CNOT-CZ process where CNOT gate is activated at the first step and CZ
at the second step, the output states predicted by the process tensor yields
extremely high fidelity and stability, averaged to $99.86\%\pm 1.1\permil$
over 20 repetitions. While the traditional quantum map method with the Markov
assumption can not accurately predict the output of the process, whose average
fidelity is only $80.25\%\pm 13.2\permil$. A second CZ-CNOT process is also
characterized for reference. Both methods can well predict the process, whose
average fidelities are $99.87\%\pm 1.0\permil$ (process tensor) and
$99.76\%\pm 1.5\permil$ (quantum map of the Markovian system).

Figure 3:  Prediction fidelities of the CNOT-CZ process. We show three
representative groups where qubit is firstly projected to three states in the
$zy$-plane ($\phi=0$), (a) $\mathcal{A}_{0}=\mathcal{P}_{z+}$ ($\theta=0$),
(b) $\mathcal{A}_{0}=\mathcal{P}_{zy+}$ ($\theta=\pi/4$), (c)
$\mathcal{A}_{0}=\mathcal{P}_{y-}$ ($\theta=\pi/2$). Prediction fidelities of
process tensor method (blue dot) and traditional quantum map method (solid
bar) are drawn on the horizontal axis for different second-step projections
$\mathcal{A}_{1}$. Theoretical result of the second method is calculated with
all gates ideal (black line). Note that one fidelity is omitted when
projections are chosen as $\mathcal{A}_{1:0}$ =
$\mathcal{P}_{z+}$-$\mathcal{P}_{z-}$ because the measurement probability is
close to zero. In other words, this trajectory is forbidden in this process.

To study how the previous trajectory of the system affects the subsequent
dynamics, partial results of the CNOT-CZ process are presented in Fig. 3. The
fidelity of predicted states are grouped by different choices of
$\mathcal{A}_{0}$: $\mathcal{P}_{z+}$[Fig. 3(a)], $\mathcal{P}_{zy+}$[Fig.
3(b)], $\mathcal{P}_{y-}$[Fig. 3(c)]. We can distinguish the non-Markovian
trajectories by checking the discrepancies of the fidelity. For example, if we
fix the second-step projection to $\mathcal{P}_{x+}$, the Markov predictions
(quantum map method) are unchanged because the reduced CZ gate is actually an
identity map on $\mathrm{\mathbf{S}}$ [left panel of Fig. 2(d)]. However, the
distances between the unchanged state ($\sigma$=$\mathcal{P}_{x+}$) and the
three experimental measured $\rho$ are different, indicating the change of the
final state. Clearly, the second-step trajectory of $\mathrm{\mathbf{S}}$
depends on the the first-step operation $\mathcal{A}_{0}$, a distinguishing
feature of non-Markovianity. In other words, the system dynamics is affected
by the history of operations on it. For the CZ-CNOT process, however, we can
not find any evident changes of the second-step quantum trajectory conditioned
on different $\mathcal{A}_{0}$. This means that the CZ-CNOT process is almost
Markovian. Note that the visibility of non-Markovianity in experiment (solid-
red bar) is lower than the ideal one (black reference line). This is because
that both $\mathrm{\mathbf{S}}$ and $\mathrm{\mathbf{E}}$ in our processor
dissipate to a larger lossy environment. For qubit $\mathrm{\mathbf{S}}$, the
measured state and the predicted state shrink to the north pole of the Bloch
sphere, leading to a closer distance in between. For qubit
$\mathrm{\mathbf{E}}$, it relaxes to the ground state during the process,
losing some memory.

###  III.0.4 Quantify the Non-Markovianity

We further quantify the non-Markovianity of the process. With recorded results
of a complete set of POVMs on each step (Sec. III.0.3), we can vary the first-
step POVM operator $\mathcal{P}(\theta,0)$ and derive the rest process tensor,
parameterized as $\mathbf{\mathcal{T}}_{P}^{2:1}(\theta)$. The azimuth angle
$\phi$ is set to zero, i.e. all states of $\mathrm{\mathbf{S}}$ after the
first projections are in the $zy$ plane.

The non-Markovian evolution of a system are held in the process tensor
$\mathbf{\mathcal{T}}_{P}$, which can be conveniently rewritten as a many-body
generalized Choi state $\Upsilon$ using the ‘Choi-Jamiołkowski’ representation
[37]. Here we choose to gauge the non-Markovianity by the von Neumann mutual
information [35] in $\Upsilon$, which can also be viewed as the distance
between $\Upsilon$ and its uncorrelated state $\Upsilon_{\mathrm{\otimes}}$.
The distance is calculated by the von Neumann relative entropy,

$\displaystyle\mathcal{D}(\Upsilon\|\Upsilon_{\mathrm{\otimes}})=\operatorname{Tr}\left[\Upsilon(\operatorname{ln}\Upsilon-\operatorname{ln}\Upsilon_{\mathrm{\otimes}})\right].$
(9)

In our experiment, the Choi state of the last-step process is denoted by
$(\Upsilon^{2:1})_{4^{1}\cdot 2\times 4^{1}\cdot 2}$, with subscripts for its
dimensions. The uncorrelated state $\Upsilon^{2:1}_{\mathrm{\otimes}}$
represents the last-two-step process with a Markov assumption [Fig. 2(c)]. We
calculate it as the tensor product of the average initial state of
$\mathrm{\mathbf{S}}$ and the Choi state of the reduced quantum map
(conditioned on the average state of $\mathrm{\mathbf{E}}$),

$\displaystyle(\Upsilon^{2:1}_{\mathrm{\otimes}})_{4\cdot 2\times 4\cdot
2}=(\Lambda^{2:1})_{4\times 4}\otimes(\rho_{1})_{2\times 2}.$ (10)

Care should be taken when deriving the restricted process tensor
$\mathbf{\mathcal{T}}_{P}^{2:1}$, because it is under-determined using only
POVMs. The solutions may not always be positive semi-define (PSD). In
contrast, on the set of complete basis including both projective measurement
and unitary control, the full process tensor is well-determined and PSD [25].
Here we use an optimization method to minimize the non-Markovianity of
$\Upsilon^{2:1}$ in the space spanned by POVMs with the PSD constraint. The
procedure to derive $\mathcal{N}^{2:1}$ is illustrated in the left panel of
Fig. 4(a). The restricted process tensor determined by projective measurements
is defined in a linear space spanned by the POVM operators
($\mathrm{\mathbf{PM}}$, dashed-blue line). It contains a smaller space of
process tensors whose Choi states are PSD (solid-orange area). In the
optimization method, we first find an initial solution $\Upsilon^{2:1}_{0}$
(black dot) using the least-square fitting, and then progressively optimize
the Choi state $\Upsilon^{2:1}$ by minimizing the distance [Eq. 9] to its
uncorrelated state (cross circle $\otimes$) [Eq. 10]. We take the minimal
distance as the measure of non-Markovianity $\mathcal{N}^{2:1}$,

$\displaystyle\mathcal{N}^{2:1}=\min_{\Upsilon^{2:1}\in{\\{PSD\\}}}\mathcal{D}(\Upsilon^{2:1}\|\Upsilon^{2:1}_{\mathrm{\otimes}}).$
(11)

The result of Eq. 11 converges to the point (black cross) that has the minimal
distance to $\Upsilon^{2:1}_{\mathrm{\otimes}}$.

Figure 4:  (a) Left panel: space where the last-two-step restricted process
tensor is defined. The arrow shows the optimization procedure, starting from
an initial value (black dot), to one point $\Upsilon^{2:1}$ (black cross) that
is closet to the uncorrelated Choi state $\Upsilon^{2:1}_{\mathrm{\otimes}}$
(cross circle). Right panel: the non-Markovianity of the last-two-step process
conditioned on different first-step POVMs parameterized with $\theta$. The
CNOT-CZ process yields the maximum non-markovianity when $\theta=\pi/2$. The
experimental result (red dashed line) is lower than the ideal case (red solid
line). Numerical simulations considering dissipation of both
$\mathrm{\mathbf{S}}$ and $\mathrm{\mathbf{E}}$ (red dotted line) accounts for
the trend. The other CZ-CNOT process is almost Markovian and memoryless (blue
lines). (b) Volumes of the accessible states of $\mathrm{\mathbf{S}}$. At the
beginning of the second step of CNOT-CZ process, POVM operators spans the
states all around the Bloch sphere (left volume). When the first-step
projection $\mathcal{A}_{0}$ is $\mathcal{P}_{z+}$, the final state of the
process is unchanged (Identity to $\mathcal{A}_{1}$), well predicted by both
methods. When $\mathcal{A}_{0}$ is $\mathcal{P}_{zy+}$ or $\mathcal{P}_{y-}$,
the process tensor $\mathbf{\mathcal{T}}_{2:1}$ and the quantum map
$\mathbf{\Lambda}_{2:1}$ generate the output of the process differently. We
color the states according to the projection angle $\theta$ of
$\mathcal{A}_{1}$.

The right panel of Fig. 4(b) shows the results of $\mathcal{N}^{2:1}$ versus
different choices of the first-step operator ${\mathcal{P}(\theta,0)}$. For
the CNOT-CZ process, it yields the maximum non-Markovianity when
$\theta$=$\frac{\pi}{2}$ (red dashed line). On the contrary, the CZ-CNOT
process does not produce significant non-Markovianity (blue dashed line)
whatever $\mathcal{A}_{0}$ we choose. As we conclude in Sec. III.0.3, the CZ-
CNOT process is almost Markovian.

Note that non-Markovianity of the CNOT-CZ process in the real quantum
processor (red dashed line) shows a lowered visibility compared with the ideal
theoretical result (red solid line). Similar to the observation in Sec.
III.0.3, this lowered visibility is due to the dissipation of both
$\mathrm{\mathbf{S}}$ and $\mathrm{\mathbf{E}}$ to a larger environment. We
numerically calculate $\mathbf{\mathcal{T}}_{P}^{2:0}$ using the non-ideal
CNOT and CZ gates that we have experimentally characterized, and then derive
the $\mathcal{N}^{2:1}$s versus $\theta$ (red dotted line), from which the
trend can be partially verified. The even lower value of $\mathcal{N}^{2:1}$
is most likely caused by other noises, such as the phase noise of a ‘mediocre’
clock when qubit sequence gets longer [48], and the photon number fluctuations
noise during the projective measurement [49, 50].

Both processes show a larger standard deviation when $\theta$ gets bigger.
This is because that we have normalized the Choi state by dividing the
probability of $\mathcal{A}_{0}$, which reduces to 0 as $\theta$ gets close to
$\pi$. $\mathcal{N}^{2:1}$s of both processes also get higher when $\theta$
approaches to $\pi$. This indicates that some memory sources are not included
in the simple two-qubit model when $\mathrm{\mathbf{S}}$ is projected to a
higher population of excited-state. Some possible candidates are the
measurement induced state transitions [51] and the crosstalk from neighboring
qubit. Nevertheless, such temporal correlation effects will be absorbed in the
process tensor, which can accurately predict the non-Markovian process (Sec.
III.0.3).

To illustrate the operational meaning of $\mathcal{N}^{2:1}$, we analyze the
cause of non-Markovianity in the CNOT-CZ process. At the beginning of the
first process, $\mathrm{\mathbf{E}}$ is reset to $|0\rangle$. After one of the
first projections $\mathcal{P}_{z+}$, $\mathcal{P}_{zy+}$, or
$\mathcal{P}_{y-}$, the state of $\mathrm{\mathbf{S}}$ changes to either
$|0\rangle$,
$\mathrm{cos}\frac{\pi}{8}|0\rangle-i\mathrm{sin}\frac{\pi}{8}|1\rangle$, or
$(|0\rangle-i|1\rangle)/{\sqrt{2}}$. Correspondingly, the state of
$\mathrm{\mathbf{SE}}$ changes after the CNOT gate to either a product state
$|00\rangle$, a partial entangled state
$\mathrm{cos}\frac{\pi}{8}|00\rangle-i\mathrm{sin}\frac{\pi}{8}|11\rangle$, or
a maximum entangled state (MES) $(|00\rangle-i|11\rangle)/\sqrt{2}$. As has
been studied in Refs. [34, 52, 11, 26], initial correlations between a system
$\mathrm{\mathbf{S}}$ and its environment $\mathrm{\mathbf{E}}$ carries
historical information and will affect its dynamics at a later time. At the
beginning of the second step, the complete set of POVMs $\mathcal{A}_{1}$
(after normalization) effectively spans the state all around the Bloch sphere,
[left volume in Fig. 4(b)]. The process tensor and the quantum map method
describe the following process differently. We present the volume changes of
the set of accessible states of $\mathrm{\mathbf{S}}$ after the CZ gate in
fig. 4(b). We first check the quantum map method. When the
$\mathrm{\mathbf{SE}}$ state is $|00\rangle$, the reduced CZ gate is
effectively an identity gate, leaving the states unchanged (left volume). When
the $\mathrm{\mathbf{SE}}$ state is partially entangled
($\mathrm{cos}\frac{\pi}{8}|00\rangle-i\mathrm{sin}\frac{\pi}{8}|11\rangle$),
some of the phase information of $\mathrm{\mathbf{S}}$ is lost statistically
after the quantum map and the volume of states shrinks toward the $z$-axis
(middle lower volume). An extreme case is when $\mathrm{\mathbf{SE}}$ is
initially in MES. Phase information of $\mathrm{\mathbf{S}}$ is completely
erased through the map, and all output states goes to the $z$-axis (right
lower volume). Quite differently, the process tensor ‘knows’ how
$\mathrm{\mathbf{S}}$ and $\mathrm{\mathbf{E}}$ are correlated after
$\mathcal{A}_{1}$, which breaks the entanglement and projects
$\mathrm{\mathbf{S}}$ to a predefined state. Thus the last-two-step process is
more accurately described with tensor operators. The volume of accessible
state predicted by the process tensor [right upper volumes in Fig. 4(b)]
differs most with that of the quantum map when $\mathrm{\mathbf{SE}}$ is
initially in MES. This, in turn, is consistent with the non-Markovianity
obtained in Fig. 4(a). We will most unlikely to confuse the last-two step
process to be Markovian, if we apply $\mathcal{P}_{y-}$ at the first step and
compare the output states with the Markov predictions. The theoretical value
(red solid line) approaching $\ln{2}$ means that the maximum probability of
not finding the process non-Markovian would be $e^{-\mathcal{N}}=0.5$ for
every single ideal experiment.

For the quantification of non-Markovianity, the data analysis code is
available in the repository of Github [53].

##  IV Conclusion

In conclusion, we experimentally quantify the non-Markovianity of a quantum
system by intervening projective measurements. The restricted process tensor
is determined with a complete set of POVMs. Compared with the traditional
quantum map method, the process tensor leads to remarkably high fidelities in
predicting the output of an open quantum process with or without memory
effect. The non-Markovianity of a subset of the process is quantified
conditioned on the choices of the first-step projection. For the CNOT-CZ
process, we unambiguously determine the exsitance of non-Markovianity and show
that the memory effect is rooted in the spacial correlation between the system
and its environment. Based on the experiment, we illustrate the operational
meaning of the non-Markovianity: as the non-Markovianity goes high, there is
an increased likelihood to find the Markov assumption wrong.

Although an ancilla qubit is used in our work to simulate the environment, the
process tensor itself is an inclusive model to represent the non-Markovian
noise stemmed from a wide range of microscopic mechanisms [54, 55, 56, 57, 58,
59, 10, 15]. The process tensor method can identify the non-Markovian noise
when the experimenter actively intervenes with the quantum evolution by either
measurement or control. It will be helpful to analyze and quantify the non-
Markovian noise environment in larger quantum processor [60]. Our work also
provides a baseline for applying POVMs during the qubits sequence. Integrated
with the process tensor, the measurement based operator can be very useful in
the real-time quantum error correction [61]. It is also interesting to explore
quantum dynamics with varying time steps, for example, to study the coherent-
to-incoherent transition of noise and the change of memory length. Since
determining a larger tensor network demands exponentially more resources, we
need useful tricks to compress the set of local operators needed. We can also
wisely choose the placement of operations to more efficiently probe systems of
interest.

###### Acknowledgements.

We thank Kavan Modi for helpful discussion. The work reported here was
supported by the National Key Research and Development Program of China (Grant
No. 2019YFA0308602, No. 2016YFA0301700), the National Natural Science
Foundation of China (Grants No. 12074336, No. 11934010, No. 11775129), the
Fundamental Research Funds for the Central Universities in China (No.
2020XZZX002-01), and the Anhui Initiative in Quantum Information Technologies
(Grant No. AHY080000). Y.Y. acknowledge the funding support from Tencent
Corporation. This work was partially conducted at the University of Science
and Technology of China Center for Micro- and Nanoscale Research and
Fabrication.

## References

  * [1] M. A. Nielsen and I. L. Chuang,  Quantum Computation and Quantum Information,  (Cambridge University Press, Cambridge, England 2010).
  * [2] H. P. Breuer and F. Petruccione,  The theory of open quantum systems,  (Oxford University Press on Demand, 2002).
  * [3] Á. Rivas and S. F. Huelga,  Open Quantum Systems,  (Springer Berlin Heidelberg, Berlin, 2012).
  * [4] A. Kossakowski,  On quantum statistical mechanics of non-Hamiltonian systems,  Reports Math. Phys. 3, 247 (1972).
  * [5] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan  Completely Positive Dynamical Semigroups of N-Level Systems,  J. Math. Phys. 17, 821 (1976).
  * [6] G. Lindblad,  On the generators of quantum dynamical semigroups,  Commun. Math. Phys. 48, 119–130 (1976).
  * [7] D. A. Lidar,  Lecture Notes on the Theory of Open Quantum Systems,  arXiv:1902.00967 (2019).
  * [8] Á. Rivas, S. F. Huelga, and M. B. Plenio,  Quantum Non-Markovianity: Characterization, Quantification and Detection,  Reports Prog. Phys. 77, 094001 (2014).
  * [9] F. A. Pollock, C. Rodríguez-Rosario, T. Frauenheim, M. Paternostro, and K. Modi,  Operational Markov Condition for Quantum Processes,  Phys. Rev. Lett. 120, 040405 (2018).
  * [10] P. J. J. O’Malley, J. Kelly, R. Barends, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. G. Fowler, I. C. Hoi, E. Jeffrey, A. Megrant, J. Mutus, C. Neill, C. Quintana, P. Roushan, D. Sank, A. Vainsencher, J. Wenner, T. C. White, A. N. Korotkov, A. N. Cleland, and J. M. Martinis,  Qubit Metrology of Ultralow Phase Noise Using Randomized Benchmarking,  Phys. Rev. Appl. 3, 044009 (2015).
  * [11] M. Ringbauer, C. J. Wood, K. Modi, A. Gilchrist, A. G. White, and A. Fedrizzi,  Characterizing Quantum Dynamics with Initial System-Environment Correlations,  Phys. Rev. Lett. 114, 090402 (2015).
  * [12] S. Yu, Y. T. Wang, Z. J. Ke, W. Liu, Y. Meng, Z. P. Li, W. H. Zhang, G. Chen, J. S. Tang, C. F. Li, and G. C. Guo,  Experimental Investigation of Spectra of Dynamical Maps and Their Relation to Non-Markovianity,  Phys. Rev. Lett. 120, 60406 (2018).
  * [13] Y.-N. Lu, Y.-R. Zhang, G.-Q. Liu, F. Nori, H. Fan, and X.-Y. Pan,  Observing Information Backflow from Controllable Non-Markovian Multichannels in Diamond,  Phys. Rev. Lett. 124, 210502 (2020).
  * [14] Y.-Q. Chen, K.-L. Ma, Y.-C. Zheng, J. Allcock, S. Zhang, and C.-Y. Hsieh,  Non-Markovian Noise Characterization with the Transfer Tensor Method,  Phys. Rev. Appl. 13, 034045 (2020).
  * [15] G. O. Samach, A. Greene, J. Borregaard, M. Christandl, D. K. Kim, C. M. McNally, A. Melville, B. M. Niedzielski, Y. Sung, D. Rosenberg, M. E. Schwartz, J. L. Yoder, T. P. Orlando, J. I.-J. Wang, S. Gustavsson, M. Kjaergaard, and W. D. Oliver,  Lindblad Tomography of a Superconducting Quantum Processor,  arXiv: 2105.02338 (2021).
  * [16] H. P. Breuer, E. M. Laine, and J. Piilo,  Measure for the Degree of Non-Markovian Behavior of Quantum Processes in Open Systems,  Phys. Rev. Lett. 103, 1 (2009).
  * [17] Á. Rivas, S. F. Huelga, and M. B. Plenio,  Entanglement and Non-Markovianity of Quantum Evolutions,  Phys. Rev. Lett. 105, 1 (2009).
  * [18] H. P. Breuer, E.-M. Laine, J. Piilo, and B. Vacchini,  Colloquium : Non-Markovian Dynamics in Open Quantum Systems,  Rev. Mod. Phys. 88, 021002 (2016).
  * [19] L. Li, M. J. W. Hall, and H. M. Wiseman,  Concepts of Quantum Non-Markovianity: A Hierarchy,  Phys. Rep. 759, 1 (2018).
  * [20] C.-F. Li, G.-C. Guo, and J. Piilo,  Non-Markovian Quantum Dynamics: What Does It Mean?,  EPL (Europhysics Lett. 127, 50001 (2019).
  * [21] C.-F. Li, G.-C. Guo, and J. Piilo,  Non-Markovian Quantum Dynamics: What Is It Good For?,  EPL (Europhysics Lett. 128, 30001 (2020).
  * [22] S. Milz, M. S. Kim, F. A. Pollock, and K. Modi,  Completely Positive Divisibility Does Not Mean Markovianity,  Phys. Rev. Lett. 123, 40401 (2019).
  * [23] F. Costa and S. Shrapnel,  Quantum Causal Modelling,  New J. Phys. 18, 063032 (2016).
  * [24] S. Milz, F. Sakuldee, F. A. Pollock, and K. Modi,  Kolmogorov Extension Theorem for (Quantum) Causal Modelling and General Probabilistic Theories,  Quantum 4, 255 (2020).
  * [25] F. A. Pollock, C. Rodríguez-Rosario, T. Frauenheim, M. Paternostro, and K. Modi,  Non-Markovian Quantum Processes: Complete Framework and Efficient Characterization,  Phys. Rev. A 97, 012127 (2018).
  * [26] S. Milz, F. A. Pollock, and K. Modi,  Reconstructing Non-Markovian Quantum Dynamics with Limited Control,  Phys. Rev. A 98, 012108 (2018).
  * [27] L. Xiang, Z. Zong, Z. Sun, Z. Zhan, Y. Fei, Z. Dong, C. Run, Z. Jia, P. Duan, J. Wu, Y. Yin, and G. Guo,  Simultaneous Feedback and Feedforward Control and Its Application to Realize a Random Walk on the Bloch Sphere in an Xmon-Superconducting-Qubit System,  Phys. Rev. Appl. 14, 014099 (2020).
  * [28] G. A. L. White, C. D. Hill, F. A. Pollock, L. C. L. Hollenberg, and K. Modi,  Demonstration of Non-Markovian Process Characterisation and Control on a Quantum Processor,  Nat. Commun. 11, 6301 (2020).
  * [29] R. Penrose,  Tensor Methods in Algebraic Geometry,  Ph.D. Thesis, Cambridge University (1956).
  * [30] R. Penrose,  Applications of negative dimensional tensors in Combinatorial Mathematics and its Applications, edited by D. Welsh,  (Academic Press, New York 1971).
  * [31] J. Biamonte and V. Bergholm,  Tensor Networks in a Nutshell,  arXiv:1708.00006 (2017).
  * [32] G. Chiribella, G. M. D’Ariano, and P. Perinotti,  Quantum Circuit Architecture,  Phys. Rev. Lett. 101, 1 (2008).
  * [33] U. Fano,  Description of States in Quantum Mechanics by Density Matrix and Operator Techniques,  Rev. Mod. Phys. 29, 74 (1957).
  * [34] A. Kuah, K. Modi, C. A. Rodruez-Rosario, and E. C. G. Sudarshan,  How State Preparation Can Affect a Quantum Experiment: Quantum Process Tomography for Open Systems,  Phys. Rev. A 76, 042113 (2007).
  * [35] V. Vedral,  The Role of Relative Entropy in Quantum Information Theory,  Rev. Mod. Phys. 74, 197 (2002).
  * [36] M.-D. Choi,  Positive Linear Maps on Complex,  Linear Algebra Appl. 10, 285 (1975).
  * [37] Geometry of quantum states: an introduction to quantum entanglement,  (Cambridge University Press, 2017).
  * [38] T. Wang, Z. Zhang, L. Xiang, Z. Jia, P. Duan, Z. Zong, Z. Sun, Z. Dong, J. Wu, Y. Yin, and G. Guo,  Experimental Realization of a Fast Controlled-Z Gate via a Shortcut to Adiabaticity,  Phys. Rev. Appl. 11, 034030 (2019).
  * [39] Z. Zhan, C. Run, Z. Zong, L. Xiang, Y. Fei, Z. Sun, Y. Wu, Z. Jia, P. Duan, J. Wu, Y. Yin, and G. Guo,  Experimental Determination of Electronic States via Digitized Shortcut-to-Adiabaticity,  arXiv:2103.06098 (2021).
  * [40] E. Jeffrey, D. Sank, J. Y. Mutus, T. C. White, J. Kelly, R. Barends, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, A. Megrant, P. J. J. O’Malley, C. Neill, P. Roushan, A. Vainsencher, J. Wenner, A. N. Cleland, and J. M. Martinis,  Fast Accurate State Measurement with Superconducting Qubits,  Phys. Rev. Lett. 112, 190504 (2014).
  * [41] J. E. Johnson, C. Macklin, D. H. Slichter, R. Vijay, E. B. Weingarten, J. Clarke, and I. Siddiqi,  Heralded State Preparation in a Superconducting Qubit,  Phys. Rev. Lett. 109, 050506 (2012).
  * [42] S. Milz, F. A. Pollock, and K. Modi,  An Introduction to Operational Quantum Dynamics,  Open Syst. Inf. Dyn. 24, 1740016 (2017).
  * [43] T. Yamamoto, M. Neeley, E. Lucero, R. C. Bialczak, J. Kelly, M. Lenander, M. Mariantoni, A. D. O’Connell, D. Sank, H. Wang, M. Weides, J. Wenner, Y. Yin, A. N. Cleland, and J. M. Martinis,  Quantum Process Tomography of Two-Qubit Controlled-Z and Controlled-NOT Gates Using Superconducting Phase Qubits,  Phys. Rev. B 82, 184515 (2010).
  * [44] M. Steffen, M. Ansmann, R. Mcdermott, N. Katz, R. C. Bialczak, E. Lucero, M. Neeley, E. M. Weig, A. N. Cleland, and J. M. Martinis,  State Tomography of Capacitively Shunted Phase Qubits with High Fidelity,  Phys. Rev. Lett. 97, 050502 (2006).
  * [45] S. T. Flammia, A. Silberfarb, and C. M. Caves,  Minimal Informationally Complete Measurements for Pure States,  Found. Phys. 35, 1985 (2005).
  * [46] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant,  Array Programming with NumPy,  Nature 585, 357 (2020).
  * [47] R. Jozsa and B. Schumacher,  A New Proof of the Quantum Noiseless Coding Theorem,  J. Mod. Opt. 41, 2343 (1994).
  * [48] R. Barends, J. Kelly, A. Megrant, A. Veitia, D. Sank, E. Jeffrey, T. C. White, J. Mutus, A. G. Fowler, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, C. Neill, P. O’Malley, P. Roushan, A. Vainsencher, J. Wenner, A. N. Korotkov, A. N. Cleland, and J. M. Martinis,  Superconducting Quantum Circuits at the Surface Code Threshold for Fault Tolerance,  Nature 508, 500 (2014).
  * [49] D. I. Schuster, A. Wallraff, A. Blais, L. Frunzio, R.-S. Huang, J. Majer, S. M. Girvin, and R. J. Schoelkopf,  Ac Stark Shift and Dephasing of a Superconducting Qubit Strongly Coupled to a Cavity Field,  Phys. Rev. Lett. 94, 123602 (2005).
  * [50] P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando, S. Gustavsson, and W. D. Oliver,  A Quantum Engineer’s Guide to Superconducting Qubits,  Appl. Phys. Rev. 6, 021318 (2019).
  * [51] D. Sank, Z. Chen, M. Khezri, J. Kelly, R. Barends, B. Campbell, Y. Chen, B. Chiaro, A. Dunsworth, A. Fowler, E. Jeffrey, E. Lucero, A. Megrant, J. Mutus, M. Neeley, C. Neill, P. J. J. O’Malley, C. Quintana, P. Roushan, A. Vainsencher, T. White, J. Wenner, A. N. Korotkov, and J. M. Martinis,  Measurement-Induced State Transitions in a Superconducting Qubit: Beyond the Rotating Wave Approximation,  Phys. Rev. Lett. 117, 190503 (2016).
  * [52] K. Modi,  Operational Approach to Open Dynamics and Quantifying Initial Correlations,  Sci. Rep. 2, 581 (2012).
  * [53] L. Xiang,  data analysis code as git repository: QUCSE,  http://github.com/xlelephant/qucse.
  * [54] J. M. Martinis, S. Nam, J. Aumentado, K. M. Lang, and C. Urbina,  Decoherence of a Superconducting Qubit Due to Bias Noise,  Phys. Rev. B 67, 094510 (2003).
  * [55] J. M. Martinis, K. B. Cooper, R. McDermott, M. Steffen, M. Ansmann, K. D. Osborn, K. Cicak, S. Oh, D. P. Pappas, R. W. Simmonds, and C. C. Yu,  Decoherence in Josephson Qubits from Dielectric Loss,  Phys. Rev. Lett. 95, 210503 (2005).
  * [56] R. C. Bialczak, R. McDermott, M. Ansmann, M. Hofheinz, N. Katz, E. Lucero, M. Neeley, A. D. O’Connell, H. Wang, A. N. Cleland, and J. M. Martinis,  1/f Flux Noise in Josephson Phase Qubits,  Phys. Rev. Lett. 99, 187006 (2007).
  * [57] D. Sank, R. Barends, R. C. Bialczak, Y. Chen, J. Kelly, M. Lenander, E. Lucero, M. Mariantoni, A. Megrant, M. Neeley, P. J. J. O’Malley, A. Vainsencher, H. Wang, J. Wenner, T. C. White, T. Yamamoto, Y. Yin, A. N. Cleland, and J. M. Martinis,  Flux Noise Probed with Real Time Qubit Tomography in a Josephson Phase Qubit,  Phys. Rev. Lett. 109, 067001 (2012).
  * [58] A. P. Sears, A. Petrenko, G. Catelani, L. Sun, H. Paik, G. Kirchmair, L. Frunzio, L. I. Glazman, S. M. Girvin, and R. J. Schoelkopf,  Photon Shot Noise Dephasing in the Strong-Dispersive Limit of Circuit QED,  Phys. Rev. B 86, 180504 (2012).
  * [59] D. Ristè, C. C. Bultink, M. J. Tiggelman, R. N. Schouten, K. W. Lehnert, and L. DiCarlo,  Millisecond Charge-Parity Fluctuations and Induced Decoherence in a Superconducting Transmon Qubit,  Nat. Commun. 4, 1913 (2013).
  * [60] M. McEwen, D. Kafri, Z. Chen, J. Atalaya, K. J. Satzinger, C. Quintana, P. V. Klimov, D. Sank, C. Gidney, A. G. Fowler, F. Arute, K. Arya, B. Buckley, B. Burkett, N. Bushnell, B. Chiaro, R. Collins, S. Demura, A. Dunsworth, C. Erickson, B. Foxen, M. Giustina, T. Huang, S. Hong, E. Jeffrey, S. Kim, K. Kechedzhi, F. Kostritsa, P. Laptev, A. Megrant, X. Mi, J. Mutus, O. Naaman, M. Neeley, C. Neill, M. Niu, A. Paler, N. Redd, P. Roushan, T. C. White, J. Yao, P. Yeh, A. Zalcman, Y. Chen, V. N. Smelyanskiy, J. M. Martinis, H. Neven, J. Kelly, A. N. Korotkov, A. G. Petukhov, and R. Barends,  Removing Leakage-Induced Correlated Errors in Superconducting Quantum Error Correction,  Nat. Commun. 12, 1761 (2021).
  * [61] L. Hu, Y. Ma, W. Cai, X. Mu, Y. Xu, W. Wang, Y. Wu, H. Wang, Y. P. Song, C. L. Zou, S. M. Girvin, L. M. Duan, and L. Sun,  Quantum Error Correction and Universal Gate Set Operation on a Binomial Bosonic Logical Qubit,  Nat. Phys. 15, 503 (2019).
  * [62] A. Strathearn, P. Kirton, D. Kilda, J. Keeling, and B. W. Lovett,  Efficient Non-Markovian Quantum Dynamics Using Time-Evolving Matrix Product Operators,  Nat. Commun. 9, 3322 (2018).


# Gravitational waves from axions annihilation through scattering amplitude

Jing Yang  Fa Peng Huang<EMAIL_ADDRESS>MOE Key Laboratory of
TianQin Mission, TianQin Research Center for Gravitational Physics & School of
Physics and Astronomy, Frontiers Science Center for TianQin, Gravitational
Wave Research Center of CNSA, Sun Yat-sen University (Zhuhai Campus), Zhuhai
519082, China

###### Abstract

We use the scattering method to calculate the gravitational wave from axions
annihilation in the axion cloud formed by the superradiance process around the
Kerr black hole. We consider axions annihilating to gravitons as a three-body
decay process and then calculate the corresponding decay width. In this
approach, we can simply obtain the radiation power of gravitational wave and
give the analytical approximate result with the spin effects of the Kerr black
hole. Our study can also provide a cross-check to the numerical results in the
traditional method.

##  I introduction

Motivated by the current situation of dark matter search, recently more and
more studies focus on ultralight bosons, such as ultralight aixons or dark
photons. Axion Peccei:1977hh  or axion-like particle Svrcek:2006yi  is one of
the most promising dark matter candidates Preskill:1982cy ; Sikivie:2006ni .
As pseudo scalar particles, they can form Bose-Einstein condensation and
further form dense axion object, like axion star or axion cluster. Especially,
the so-called axion cloud can be formed from the superradiance process around
Kerr black hole (BH) Penrose:1969pc ; zeld ; Starobinsky:1973aij ;
Detweiler:1980uk . Once the axion cloud is formed, abundant signals could be
produced. For example, the axion cloud can be a promising gravitational wave
(GW) source Brito:2015oca ; Arvanitaki:2014wva ; Arvanitaki:2016qwi ;
Baryakhtar:2017ngi ; Dev:2016hxv ; Xie:2022uvp , which can be used to explore
the properties of axion particles Brito:2017zvb . To unravel the axion
properties from GW signals, it is crucial to obtain a reliable prediction of
the GW radiation power from the axion cloud. Ref. Arvanitaki:2010sy  studied
the GW radiation power in the flat space approximation. Ref. Brito:2014wla
calculated the GW radiation power of axions annihilating to gravitons in the
Schwarzschild metric. Refs. Brito:2017zvb Yoshino:2013ofa  calculated the GW
radiation power in the Kerr metric. However, the traditional method is
complicated in the Kerr metric. To simplify the calculation and cross-check
this important result, we use the method of scattering amplitude to calculate
power radiation generated by the axions annihilation process.

In Sec. II, we simply review the basic concepts about axion cloud from
superradiance, then we give the relevant Feynman rules for axions annihilating
to gravitons and derive the formulations about decay width. The radiation
power of gravitational waves from axions annihilating to gravitons is
calculated in Sec. III where we give both analytical and numerical results. In
Sec. IV, we derive the formula used in the flat space approximation through
the scattering amplitude method in the classical source limit. Concise
conclusions are given in Sec. V.

##  II Scattering amplitude of axions annihilating to gravitons around Kerr
black hole

###  II.1 Axion cloud from superradiance around Kerr black hole

When the Compton wavelength of the axion boson becomes comparable with the
Kerr BH size, axion cloud could form by extracting the energy and angular
momentum of the Kerr BH through the so-called superradiance process. In
analogy with the atom, the axion cloud and the Kerr BH can form the
“gravitational atom”. This formation of axion cloud can be understood from the
particle or wave aspect. From the perspective of particle, the axion cloud can
be formed due to the Penrose process Penrose:1969pc  of Kerr BH since there
exists an ergosphere between the out horizon and out static limit surface.
From the perspective of wave, we could directly solve the Klein-Gordon
equation in the Kerr metric. At large distance, the radial solution is just
confluent hypergeometric function. While, at small distance (close to the
outer event horizon), the radial solution can be expressed as the
hypergeometric function. Matching them together, one can find the exponential
growth solution when the Compton length is larger than the BH size
Detweiler:1980uk . Once the axion cloud or gravitational atom is formed, they
can emit GW signals from axions annihilating to gravitons or energy-level
transitions of gravitational atom. In this work, we focus on the GW signals
from axions annihilating to gravitons.

###  II.2 Relevant Feynman rules for axions annihilating to gravitons

The Einstein-Hilbert action in the general relativity is

${\cal{S}}=\int d^{4}x{\cal{L}}=\int
d^{4}x\sqrt{-g}(\frac{2R}{\kappa^{2}}+{\cal L}_{\text{matter}}),$ (1)

where the gravitational coupling $\kappa^{2}\equiv 32\pi G=32\pi/M_{pl}^{2}$
and Riemann tensor
$R^{\mu}_{\nu\alpha\beta}\equiv\partial_{\alpha}\Gamma^{\mu}_{\nu\beta}-\partial_{\beta}\Gamma^{\mu}_{\nu\alpha}+\Gamma^{\mu}_{\sigma\alpha}\Gamma^{\sigma}_{\nu\beta}-\Gamma^{\mu}_{\sigma\beta}\Gamma^{\sigma}_{\nu\alpha}$.
We expand the metric as

$g_{\mu\nu}=\eta_{\mu\nu}+\kappa h_{\mu\nu},$ (2)

where $\eta_{\mu\nu}=\text{diag}\\{1,-1,-1,-1\\}$ is the flat space-time
metric and $h_{\mu\nu}$ is the graviton field. Then we have

$\displaystyle g^{\mu\nu}$ $\displaystyle=$ $\displaystyle\eta^{\mu\nu}-\kappa
h^{\mu\nu}+\kappa^{2}h^{\mu}_{\alpha}h^{\alpha\nu}+...,$
$\displaystyle\sqrt{-g}$ $\displaystyle=$ $\displaystyle 1+\frac{1}{2}\kappa
h-\frac{1}{4}\kappa^{2}h^{\alpha}_{\beta}h^{\beta}_{\alpha}+\frac{1}{8}\kappa^{2}h^{2}+...,$
(3)

where $h=\eta^{\mu\nu}h_{\mu\nu}$ and we use $\eta^{\mu\nu}$ to raise and
lower indices. The relevant Lagrangian of axion field is ${\cal
L}=\frac{1}{2}g_{\mu\nu}\partial^{\mu}\Phi\partial^{\nu}\Phi-\frac{1}{2}m_{a}^{2}\Phi^{2}$,
where $m_{a}$ is the axion mass. After expanding the Einstein-Hilbert action
with $h_{\mu\nu}$, one can get the following Feynman rules for axions
annihilating to gravitons in the harmonic gauge DeWitt:1967uc Bjerrum-
Bohr:2004qcf .

$\displaystyle=\frac{i}{q^{2}-m_{a}^{2}+i\epsilon}$

$\displaystyle=\frac{i{\cal P}^{\alpha\beta\gamma\delta}}{q^{2}+i\epsilon}$

where

${\cal
P}^{\alpha\beta\gamma\delta}=\frac{1}{2}\left(\eta^{\alpha\gamma}\eta^{\beta\delta}+\eta^{\beta\gamma}\eta^{\alpha\delta}-\eta^{\alpha\beta}\eta^{\gamma\delta}\right).$
(4)

$\displaystyle=\tau_{1}^{\mu\nu}(p,p^{\prime},m_{a})$

where

$\tau_{1}^{\mu\nu}(p,p^{\prime},m_{a})=-\frac{i\kappa}{2}\Big{(}p^{\mu}p^{\prime\nu}+p^{\nu}p^{\prime\mu}-\eta^{\mu\nu}\left((p\cdot
p^{\prime})-m_{a}^{2}\right)\Big{)}.$ (5)

$\displaystyle=\tau_{2}^{\eta\lambda\rho\sigma}(p,p^{\prime},m_{a})$

where

$\displaystyle\tau_{2}^{\eta\lambda\rho\sigma}(p,p^{\prime})$
$\displaystyle={i\kappa^{2}}\bigg{(}\Big{(}I^{\eta\lambda\alpha\delta}{I}^{\rho\sigma\beta}_{\
\ \ \
\delta}-\frac{1}{4}\left(\eta^{\eta\lambda}I^{\rho\sigma\alpha\beta}+\eta^{\rho\sigma}I^{\eta\lambda\alpha\beta}\right)\Big{)}(p_{\alpha}p^{\prime}_{\beta}+p^{\prime}_{\alpha}p_{\beta})$
(6)
$\displaystyle-\frac{1}{2}\Big{(}I^{\eta\lambda\rho\sigma}-\frac{1}{2}\eta^{\eta\lambda}\eta^{\rho\sigma}\Big{)}\Big{(}(p\cdot
p^{\prime})-m_{a}^{2}\Big{)}\bigg{)},$

with

$\displaystyle{\cal I}_{\mu\nu\alpha\beta}$ $\displaystyle=$
$\displaystyle{1\over
2}(\eta_{\mu\alpha}\eta_{\nu\beta}+\eta_{\nu\alpha}\eta_{\mu\beta}).$ (7)

$\displaystyle={\tau_{3}}_{\alpha\beta\gamma\delta}^{\mu\nu}(k,q)$

where

$\displaystyle{\tau_{3}}_{\alpha\beta\gamma\delta}^{\mu\nu}(k,q)$
$\displaystyle=-\frac{i\kappa}{2}\times\bigg{(}{\cal
P}_{\alpha\beta\gamma\delta}\Big{(}k^{\mu}k^{\nu}+(k-q)^{\mu}(k-q)^{\nu}+q^{\mu}q^{\nu}-\frac{3}{2}\eta^{\mu\nu}q^{2}\Big{)}$
(8) $\displaystyle+2q_{\lambda}q_{\sigma}\Big{(}I_{\alpha\beta}^{\ \ \
\sigma\lambda}I_{\gamma\delta}^{\ \ \ \mu\nu}+I_{\gamma\delta}^{\ \ \
\sigma\lambda}I_{\alpha\beta}^{\ \ \ \mu\nu}-I_{\alpha\beta}^{\ \ \
\mu\sigma}I_{\gamma\delta}^{\ \ \ \nu\lambda}-I_{\gamma\delta}^{\ \ \
\mu\sigma}I_{\alpha\beta}^{\ \ \ \nu\lambda}\Big{)}$
$\displaystyle+\Big{(}q_{\lambda}q^{\mu}\big{(}\eta_{\alpha\beta}I_{\gamma\delta}^{\
\ \ \nu\lambda}+\eta_{\gamma\delta}I_{\alpha\beta}^{\ \ \
\nu\lambda}\big{)}+q_{\lambda}q^{\nu}\big{(}\eta_{\alpha\beta}I_{\gamma\delta}^{\
\ \ \mu\lambda}+\eta_{\gamma\delta}I_{\alpha\beta}^{\ \ \ \mu\lambda}\big{)}$
$\displaystyle-q^{2}\left(\eta_{\alpha\beta}I_{\gamma\delta}^{\ \ \
\mu\nu}-\eta_{\gamma\delta}I_{\alpha\beta}^{\ \ \
\mu\nu}\right)-\eta^{\mu\nu}q_{\sigma}q_{\lambda}\left(\eta_{\alpha\beta}I_{\gamma\delta}^{\
\ \ \sigma\lambda}+\eta_{\gamma\delta}I_{\alpha\beta}^{\ \ \
\sigma\lambda}\right)\Big{)}$
$\displaystyle+\Big{(}2q_{\lambda}\big{(}I_{\alpha\beta}^{\ \ \
\lambda\sigma}I_{\gamma\delta\sigma}^{\ \ \ \
\nu}(k-q)^{\mu}+I_{\alpha\beta}^{\ \ \ \lambda\sigma}I_{\gamma\delta\sigma}^{\
\ \ \ \mu}(k-q)^{\nu}-I_{\gamma\delta}^{\ \ \
\lambda\sigma}I_{\alpha\beta\sigma}^{\ \ \ \ \nu}k^{\mu}-I_{\gamma\delta}^{\ \
\ \lambda\sigma}I_{\alpha\beta\sigma}^{\ \ \ \ \mu}k^{\nu}\big{)}$
$\displaystyle+q^{2}\left(I_{\alpha\beta\sigma}^{\ \ \ \
\mu}I_{\gamma\delta}^{\ \ \ \nu\sigma}+I_{\alpha\beta}^{\ \ \
\nu\sigma}I_{\gamma\delta\sigma}^{\ \ \ \
\mu}\right)+\eta^{\mu\nu}q_{\sigma}q_{\lambda}\big{(}I_{\alpha\beta}^{\ \ \
\lambda\rho}I_{\gamma\delta\rho}^{\ \ \ \ \sigma}+I_{\gamma\delta}^{\ \ \
\lambda\rho}I_{\alpha\beta\rho}^{\ \ \ \ \sigma}\big{)}\Big{)}$
$\displaystyle+\Big{(}(k^{2}+(k-q)^{2})\big{(}I_{\alpha\beta}^{\ \ \
\mu\sigma}I_{\gamma\delta\sigma}^{\ \ \ \ \nu}+I_{\gamma\delta}^{\ \ \
\mu\sigma}I_{\alpha\beta\sigma}^{\ \ \ \ \nu}-\frac{1}{2}\eta^{\mu\nu}{\cal
P}_{\alpha\beta\gamma\delta}\big{)}$ $\displaystyle-\left(I_{\gamma\delta}^{\
\ \ \mu\nu}\eta_{\alpha\beta}k^{2}+I_{\alpha\beta}^{\ \ \
\mu\nu}\eta_{\gamma\delta}(k-q)^{2}\right)\Big{)}\bigg{)}.$

###  II.3 Formulations about decay width

The annihilation process of two axions in the background field can be
understood as a decay process of a three-body bound state which is composed of
two axions and a BH. If the “fine structure constant” $\alpha=GM_{b}m_{a}$ of
the gravitational system is small where $M_{b}$ is BH mass, then the internal
motions are slow and we can treat the bound state as a non-relativistic
system. Thus we can approximately expand the bound state which is denoted as
$|B\rangle$ in states of the free field theory in the center-of-mass (COM)
frame as

$|B\rangle=\sqrt{2M}\iint\frac{d^{3}\boldsymbol{p_{1}}}{(2\pi)^{3}}\frac{d^{3}\boldsymbol{p_{2}}}{(2\pi)^{3}}\Psi_{1}(\boldsymbol{p_{1}})\Psi_{2}(\boldsymbol{p_{2}})\frac{1}{\sqrt{2E_{1}}\sqrt{2E_{2}}\sqrt{2E_{b}}}|\boldsymbol{p_{b}},\boldsymbol{p_{1}},\boldsymbol{p_{2}}\rangle,$
(9)

where $\boldsymbol{p_{b}}$ is the momentum of BH, $\boldsymbol{p_{1}}$,
$\boldsymbol{p_{2}}$ are the momenta of two axions satisfying
$\boldsymbol{p_{b}}=-\boldsymbol{p_{1}}-\boldsymbol{p_{2}}$ in the COM frame,
$\Psi_{1}$($\Psi_{2}$) is the wave function of bounded axion in momentum
space, and $M$ is the mass of the bound state. $E_{1}$, $E_{2}$, $E_{b}$ are
the energy of axions and BH respectively.

The decay with of the bound state in the COM frame is

$\displaystyle\Gamma$ $\displaystyle\simeq$
$\displaystyle\iint\frac{d^{3}\boldsymbol{k}}{(2\pi)^{3}}\frac{d^{3}\boldsymbol{k_{b}}}{(2\pi)^{3}}\frac{1}{2|\boldsymbol{k}|}\frac{1}{2E_{b}^{\prime}}\big{|}\iint\frac{d^{3}\boldsymbol{p_{1}}}{(2\pi)^{3}}\frac{d^{3}\boldsymbol{p_{2}}}{(2\pi)^{3}}\Psi_{1}(\boldsymbol{p_{1}})\Psi_{2}(\boldsymbol{p_{2}})\frac{1}{\sqrt{2E_{1}}\sqrt{2E_{2}}\sqrt{2E_{b}}}$
(10) $\displaystyle{\cal
M}(\boldsymbol{p_{b}},\boldsymbol{p_{1}},\boldsymbol{p_{2}}\rightarrow\boldsymbol{k},\boldsymbol{k_{b}})\big{|}^{2}(2\pi)^{4}\delta^{4}(p_{1}+p_{2}+p_{b}-k-k_{b}),$

where ${\cal
M}(\boldsymbol{p_{b}},\boldsymbol{p_{1}},\boldsymbol{p_{2}}\rightarrow\boldsymbol{k},\boldsymbol{k_{b}})$
is the scattering amplitude of
$|\boldsymbol{p_{b}},\boldsymbol{p_{1}},\boldsymbol{p_{2}}\rangle\rightarrow|\boldsymbol{k},\boldsymbol{k_{b}}\rangle$
process, $\boldsymbol{k}$ and $\boldsymbol{k_{b}}$ are the momenta of graviton
and BH respectively in the final state, $E_{b}^{\prime}$ is the energy of BH
in the final state. There are four Feynman diagrams relevant to the
annihilation process $a+a+b\rightarrow g+b$ at tree level as shown in Fig. 1
where $b$ and $g$ represent the BH and graviton respectively.

Figure 1: The Feynman diagrams for the scattering amplitude of
$|\boldsymbol{p_{b}},\boldsymbol{p_{1}},\boldsymbol{p_{2}}\rangle\rightarrow|\boldsymbol{k},\boldsymbol{k_{b}}\rangle$
process. The blue and black dashed lines depict the axion and BH respectively
while the red double-wavy line represents the graviton.

We can approximately take the COM frame as the BH frame due to $M_{b}\gg
m_{a}$, which means we can take the static source limit

$E_{b}\simeq E_{b}^{\prime}\simeq M_{b},\ \ \ \boldsymbol{p_{b}}=0.$ (11)

Then we have

$\displaystyle\Gamma$ $\displaystyle\simeq$
$\displaystyle\int\frac{d^{3}\boldsymbol{k}}{(2\pi)^{3}}\frac{1}{2|\boldsymbol{k}|}\frac{1}{2M_{b}}\big{|}\iint\frac{d^{3}\boldsymbol{p_{1}}}{(2\pi)^{3}}\frac{d^{3}\boldsymbol{p_{2}}}{(2\pi)^{3}}\Psi_{1}(\boldsymbol{p_{1}})\Psi_{2}(\boldsymbol{p_{2}})\frac{1}{\sqrt{2E_{1}}\sqrt{2E_{2}}\sqrt{2M_{b}}}$
(12) $\displaystyle{\cal
M}(\boldsymbol{p_{b}}=0,\boldsymbol{p_{1}},\boldsymbol{p_{2}}\rightarrow\boldsymbol{k},\boldsymbol{k_{b}}=-\boldsymbol{k})\big{|}^{2}2\pi\delta(E_{1}+E_{2}-|\boldsymbol{k}|),$

where ${\cal
M}(\boldsymbol{p_{b}}=0,~{}\boldsymbol{p_{1}},\boldsymbol{p_{2}}\rightarrow\boldsymbol{k},\boldsymbol{k_{b}}=-\boldsymbol{k})$
is just the scattering amplitude of two axions annihilating to a graviton in
the static classical source background.

Figure 2: The Feynman diagram for the scattering amplitude of
$|\boldsymbol{p_{b}},\boldsymbol{p_{1}},\boldsymbol{p_{2}}\rangle\rightarrow|\boldsymbol{k},\boldsymbol{k_{b}}\rangle$
process (left panel). We can take the classical static source limit due to
$M_{b}\gg m_{a}$ which means the gravitational field generated by a BH is
considered as a background field (right panel).

We can decompose the scattering amplitude of
$\langle\boldsymbol{k},\boldsymbol{k_{b}}|\boldsymbol{p_{b}},\boldsymbol{p_{1}},\boldsymbol{p_{2}}\rangle$
as

$\displaystyle i{\cal
M}(\boldsymbol{p_{b}},\boldsymbol{p_{1}},\boldsymbol{p_{2}}\rightarrow\boldsymbol{k},\boldsymbol{k_{b}})=$

$\displaystyle\times$

$\ \ \ =i{\cal
M}_{1}^{\mu\nu}(\boldsymbol{p_{1}},\boldsymbol{p_{2}},\boldsymbol{k})\times{\cal
M}_{2\mu\nu}(\boldsymbol{p_{b}},\boldsymbol{k_{b}}),\ \ \ \ \ \ $ (13)

where

$\displaystyle(2\pi)^{4}\delta^{4}(k_{b}-p_{b}-q){\cal M}_{2\mu\nu}$
$\displaystyle=$ $\displaystyle\frac{\kappa{\cal P}_{\mu\nu}^{\alpha\beta}\int
d^{4}xe^{-iq\cdotp x}\langle\boldsymbol{k_{b}}|{\cal
T}_{b\alpha\beta}|\boldsymbol{p_{b}}\rangle}{2(k_{b}-p_{b})^{2}}$ (14)
$\displaystyle=$ $\displaystyle\frac{\kappa\int d^{4}xe^{-iq\cdotp
x}\langle\boldsymbol{k_{b}}|{\cal T}_{b\mu\nu}-\eta_{\mu\nu}{\cal
T}_{b}/2|\boldsymbol{p_{b}}\rangle}{2(k_{b}-p_{b})^{2}}.$

Here ${\cal T}_{b\mu\nu}$ is energy-momentum tensor of the BH’s field and
${\cal T}_{b}={\cal T}^{\mu}_{b\mu}$. We define ${\cal T}_{b\mu\nu}={\cal
T}^{\prime}_{b\mu\nu}e^{i(k_{b}-p_{b})\cdotp x}$ because
$\langle\boldsymbol{k_{b}}|{\cal T}_{b\alpha\beta}|\boldsymbol{p_{b}}\rangle$
has a factor $e^{i(k_{b}-p_{b})\cdotp x}$. Then we find

${\cal M}_{2\mu\nu}=\frac{\kappa\langle\boldsymbol{k_{b}}|{\cal
T}^{\prime}_{b\mu\nu}-\eta_{\mu\nu}{\cal
T}^{\prime}_{b}/2|\boldsymbol{p_{b}}\rangle}{2(k_{b}-p_{b})^{2}}.$ (15)

While we need to degenerate the BH from a quantum of quantum field into a
classical object and the finite volume effect can be considered as that ${\cal
T}_{b\mu\nu}$ has continuous spectra in coordinate space,

${\cal T}_{b\mu\nu}=\int d^{3}\boldsymbol{y}\tilde{{\cal
T}}_{b\mu\nu}(\boldsymbol{y})e^{i(k_{b}-p_{b})\cdotp
x}e^{-i(\boldsymbol{k_{b}-p_{b}})\cdotp\boldsymbol{y}},$ (16)

where $\boldsymbol{y}$ is the position relative to $\boldsymbol{x}$. Thus the
relation between energy-momentum tensor of BH $T_{b\mu\nu}$ and $\tilde{{\cal
T}}_{b\mu\nu}$ is

$\langle\boldsymbol{p}_{b}|\tilde{{\cal
T}}_{b\mu\nu}(\boldsymbol{y})|\boldsymbol{p}_{b}\rangle\xrightarrow{\text{classical
static limit}}2M_{b}T_{b\mu\nu}(\boldsymbol{y}),$ (17)

where the factor $2M_{b}$ comes from the normalization of the state
$\langle\boldsymbol{k}|\boldsymbol{p}\rangle=2E_{\boldsymbol{p}}(2\pi)^{3}\delta^{3}(\boldsymbol{p}-\boldsymbol{k})$
with $E_{\boldsymbol{p}}=\sqrt{\boldsymbol{p}^{2}+M_{b}^{2}}$ in the quantum
field theory. In the static limit, $|\boldsymbol{k_{b}}-\boldsymbol{p_{b}}|\ll
M_{b}$, we have

$\langle\boldsymbol{k}_{b}|\tilde{{\cal
T}}_{b\mu\nu}(\boldsymbol{y})|\boldsymbol{p}_{b}\rangle\approx\langle\boldsymbol{p}_{b}|\tilde{{\cal
T}}_{b\mu\nu}(\boldsymbol{y})|\boldsymbol{p}_{b}\rangle.$ (18)

Thus

$\displaystyle{\cal M}_{2\mu\nu}$ $\displaystyle=$ $\displaystyle\frac{-\kappa
M_{b}\int
d^{3}\boldsymbol{y}e^{-i(\boldsymbol{k_{b}-p_{b}})\cdotp\boldsymbol{y}}(T_{b\mu\nu}(\boldsymbol{y})-\eta_{\mu\nu}T(\boldsymbol{y})/2)}{(\boldsymbol{k}_{b}-\boldsymbol{p}_{b})^{2}}$
(19) $\displaystyle=$ $\displaystyle\frac{-\kappa
M_{b}(T_{b\mu\nu}(\boldsymbol{q})-\eta_{\mu\nu}T(\boldsymbol{q})/2)}{\boldsymbol{q}^{2}}.$

Note that the Einstein equation in the weak-field approximation is

$\Box h_{\mu\nu}=-\frac{\kappa}{2}(T_{\mu\nu}-\eta_{\mu\nu}T/2),$ (20)

if the source is static,we have

$\nabla^{2}h_{\mu\nu}=\frac{\kappa}{2}(T_{\mu\nu}-\eta_{\mu\nu}T/2),$ (21)

then the solution reads

$h_{\mu\nu}(\boldsymbol{x})=-\frac{\kappa}{8\pi}\int
d^{3}\boldsymbol{y}\frac{(T_{\mu\nu}(\boldsymbol{y})-\eta_{\mu\nu}T(\boldsymbol{y})/2)}{|\boldsymbol{x}-\boldsymbol{y}|}.$
(22)

The Fourier transformation of $h_{\mu\nu}(\boldsymbol{x})$ in momentum space
is

$\tilde{h}_{\mu\nu}(\boldsymbol{q})=-\frac{\kappa}{2}\frac{T_{b\mu\nu}(\boldsymbol{q})-\eta_{\mu\nu}T(\boldsymbol{q})/2}{\boldsymbol{q}^{2}}.$
(23)

Substituting Eq. 23 into Eq. 19, we get

${\cal M}_{2\mu\nu}=2M_{b}\tilde{h}_{\mu\nu}(\boldsymbol{q}),$ (24)

where $\tilde{h}_{\mu\nu}(\boldsymbol{q})$ is the field generated by the Kerr
BH in momentum space. Thus the decay width $\Gamma$ can be simplified to

$\Gamma=\int\frac{d^{3}\boldsymbol{k}}{(2\pi)^{3}}\frac{\pi\delta(E_{1}+E_{2}-|\boldsymbol{k}|)}{|\boldsymbol{k}|}\big{|}\iint\frac{d^{3}\boldsymbol{p_{1}}}{(2\pi)^{3}}\frac{d^{3}\boldsymbol{p_{2}}}{(2\pi)^{3}}\frac{\Psi_{1}(\boldsymbol{p_{1}})\Psi_{2}(\boldsymbol{p_{2}})}{\sqrt{2E_{1}}\sqrt{2E_{2}}}{\cal
M}_{1}^{\mu\nu}(\boldsymbol{p_{1}},\boldsymbol{p_{2}},\boldsymbol{k})\tilde{h}_{\mu\nu}(\boldsymbol{q})\big{|}^{2}.$
(25)

The gravitational background field of the Kerr BH in the weak-field
approximation is

$h_{\mu\nu}(\boldsymbol{x})=h_{\mu\nu,0}(\boldsymbol{x})+h_{\mu\nu,J}(\boldsymbol{x}),$
(26) $h_{\mu\nu,0}(\boldsymbol{x})=N_{\mu\nu}\frac{\kappa M_{b}}{16\pi r},$
(27) $h_{\mu\nu,J}(\boldsymbol{x})=-N^{i}_{\mu\nu}\frac{\kappa}{16\pi
r^{3}}\epsilon^{ikm}x^{k}J^{m},$ (28)

where

$N_{\mu\nu}=\eta_{\mu\nu}-2\eta_{\mu 0}\eta_{\nu 0},$ (29)
$N^{i}_{\mu\nu}=\eta^{0}_{\mu}\eta^{i}_{\nu}+\eta^{0}_{\nu}\eta^{i}_{\mu},\
i=1,2,3.$ (30)

Then the gravitational background field of the Kerr BH in the momentum space
$\tilde{h}_{\mu\nu}(\boldsymbol{q})$ can be derived as (Note that the
direction of $\boldsymbol{q}$ is incoming to the Kerr BH in the Feynman
diagrams which is opposite to the usual Fourier transform convention where the
direction of momentum is outgoing.)

$\tilde{h}_{\mu\nu}(\boldsymbol{q})=\tilde{h}_{\mu\nu,0}(\boldsymbol{q})+\tilde{h}_{\mu\nu,J}(\boldsymbol{q}),$
(31) $\tilde{h}_{\mu\nu,0}(\boldsymbol{q})=N_{\mu\nu}\frac{\kappa
M_{b}}{4\boldsymbol{q}^{2}},$ (32)
$\tilde{h}_{\mu\nu,J}(\boldsymbol{q})=N^{i}_{\mu\nu}\frac{-i\kappa\epsilon^{ikm}q^{k}J^{m}}{4\boldsymbol{q}^{2}}.$
(33)

And ${\cal M}_{1}^{\mu\nu}$ is the sum of the corresponding four Feynman
diagrams in Fig. 1,

${\cal M}_{1,a}^{s,\mu\nu}=\tau_{1}^{\rho\sigma}(p_{1},-p_{2})\frac{i{\cal
P}^{\gamma\delta}_{\rho\sigma}}{(p_{1}+p_{2})^{2}}\tau^{\mu\nu}_{3,\alpha\beta\gamma\delta}(-k,q)\epsilon^{s,\alpha\beta}\
\ \text{for Fig.1(a)},$ (34) ${\cal
M}_{1,b}^{s,\mu\nu}=\tau_{2,\alpha\beta}^{\mu\nu}(p_{1},-p_{2})\epsilon^{s,\alpha\beta}\
\ \text{for Fig.1(b)},$ (35) ${\cal
M}_{1,c}^{s,\mu\nu}=\tau_{1,\alpha\beta}(p_{1},p_{1}-k)\frac{i}{(p_{1}-k)^{2}-m_{a}^{2}}\tau^{\mu\nu}_{1}(p_{1}-k,-p_{2})\epsilon^{s,\alpha\beta}\
\ \text{for Fig.1(c)},$ (36) ${\cal
M}_{1,d}^{s,\mu\nu}=\tau_{1}^{\mu\nu}(p_{1},k-p_{2})\frac{i}{(k-p_{2})^{2}-m_{a}^{2}}\tau_{1,\alpha\beta}(k-p_{2},-p_{2})\epsilon^{s,\alpha\beta}\
\ \text{for Fig.1(d)}.$ (37)

Thus we get

$M_{1}^{\mu\nu}=\sum_{s=\pm}(M_{1,a}^{s,\mu\nu}+M_{1,b}^{s,\mu\nu}+M_{1,c}^{s,\mu\nu}+M_{1,d}^{s,\mu\nu}).$
(38)

The polarization tensor of the physical graviton can be described in terms of
a direct product of unit spin polarization vectors Bjerrum-Bohr:2014lea

$\displaystyle\epsilon^{+}_{\mu\nu}$ $\displaystyle=$
$\displaystyle\epsilon^{+}_{\mu}\epsilon^{+}_{\nu},\ \ \ \text{for
helicity}=+2,$ $\displaystyle\epsilon^{-}_{\mu\nu}$ $\displaystyle=$
$\displaystyle\epsilon^{-}_{\mu}\epsilon^{-}_{\nu},\ \ \ \text{for
helicity}=-2.$ (39)

In the harmonic gauge, they satisfy

$\eta^{\mu\nu}\epsilon^{+}_{\mu}\epsilon^{+}_{\nu}=\eta^{\mu\nu}\epsilon^{-}_{\mu}\epsilon^{-}_{\nu}=0,\
\ k^{\mu}\epsilon^{\pm}_{\mu}=0.$ (40)

We take the spin direction of the BH as the $z$-axis, then the unit
polarization vectors take the form

$\displaystyle\epsilon^{+}_{\mu}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}(0,-\cos{\theta}\cos{\phi}+i\sin\phi,-\cos\theta\sin\phi-i\cos\phi,\sin\theta),$
$\displaystyle\epsilon^{-}_{\mu}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}(0,\cos{\theta}\cos{\phi}+i\sin\phi,\cos\theta\sin\phi-i\cos\phi,-\sin\theta),$
(41)

where $\theta$ and $\phi$ are polar and azimuth angles of the momentum
$\boldsymbol{k}$ in the spherical coordinate.

##  III radiation power of Gravitational wave

The superradiant rate of different energy level is estimated as
Arvanitaki:2010sy

$\Gamma_{nlm}=2m_{a}\alpha^{4l+4}(GM_{b}+\sqrt{G^{2}M_{b}^{2}-J^{2}/M_{b}^{2}})\left(\frac{m}{2G^{2}M_{b}^{3}}\frac{J}{1+\sqrt{1-(J/GM_{b}^{2})^{2}}}-m_{a}\right)C_{nlm},$
(42)

where

$\displaystyle C_{nlm}$ $\displaystyle=$
$\displaystyle\frac{2^{4l+2}(3l+2)!}{(2l+2)^{2l+4}(l+1)!}\left(\frac{l!}{(2l)!(2l+1)!}\right)^{2}\prod_{j=1}^{l}\Bigg{(}j^{2}\left(1-\frac{J^{2}}{G^{2}M_{b}^{4}}\right)$
(43) $\displaystyle+$ $\displaystyle
4\left(GM_{b}+\sqrt{G^{2}M_{b}^{2}-J^{2}/M_{b}^{2}}\right)^{2}\left(\frac{m}{2G^{2}M_{b}^{3}}\frac{J}{1+\sqrt{1-(J/GM_{b}^{2})^{2}}}-m_{a}\right)^{2}\Bigg{)}.$

Thus the fastest superradiant level of axions in the cloud is $(n=2,l=1,m=1)$
where $n$, $l$, $m$ are the principal, orbital and magnetic quantum number.

In the following calculations, we can only consider the axions annihilating in
the $(n=2,l=1,m=1)$ state in the cloud, and the wave function of the
$(n=2,l=1,m=1)$ state can be solved by Schrödinger equation of the
gravitational atom,

$\Psi(\boldsymbol{x})=\frac{\sqrt{6}r}{12a_{0}^{5/2}}e^{-r/2a_{0}}\sqrt{\frac{3}{8\pi}}\sin\theta
e^{i\phi},$ (44)

where $a_{0}=1/(\alpha m_{a})$ is the Bohr radius of the gravitational atom
and the wave function in momentum space can be derived as

$\Psi(\boldsymbol{p})=-i2\sqrt{\pi}a_{0}^{3/2}\frac{a_{0}(p_{x}+ip_{y})}{(a_{0}^{2}\boldsymbol{p}^{2}+1/4)^{3}}.$
(45)

The axions annihilating to gravitons in the system of BH and axions cloud can
be considered as the process of $N(N-1)/2$ three-body bound states decay where
$N$ is the total number of axions in the state $(n=2,l=1,m=1)$. The radiation
power can be calculated by

$\displaystyle P$ $\displaystyle=$
$\displaystyle-\frac{dE}{dt}=\frac{N(N-1)}{2}(M-M_{b})\Gamma\simeq\frac{N^{2}}{2}(M-M_{b})\int\frac{d^{3}\boldsymbol{k}}{(2\pi)^{3}}\frac{\pi\delta(E_{1}+E_{2}-|\boldsymbol{k}|)}{|\boldsymbol{k}|}$
(46)
$\displaystyle\big{|}\iint\frac{d^{3}\boldsymbol{p_{1}}}{(2\pi)^{3}}\frac{d^{3}\boldsymbol{p_{2}}}{(2\pi)^{3}}\frac{\Psi_{1}(\boldsymbol{p_{1}})\Psi_{2}(\boldsymbol{p_{2}})}{\sqrt{2E_{1}}\sqrt{2E_{2}}}{\cal
M}_{1}^{\mu\nu}(\boldsymbol{p_{1}},\boldsymbol{p_{2}},\boldsymbol{k})\tilde{h}_{\mu\nu}(\boldsymbol{q})\big{|}^{2},$

where $M-M_{b}=2m_{a}(1-\alpha^{2}/8)$ for two axions annihilating from the
bound state $(n=2,l=1,m=1)$.

For the sake of $E_{1}\gg|\boldsymbol{p_{1}}|$ and
$E_{2}\gg|\boldsymbol{p_{2}}|$, we expand the propagators in the non-
relativistic limit to third order as

$\displaystyle\frac{1}{(p_{1}+p_{2})^{2}}$ $\displaystyle=$
$\displaystyle\frac{1}{2m_{a}^{2}+2E_{1}E_{2}-2\boldsymbol{p_{1}}\cdot\boldsymbol{p_{2}}}=\frac{1}{2}\Big{(}\frac{1}{m_{a}^{2}+E_{1}E_{2}}+\frac{\boldsymbol{p_{1}}\cdot\boldsymbol{p_{2}}}{(m_{a}^{2}+E_{1}E_{2})^{2}}$
(47) $\displaystyle+$
$\displaystyle\frac{(\boldsymbol{p_{1}}\cdot\boldsymbol{p_{2}})^{2}}{(m_{a}^{2}+E_{1}E_{2})^{3}}+\cdots\Big{)},$
$\displaystyle\frac{1}{(p_{i}-k)^{2}-m_{a}^{2}}=-\frac{1}{2}\frac{1}{E_{i}|\boldsymbol{k}|-\boldsymbol{p_{i}}\cdot\boldsymbol{k}}=-\frac{1}{2}\Big{(}\frac{1}{E_{i}|\boldsymbol{k}|}+\frac{\boldsymbol{p_{i}}\cdot\boldsymbol{k}}{E_{i}^{2}|\boldsymbol{k}|^{2}}+\frac{(\boldsymbol{p_{i}}\cdot\boldsymbol{k})^{2}}{E_{i}^{3}|\boldsymbol{k}|^{3}}+\cdots\Big{)},i=1,2.$
(48)
$\displaystyle\frac{1}{\boldsymbol{q}^{2}}=\frac{1}{(E_{1}+E_{2})^{2}}\Big{(}1-\frac{(\boldsymbol{p_{1}}+\boldsymbol{p_{2}})^{2}}{(E_{1}+E_{2})^{2}}+\frac{2\boldsymbol{k}\cdot(\boldsymbol{p_{1}}+\boldsymbol{p_{2}})}{(E_{1}+E_{2})^{2}}+\cdots\Big{)}.$
(49)

In order to get analytical form of the approximate results, we factorize the
integrand of $|\boldsymbol{p_{1}}|$ and $|\boldsymbol{p_{2}}|$ into the
product of the steep and smooth parts,

$\displaystyle\iint\frac{d^{3}\boldsymbol{p_{1}}}{(2\pi)^{3}}\frac{d^{3}\boldsymbol{p_{2}}}{(2\pi)^{3}}\frac{\Psi_{1}(\boldsymbol{p_{1}})\Psi_{2}(\boldsymbol{p_{2}})}{\sqrt{2E_{1}}\sqrt{2E_{2}}}{\cal
M}_{1}^{\mu\nu}(\boldsymbol{p_{1}},\boldsymbol{p_{2}},\boldsymbol{k})\tilde{h}_{\mu\nu}(\boldsymbol{q})=\iiiint
d|\boldsymbol{p_{1}}|d|\boldsymbol{p_{2}}|d\Omega_{1}d\Omega_{2}$ (50)
$\displaystyle\sin\theta_{1}(\cos\phi_{1}+i\sin\phi_{1})\sin\theta_{2}(\cos\phi_{2}+i\sin\phi_{2})I_{1}(|\boldsymbol{p_{1}}|,|\boldsymbol{p_{2}}|)I_{2}(\boldsymbol{p_{1}},\boldsymbol{p_{2}}),$

where $I_{1}$ and $I_{2}$ are the steep and smooth functions of
$|\boldsymbol{p_{1}}|$ and $|\boldsymbol{p_{2}}|$ respectively,

$\displaystyle I_{1}$ $\displaystyle=$
$\displaystyle\frac{a_{0}^{8}|\boldsymbol{p_{1}}|^{3}|\boldsymbol{p_{2}}|^{3}}{(a_{0}^{2}|\boldsymbol{p_{1}}|^{2}+1/4)^{3}(a_{0}^{2}|\boldsymbol{p_{2}}|^{2}+1/4)^{3}},$
(51) $\displaystyle I_{2}$ $\displaystyle=$
$\displaystyle-\frac{1}{a_{0}^{3}}\frac{1}{(2\pi)^{5}\sqrt{E_{1}E_{2}}}{\cal
M}_{1}^{\mu\nu}(\boldsymbol{p_{1}},\boldsymbol{p_{2}},\boldsymbol{k})\tilde{h}_{\mu\nu}(\boldsymbol{q}).$

The steep function $I_{1}$ is centered at
$|\boldsymbol{p_{1}}|=|\boldsymbol{p_{2}}|=1/(\sqrt{2}a_{0})$, thus we can
make the approximation:

$I_{1}\approx(\frac{8}{3\pi})^{2}\delta(|\boldsymbol{p_{1}}|-1/(\sqrt{2}a_{0}))\delta(|\boldsymbol{p_{2}}|-1/(\sqrt{2}a_{0})),$
(52)

where $(8/3\pi)^{2}$ is the normalized constant. After calculating the
integral, we finally get an analytical approximate form of the radiation
power,

Figure 3: The power of radiation with different “fine structure constant” in
the Kerr BH’s background field. The red line represents the numerical result
through the scattering amplitude method in our study. We compare this result
with the numerical result from Brito Brito:2017zvb  (black line), the
analytical results from Arvanitaki Arvanitaki:2010sy  (green line), and Brito
Brito:2014wla  (blue line).

$\displaystyle P$ $\displaystyle=$ $\displaystyle\frac{2.436\times
10^{-38}(M_{a}/\text{GeV})^{2}\alpha^{15}(2.977\times 10^{38}-3.721\times
10^{37}\alpha^{2})^{2}}{(M_{b}/\text{GeV})^{6}(2+\alpha^{2})^{11}(4+\alpha^{2})^{4}}\bigg{[}(M_{b}/\text{GeV})^{4}$
$\displaystyle\bigg{(}\frac{448.0}{\alpha}+2.583\times
10^{3}\alpha+6.828\times 10^{3}\alpha^{3}+1.094\times
10^{4}\alpha^{5}+1.182\times 10^{4}\alpha^{7}$ $\displaystyle+$ $\displaystyle
9.086\times 10^{3}\alpha^{9}+5.092\times 10^{3}\alpha^{11}+2.097\times
10^{3}\alpha^{13}+629.4\alpha^{15}+134.4\alpha^{17}$ $\displaystyle+$
$\displaystyle
19.36\alpha^{19}+1.691\alpha^{21}+0.068\alpha^{23}\bigg{)}+J(M_{b}/\text{GeV})^{2}\sqrt{2+\alpha^{2}}(-1.257\times
10^{41}$ $\displaystyle-$ $\displaystyle 6.915\times
10^{41}\alpha^{2}-1.746\times 10^{42}\alpha^{4}-2.675\times
10^{42}\alpha^{6}-2.771\times 10^{42}\alpha^{8}$ $\displaystyle-$
$\displaystyle 2.047\times 10^{42}\alpha^{10}-1.105\times
10^{42}\alpha^{12}-4.392\times 10^{41}\alpha^{14}-1.277\times
10^{41}\alpha^{16}$ $\displaystyle-$ $\displaystyle 2.648\times
10^{40}\alpha^{18}-3.719\times 10^{39}\alpha^{20}-3.176\times
10^{38}\alpha^{22}-1.248\times 10^{37}\alpha^{24})$ $\displaystyle+$
$\displaystyle J^{2}\alpha(1.764\times 10^{79}+1.012\times
10^{80}\alpha^{2}+2.687\times 10^{80}\alpha^{4}+4.378\times 10^{80}\alpha^{6}$
$\displaystyle+$ $\displaystyle 4.886\times 10^{80}\alpha^{8}+3.944\times
10^{80}\alpha^{10}+2.370\times 10^{80}\alpha^{12}+1.074\times
10^{80}\alpha^{14}$ $\displaystyle+$ $\displaystyle 3.668\times
10^{79}\alpha^{16}+9.330\times 10^{78}\alpha^{18}+1.717\times
10^{78}\alpha^{20}+2.166\times 10^{77}\alpha^{22}$ $\displaystyle+$
$\displaystyle 1.677\times 10^{76}\alpha^{24}+6.023\times
10^{74}\alpha^{26})\bigg{]}\text{GeV}^{2},$

where $M_{a}$ is the mass of the axion cloud and we take the number of axions
in the $(n=2,l=1,m=1)$ state as $N\approx M_{a}/m_{a}$. Since $\alpha<1$, we
can have the further approximate form

$\displaystyle P$ $\displaystyle=$
$\displaystyle\frac{(M_{a}/\text{GeV})^{2}\alpha^{14}}{(M_{b}/\text{GeV})^{6}(2+\alpha^{2})^{11}(4+\alpha^{2})^{4}}\bigg{[}(M_{b}/\text{GeV})^{4}(9.671\times
10^{41}+5.577\times 10^{42}\alpha^{2}$ $\displaystyle+$ $\displaystyle
1.474\times 10^{43}\alpha^{4}+2.361\times
10^{43}\alpha^{6})+J(M_{b}/\text{GeV})^{2}\alpha(-3.839\times 10^{80}$
$\displaystyle-$ $\displaystyle-2.111\times 10^{81}\alpha^{2}-5.329\times
10^{81}\alpha^{4}-8.165\times 10^{81}\alpha^{8})+J^{2}\alpha^{2}(3.809\times
10^{118}$ $\displaystyle+$ $\displaystyle 2.184\times
10^{119}\alpha^{2}+5.799\times 10^{119}\alpha^{4}+9.450\times
10^{119}\alpha^{6})\bigg{]}~{}\text{GeV}^{2}.$

Figure 4: Left: The ratio of radiation power with different spin of BH $J$ to
the spinless case $J=0$. The blue, orange, green, red lines correspond to
$J=$0.05, 0.10, 0.15, 0.20$GM_{b}^{2}$ respectively. Right: The ratio of
analytical approximate results in Eqs. III (black) and III (orange) to
numerical result at $J=0.1GM_{b}^{2}$.

We also calculate the numerical result of radiation power through the Monte
Carlo method. The radiation power of GW with different “fine structure
constant” is shown in Fig. 3 where we take spin of the BH
$J=4\frac{GM_{b}^{2}\alpha(1-\alpha^{2}/8)}{1+4\alpha^{2}(1-\alpha^{2}/8)^{2}}$
as in Ref. Brito:2017zvb . We compare the result through the scattering
amplitude method with the traditional method from Brito Brito:2017zvb
Brito:2014wla  and Arvanitaki Arvanitaki:2010sy , our result is a good
approximation to numerical result from Brito Brito:2017zvb  at small $\alpha$
as expected. The analytic result derived from Brito Brito:2014wla  is obtained
in the background of spinless BH which correspond to take
$\tilde{h}_{\mu\nu,J}=0$ (Eq. 31) approximation in the scattering amplitude
method. While the analytical result from Arvanitaki Arvanitaki:2010sy  is
around two to five orders of magnitude smaller than the numerical result since
it is calculated in the flat space. In fact, the formula used in the
traditional method in the flat space background is equivalent to the formula
of radiation power generated by the classical source. We postpone the detailed
discussion in Sec. IV.

The ratio of the magnitude of different BH spin angular momentum $J$ to the
spinless case $J=0$ is shown in the left panel of Fig. 4. We can found the
radiation power decreases with the increase of the BH spin angular momentum
when spin is small. To verify the validity of Eqs. III and III, we also plot
the ratio of the analytical approximate results to the numerical result where
we take $J=0.1GM_{b}^{2}$ in the right panel of Fig. 4. They are good
approximations within the range of $\alpha$ and $J$ considered in this work.

##  IV Discussion

We note that the traditional method used in Ref. Arvanitaki:2010sy  is based
on the formula for the GW radiation power,

$\frac{dP}{d\Omega}=\frac{G\omega^{2}}{\pi}T_{ij}^{TT*}(\omega,\boldsymbol{k})T_{ij}^{TT}(\omega,\boldsymbol{k}),$
(55)

where

$T_{ij}(\omega,k)\equiv\int
d^{3}\boldsymbol{x}T_{ij}(\omega,\boldsymbol{x})e^{-i\boldsymbol{k}\cdotp\boldsymbol{x}},$
(56)

with $|k|=\omega$ and the superscript $TT$ means that we take the transverse-
traceless part,

$T_{ij}^{TT}\equiv(P_{ii^{\prime}}P_{jj^{\prime}}-\frac{1}{2}P_{ij}P_{i^{\prime}j^{\prime}})T_{i^{\prime}j^{\prime}},$
(57)

where

$P_{ij}=\delta_{ij}-\frac{k_{i}k_{j}}{\boldsymbol{k}^{2}}.$ (58)

Figure 5: The Feynman diagram for scattering amplitude of
$|\boldsymbol{p_{1}},\boldsymbol{p_{2}}\rangle\rightarrow|\boldsymbol{k}\rangle$
process in classical source limit.

We now derive this formula from field theory in classical limit. Since this
formula is applicable in both transition and annihilation processes, we just
derive it from transition process for simplicity. We can consider a classical
source emits a graviton as a transition process
$|p_{1}\rangle\rightarrow|p_{2},k\rangle$ where $k$ is the momentum of emitted
graviton, and $p_{1}$, $p_{2}$ are the momenta of the source in initial and
final states. The decay width can be calculated similarly as Eq. 10,

$\Gamma\simeq\frac{1}{2M_{j}}\int\frac{d^{3}\boldsymbol{k}}{(2\pi)^{3}}\frac{1}{2|\boldsymbol{k}|}\frac{1}{2M_{j}}\big{|}{\cal
M}(p_{1}\rightarrow p_{2},k)|^{2}2\pi\delta(E_{1}-E_{2}-|\boldsymbol{k}|),$
(59)

where $M_{j}$ is the mass of the classical source. And the scattering
amplitude ${\cal M}(p_{1}\rightarrow p_{2},k)$ is just that we replace the
propagator of graviton and $M_{b}$ with polarization tensor of physical
graviton and $M_{j}$ respectively in $-i{\cal M}_{2}^{\mu\nu}$. As shown in
Fig. 5, in classical source limit,

${\cal M}(p_{1}\rightarrow p_{2},k)=-i\kappa
M_{j}T_{\mu\nu}(\omega,\boldsymbol{k})\epsilon^{\mu\nu},$ (60)

where $\omega=|\boldsymbol{k}|$, then we get

$\Gamma=\sum_{s=\pm}\int\frac{d\omega
d\Omega\omega\kappa^{2}}{32\pi^{2}}|T_{\mu\nu}(\omega,\boldsymbol{k})\epsilon^{s,\mu\nu}|^{2}\delta(E_{1}-E_{2}-\omega),$
(61)

the polarization sum
$\sum_{s_{1},s_{2}=\pm}\epsilon^{s_{1},\mu\nu*}\epsilon^{s_{2},\alpha\beta}$can
be calculated by using Eq. II.3. Finally we get

$\Gamma=\int
d\Omega\frac{G\omega_{12}}{\pi}T_{ij}^{TT*}(\omega_{12},\boldsymbol{k})T_{ij}^{TT}(\omega_{12},\boldsymbol{k}),$
(62)

where $\omega_{12}=E_{1}-E_{2}$, thus the radiation power is

$\frac{dP}{d\Omega}=\frac{G\omega_{12}^{2}}{\pi}T_{ij}^{TT*}(\omega_{12},\boldsymbol{k})T_{ij}^{TT}(\omega_{12},\boldsymbol{k}).$
(63)

This formula is only applicable in classical source limit in flat space. But
for axions, process $|p_{1}\rangle\rightarrow|p_{2},k\rangle$ is kinematically
forbidden due to the on-shell requirement in flat space. Thus the flat space
approximation in Ref. Arvanitaki:2010sy  is equivalent to the classical source
limit in which the matter has no on-shell requirement. For the real physical
process, we must consider the effect of the background field. We need to take
the Kerr BH into account for this process as the Feynman diagrams in Fig. 1
and then treat it as a classical source. This is the method we have used in
this work.

##  V conclusion

Using the method of scattering amplitude, we have calculated the GW radiation
power of the axions annihilation process in the axion cloud around Kerr BH.
Compared with the traditional method for calculating the GW radiation power,
the scattering amplitude method has the following advantages:

1\. The calculation is simple and straightforward. The Hamiltonian of the
traditional method is $H=H_{0}+H_{I}$ with
$H_{0}=H_{0,a}+H_{bf,a},H_{I}=H_{int,a}$, where $H_{0,a}$ and $H_{int,a}$
represent the free Hamiltonian and gravitational self-interaction of axion
field. $H_{bf,a}$ represents the interaction between the background and axion
field. The traditional method is equivalent to derive the eigenstates of
$H_{0}$ first and then calculate the transition rate of different eigenstates
under perturbation $H_{I}$. While for the method of scattering amplitude, we
take the form of Hamiltonian as $H_{0}=H_{0,a}+H_{0,b}$,
$H_{I}=H_{a,b}+H_{int,a}$, where $H_{0,b}$ is the free Hamiltonian of BH and
$H_{a,b}$ represents the interaction between Kerr BH and axion particle.
Although the perturbation method of scattering amplitude is valid for small
$\alpha$, the numerical calculation is much easier than the traditional
method.

2\. It is natural and easy to include the Kerr metric effects of the Kerr BH.
We can divide the calculation of scattering amplitude into two parts as in Eq.
13, where $M_{1\mu\nu}$ stands for axions annihilation process and
$M_{2\mu\nu}$ corresponds to metric effects induced by Kerr BH. The metric
effects of Kerr BH act as a background field
$\tilde{h}_{\mu\nu}=\tilde{h}_{\mu\nu,0}+\tilde{h}_{\mu\nu,J}$ which is
included in $M_{2\mu\nu}$ as in Eq. 24. Since the scattering method is valid
for small $\alpha$, we only expand the metric up to first non-vanish terms in
the weak-field approximation. And note that only the first-order term of BH
spin $J$ in Eq. 28 is retained. In principle, when the spin of Kerr BH is
large, higher power terms of spin $J$ needs to be included in Eq. 28 for a
more accurate result and we can get these terms by expanding Kerr metric
directly. For simplicity, we only retain linear terms of $J$ which is a good
approximation for the small BH spin.

3\. The microscopic physical process is intuitive. The annihilation process of
two axions can be considered as a decay process of a three-body bound state.
For this decay process, there are four Feynman diagrams corresponding to
different scattering channels between particles. We can naturally extend our
calculation to four, six, or more axions annihilation processes.

4\. It is clear and simple to demonstrate the analytic approximation formulae.
Since the wave function of the axion state $(n=2,l=1,m=1)$ is centered at
$p\sim\alpha m_{a}/2$, we can get an analytic approximate formula with the
spin effect from the phase integration of scattering amplitude. While the spin
effect is hard to be included in the analytic approximation in the traditional
method.

We have cross-checked and confirmed the numerical results from the traditional
method Brito:2017zvb  which is important for the GW experiments and axion
search. More precise calculations, implications for GW or axion detection, and
more broad applications are working in progress.

###### Acknowledgements.

This work was supported by the National Natural Science Foundation of China
(NNSFC) under Grant No. 12205387.

## References

  * (1) R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440-1443 (1977) doi:10.1103/PhysRevLett.38.1440; S. Weinberg, Phys. Rev. Lett. 40, 223-226 (1978) doi:10.1103/PhysRevLett.40.223; F. Wilczek, Phys. Rev. Lett. 40, 279-282 (1978) doi:10.1103/PhysRevLett.40.279; J. E. Kim, Phys. Rev. Lett. 43, 103 (1979) doi:10.1103/PhysRevLett.43.103; M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, Nucl. Phys. B 166, 493-506 (1980) doi:10.1016/0550-3213(80)90209-6; A. R. Zhitnitsky, Sov. J. Nucl. Phys. 31, 260 (1980); M. Dine, W. Fischler and M. Srednicki, Phys. Lett. B 104, 199-202 (1981) doi:10.1016/0370-2693(81)90590-6; P. Sikivie, Rev. Mod. Phys. 93, no.1, 015004 (2021) doi:10.1103/RevModPhys.93.015004 [arXiv:2003.02206 [hep-ph]].
  * (2) P. Svrcek and E. Witten, JHEP 06, 051 (2006) doi:10.1088/1126-6708/2006/06/051 [arXiv:hep-th/0605206 [hep-th]].
  * (3) J. Preskill, M. B. Wise and F. Wilczek, Phys. Lett. B 120, 127-132 (1983) doi:10.1016/0370-2693(83)90637-8; M. Dine and W. Fischler, Phys. Lett. B 120, 137-141 (1983) doi:10.1016/0370-2693(83)90639-1
  * (4) P. Sikivie, Lect. Notes Phys. 741, 19-50 (2008) doi:10.1007/978-3-540-73518-2_2 [arXiv:astro-ph/0610440 [astro-ph]].
  * (5) R. Penrose, Riv. Nuovo Cim. 1, 252-276 (1969) doi:10.1023/A:1016578408204
  * (6) Zel’Dovich, Ya. B. Journal of Experimental and Theoretical Physics Letters, Vol. 14, (1971)
  * (7) A. A. Starobinsky, Sov. Phys. JETP 37, no.1, 28-32 (1973)
  * (8) S. L. Detweiler, Phys. Rev. D 22, 2323-2326 (1980) doi:10.1103/PhysRevD.22.2323
  * (9) R. Brito, V. Cardoso and P. Pani, Physics,” Lect. Notes Phys. 906, pp.1-237 (2015) 2020, ISBN 978-3-319-18999-4, 978-3-319-19000-6, 978-3-030-46621-3, 978-3-030-46622-0 doi:10.1007/978-3-319-19000-6 [arXiv:1501.06570 [gr-qc]].
  * (10) A. Arvanitaki, M. Baryakhtar and X. Huang, Phys. Rev. D 91, no.8, 084011 (2015) doi:10.1103/PhysRevD.91.084011 [arXiv:1411.2263 [hep-ph]].
  * (11) A. Arvanitaki, M. Baryakhtar, S. Dimopoulos, S. Dubovsky and R. Lasenby, Phys. Rev. D 95, no.4, 043001 (2017) doi:10.1103/PhysRevD.95.043001 [arXiv:1604.03958 [hep-ph]].
  * (12) M. Baryakhtar, R. Lasenby and M. Teo, Phys. Rev. D 96, no.3, 035019 (2017) doi:10.1103/PhysRevD.96.035019 [arXiv:1704.05081 [hep-ph]].
  * (13) P. S. B. Dev, M. Lindner and S. Ohmer, Phys. Lett. B 773, 219-224 (2017) doi:10.1016/j.physletb.2017.08.043 [arXiv:1609.03939 [hep-ph]].
  * (14) N. Xie and F. P. Huang, [arXiv:2207.11145 [hep-ph]].
  * (15) R. Brito, S. Ghosh, E. Barausse, E. Berti, V. Cardoso, I. Dvorkin, A. Klein and P. Pani, Phys. Rev. D 96, no.6, 064050 (2017) doi:10.1103/PhysRevD.96.064050 [arXiv:1706.06311 [gr-qc]].
  * (16) A. Arvanitaki and S. Dubovsky, Phys. Rev. D 83, 044026 (2011) doi:10.1103/PhysRevD.83.044026 [arXiv:1004.3558 [hep-th]].
  * (17) R. Brito, V. Cardoso and P. Pani, Class. Quant. Grav. 32, no.13, 134001 (2015) doi:10.1088/0264-9381/32/13/134001 [arXiv:1411.0686 [gr-qc]].
  * (18) H. Yoshino and H. Kodama, PTEP 2014, 043E02 (2014) doi:10.1093/ptep/ptu029 [arXiv:1312.2326 [gr-qc]].
  * (19) B. S. DeWitt, Phys. Rev. 162, 1239-1256 (1967) doi:10.1103/PhysRev.162.1239
  * (20) N. E. Bjerrum-Bohr, [arXiv:hep-th/0410097 [hep-th]].
  * (21) N. E. J. Bjerrum-Bohr, B. R. Holstein, L. Planté and P. Vanhove, Phys. Rev. D 91, no.6, 064008 (2015) doi:10.1103/PhysRevD.91.064008 [arXiv:1410.4148 [gr-qc]].