# Performance of quantum registers in diamond in the presence of spin
impurities

Dominik Maile  Joachim Ankerhold  Institute for Complex Quanten Systems and
IQST, Universität Ulm, D-89069 Ulm, Germany

(August 27, 2024)

###### Abstract

The Nitrogen Vacancy Center in diamond coupled to addressable surrounding
nuclear spins forms a versatile building block for future quantum
technologies. While previous activities focused on sensing with only a single
or very few spins in operation, recently multi-qubit registers have been
successfully implemented for quantum information processing. Further progress
requires a detailed understanding of the performance of quantum protocols for
consecutive gate operations and thus, beyond established treatments for
relaxation and dephasing. Here, we provide such an analysis for a small spin
registers with up to four spins built out of NV and environmental constituents
in presence of ensembles of interacting impurity spins. Adapting a cluster
correlation expansion, we predict coherence properties as well as fidelities
for GHZ- and Bell-gate operations. The influence of the volume density and the
geometry of the spin-bath consisting of substitutional nitrogen atoms are also
taken into account.

##  I Introduction

Decoherence and relaxation are prime obstacles for the usability of quantum
mechanical attributes of a physical system. Over the last decades, theoretical
and experimental physicists proposed a zoo of different platforms to overcome
the fragility of information encoded in quantum superposition states. Although
several small quantum processors have been built from different platforms -
superconducting Josephson arrays and arrays of trapped ions being the most
prominent ones - each of them is subject to different limitations with respect
to scalability and coherence times. In this respect, the Nitrogen Vacancy
center (NV) in diamond is in principle a promising building block for quantum
technologies offering long coherence times and efficient optical readout even
at ambient conditions Doherty  _et al._ (2013); Pezzagna and Meijer (2021);
Wolfowicz  _et al._ (2021). While coherent processes with a single NV-center
were already realized almost two decades ago Jelezko  _et al._ (2004a) and are
now routinely exploited for sensing purposes, scaling up to larger NV-
registers remains a challenge Dolde  _et al._ (2013); Haruyama  _et al._
(2019).

Figure 1: (a) NV-based setting of a quantum register studied in this work. The
bare register (encircled) consists of one $\text{NV}^{-}$ (with one electron
spin and one nitrogen nuclear spin) and two ${}^{13}C$ atoms. The
$\text{NV}^{-}$ electron spin is located in the center of the coordinate
system and the coordinates of the ${}^{13}C$ spins are $(0.87,0,0.19)$nm and
$(0.56,0.77,0.31)$nm. Spins outside of the circle represent the spin-bath
consisting of nitrogen impurities (P1). (b) General form of the pulse
structure used to initiate operations on the register. We assume each pulse to
be an instantaneous rotation on the Bloch sphere of the $\text{NV}^{-}$ given
by the angles $\theta_{i},\phi_{i}$. Between pulses, register+spin-bath evolve
unitarily for times $t_{i}$. Parameters $\theta_{i},\phi_{i},t_{i}$ are found
by optimizing the fidelity with respect to a specific target state.

To be able to coherently couple individual NV-centers via their dipole
interaction efficiently, they have to be located at a distance of a few tens
of nanometers Jakobi  _et al._ (2016). Since this distance dependence cannot
be controlled with sufficient precision yet, the assembling of larger
registers relies on stochastic processes during the diamond growth or
implementation of nitrogen Groot-Berning  _et al._ (2021); Haruyama  _et al._
(2019); Jakobi  _et al._ (2016); Yamamoto  _et al._ (2013). Given this
situation, the focus has been laid on other magnetic impurities found in
diamond, predominantly nitrogen atoms (P1-impurites) and carbon ${}^{13}C$
isotopes. While P1-impurities yield an additional electron spin in the
lattice, the ${}^{13}C$ isotopes have a nuclear spin. Nitrogen is a necessary
ingredient to obtain NV-centers in the first place, either during the diamond
growth process or by implantation. Further, there is evidence, that nitrogen
is also an electron donor turning the uncharged $\text{NV}^{(0)}$ to the
technologically useful negatively charged $\text{NV}^{-}$ making close
P1-impurities unavoidable Doherty  _et al._ (2013). Overall, the creation
yield of these nitrogen impurities to pair with vacancies and to become NV-
centers, plays a crucial role. Too many remaining pure nitrogen impurities
form a spin-bath and lead to fast decoherence for the register as the
surrounding electron spins couple strongly to the NV-center Bauch  _et al._
(2020). Also additional ${}^{13}C$ atoms have a two-sided influence on the
performance of a diamond based spin register: Depending on their abundance and
individual position, they may act as a spin-bath or as addressable qubits. In
fact, the ${}^{13}C$’s which are strongly coupled to the $\text{NV}^{-}$ can
be initialized to the desired state and are addressable via either radio-
frequency drives or via the dipole interaction with the $\text{NV}^{-}$
Jelezko  _et al._ (2004b); Childress  _et al._ (2006); Neumann  _et al._
(2008); Jacques  _et al._ (2009). The long coherence properties of the nuclear
spins also lead to the proposal to use them as quantum memory or quantum
repeaters for quantum information networks Bernien  _et al._ (2013); Hensen
_et al._ (2015); Bradley  _et al._ (2021); Pompili  _et al._ (2021).

Based on the idea of using several ${}^{13}C$’s as being part of the register,
quantum gates and small quantum algorithms have successfully been implemented
experimentally Hegde  _et al._ (2020); Zhang  _et al._ (2020); Bradley  _et
al._ (2021). These algorithms can be initiated by local control of the full
register: MW-driving the $\text{NV}^{-}$ and additional radiofrequency drives
for the ${}^{13}C$. Here, only the multi-qubit processes rely on the dipole
interaction between the ${}^{13}C$ and the $\text{NV}^{-}$. Most
spectacularly, the coherent control of a quantum register consisting of a
single NV-center and ten nuclear spins (including the nuclear spin of the NV)
was recently realized in such a way Bradley  _et al._ (2019). On the other
hand, a minimal control scheme uses only the MW-drive on the $\text{NV}^{-}$
to perform processes on the full register. Here, also one-qubit gates on the
${}^{13}C$’s are mediated by the dipole interaction with the $\text{NV}^{-}$
Hegde  _et al._ (2020); Zhang  _et al._ (2020).

The main infidelity of NV-registers for processes on the timescale of several
microseconds and especially at lower temperature can be attributed to the
surrounding spin-bath Bar-Gill  _et al._ (2012, 2013), where details strongly
depend on the abundance and the specific distribution of impurity spins.
Further progress towards quantum computing architectures thus necessitates a
detailed understanding of decoherence induced by these environmental degrees
of freedom beyond mean field approximations and their impact on gate
operations beyond single-qubit processes. These are the main topics of the
present work.

From the theory side, a well-explored setting consists of a central spin
embedded in a spin-bath. de Sousa and Das Sarma (2003); Yao  _et al._ (2006);
Witzel and Das Sarma (2006); Breuer (2004); Dobrovitski  _et al._ (2008); Zhao
_et al._ (2012); Yang and Liu (2008); Hall  _et al._ (2014); Myers  _et al._
(2017); Onizhuk  _et al._ (2021); Kwiatkowski and Cywiński (2018); Bragar and
Cywiński (2015); Ivády (2020); Norambuena  _et al._ (2018); Cambria  _et al._
(2021). Starting initially from a dedicated quantum state, the main body of
work has considered pure dephasing processes only, particularly relevant for
sensing applications. Here, a microscopic approach that takes a specific
configuration of the spin-bath into account, uses a ”cluster correlation
expansion” (CCE) to subsequently include higher correlations within the bath
Witzel and Das Sarma (2006); Zhao  _et al._ (2012). This technique was
recently generalized to the ”generalized cluster correlation expansion” (gCCE)
which can be used to calculate the full density matrix of a system also
capturing relaxation effects Onizhuk  _et al._ (2021).

In this work, we go substantially beyond these previous achievements in that
we adapt the gCCE to simulate the dynamical performance of a complete spin
register in diamond in the presence of interacting magnetic impurities. This
way, we also demonstrate the relevance of environmental spin-spin correlations
for multi-qubit register and sequences of gate operations, and provide
fidelity maps for four qubit quantum processes.

More specifically, we analyze a quantum register consisting of one
$\text{NV}^{-}$ (containing one electron spin and the nuclear spin of the
nitrogen) and two ${}^{13}C$’s embedded in a spin-bath consisting of
P1-impurity spins, see Fig. 1a. Further, we assume a purified diamond setup
111In natural diamond the abundance $\eta_{{}^{13}C}=1.1\%$ yielding to a
reasonable probability to find strongly coupled ${}^{13}C$’s ($\sim
100\text{MHz}$). Purified diamond (e.g. $\eta_{{}^{13}C}\approx 0.01\%$) leads
to a sparse ${}^{13}C$ bath with better coherence properties for the NV spin.,
where two ${}^{13}C$ are interacting individually with the $\text{NV}^{-}$ via
dipole couplings with several hundreds of kilohertz. Since we calculate the
full evolution of the register starting from a chosen initialized quantum
state and driven by pulse sequences as described in Fig. 1b, i.e. far beyond
the investigation of pure decoherence or relaxation phenomena, we obtain a
comprehensive picture of realistic quantum operations.

We structure the article as follows: In Sec. II, we start by introducing the
full Hamiltonian of the register and the spin-bath. Then, we describe our
theoretical modelling of the pulse structure in detail. Next, we a give short
introduction into the gCCE. In Sec. III, we first truncate the full register
to a two qubit register and benchmark the gCCE with exact diagonalization for
a quantum process in a spin-bath consisting of nine impurities. Then, we study
the performance of different quantum processes on a fully connected four qubit
register in presence of different spin-baths. We conclude in Sec. IV and give
an outlook how to overcome current limitations.

##  II Theory and Methods

###  II.1 The Hamiltonian of the system

We consider a spin register composed out of the $\text{NV}^{-}$ electron spin,
the nuclear spin of the nitrogen atom and two nuclear ${}^{13}\text{C}$ spins.
In principle, P1-impurities at a distance $r<30\,$nm from the NV yield
interaction strengths comparable to the ${}^{13}C$ nuclear spins which,
however, would destroy the register in the fist place. Accordingly, we assume
a situation, where residual bath-spins reside outside a central core, see Fig.
1a. The full Hamiltonian thus reads

$\displaystyle H=H_{R}+H_{B},$ (1)

where $H_{B}$ is the Hamiltonian for the environment and $H_{R}$ the one for
the register containing three parts

$\displaystyle H_{R}=H_{NV}+H_{N}+\sum_{k}H^{(k)}_{{}^{13}\text{C}}.$ (2)

For the purpose of this work, the bare $\text{NV}^{-}$ center can be readily
described by the model Hamiltonian

$\displaystyle H_{NV}=D_{gs}\,S_{z}^{2}-g_{e}\,B_{z}S_{z},$ (3)

where $D_{gs}=2.88\,$GHz is the zero field splitting of the triplet state,
$g_{e}=-2.806\,$MHz/G is the gyromagnetic ratio of the electron spin, and
$B_{z}$ is an external magnetic field. For the results shown below, we
consider $B_{z}=148\,$G as in Hegde  _et al._ (2020). We are interested in the
effect of the environment on the quantum processes initiated by a microwave
drive of the $\text{NV}^{-}$ spin. As we will show later, these processes take
several microseconds. In this regime the longitudinal relaxation to the spin
environment is irrelevant. This justifies a truncation of the $\text{NV}^{-}$
triplet state to a two - level system built out of the $m_{s}=0$ and the
$m_{s}=-1$ state. In the following, we denote the matrices of the
$\text{NV}^{-}$ spin

$\displaystyle\\!\\!S_{z}\\!=\\!\left(\\!\\!\begin{array}[]{cc}0&0\\\
0&-1\end{array}\\!\\!\right),\,S_{x}\\!=\\!\frac{1}{\sqrt{2}}\left(\\!\\!\begin{array}[]{cc}0&1\\\
1&0\end{array}\\!\\!\right),\,S_{y}\\!=\\!\frac{1}{\sqrt{2}i}\left(\\!\\!\begin{array}[]{cc}0&1\\\
-1&0\end{array}\\!\\!\right)\\!.$ (10)

Further, the Hamiltonian encoding the triplet nuclear spin state of the
nitrogen ($I_{N}=1$) reads

$\displaystyle
H_{N}=P_{gs}\;I_{N,z}^{2}-g_{n}B_{z}I_{N,z}+\boldsymbol{S}\bar{N}\boldsymbol{I}_{N}$
(11)

where $P_{gs}=-5.08\text{ MHz}$ is the zero field splitting,
$g_{n}=0.308\,$kHz/G the gyromagnetic ratio and $\bar{N}$ the hyperfine
coupling tensor between the nitrogen nuclear spin and the $\text{NV}^{-}$
electron spin. The bold quantities denote the full spin vectors. The
${}^{13}C$’s have nuclear spin $I_{C}=1/2$ and are described by the
Hamiltonian

$\displaystyle
H^{(k)}_{{}^{13}C}=-g_{C}B_{z}I^{(k)}_{C,z}+\boldsymbol{S}\bar{M}^{(k)}\boldsymbol{I}^{(k)}_{C}+\boldsymbol{I}^{(k)}_{C}\bar{L}^{(k)}\boldsymbol{I}_{N},$
(12)

where $g_{C}=1.071\,$kHz/G is the gyromagnetic ratio, $\bar{M}^{(k)}$ the
hyperfine coupling between the $\text{NV}^{-}$ and the $k$-th ${}^{13}C$ and
$\bar{L}^{(k)}$ is the hyperfine coupling between the $k$-th nuclear spin of
the ${}^{13}C$ and the nitrogen nuclear spin. All dipole interaction tensors
are further discussed in App. A. We note that the dipole interactions between
the nuclear spins are irrelevant on the timescales we discuss in this article.
However, we take them into account in our simulation since they might become
important for long sequences of quantum operations we plan to simulate in the
future.

This register is assumed to be embedded in a finite size spin-bath consisting
of P1-impurities. The latter are substitutional nitrogen impurities carrying
an extra electron spin. Due to the large gyromagnetic ratio of the electron
spins such P1-impurities couple strongly to the spin register and other
magnetic impurities. We emphasize that the quantum register is embedded in a
common spin-bath with no simplification on the spin-spin interactions. The
general Hamiltonian for the environment thus reads

$\displaystyle H_{B}=$
$\displaystyle\sum_{j}\left(-g_{e}B_{z}\mathcal{E}^{(j)}_{z}+\boldsymbol{S}\bar{K}^{(j)}\boldsymbol{\mathcal{E}}^{(j)}+\boldsymbol{I}_{N}\bar{A}^{(j)}\boldsymbol{\mathcal{E}}^{(j)}\right.$
(13)
$\displaystyle\left.+\sum_{k}\boldsymbol{I}^{(k)}_{C}\bar{F}^{(k,j)}\boldsymbol{\mathcal{E}}^{(j)}\right)+\sum_{ij}\boldsymbol{\mathcal{E}}^{(i)}\bar{G}^{(ij)}\boldsymbol{\mathcal{E}}^{(j)},$

where $\bar{K}^{(j)}$ and $\bar{A}^{(j)}$ are the dipole interaction tensor
between the electron spin originating from the P1-impurity (denoted by
$\boldsymbol{\mathcal{E}}$) and the $\text{NV}^{-}$ electron spin and its
nuclear spin respectively. $\bar{F}^{(k,j)}$ is the dipole coupling tensor
between the ${}^{13}C$ nuclear spins of the register and the spin-bath and
$\bar{G}^{(i,j)}$ the interaction within the spin environment.

Figure 2: (a) Influence of the pulse sequence on the $I_{x}$ component of the
${}^{13}C$ nuclear spin. The vertical dotted line shows the end of the gate
operation $t_{e}$. (b) $z$-component of the ${}^{13}C$ spin. (c) Logarithmic
negativity as a measure of the entanglement in the system. The deviation of
the dissipative evolution from the unitary is a measure for the fidelity of
the process. The deviation grows for longer times. Each mean field evolution
was averaged over $M=450$ bath-states. The grey line in each plot shows the
$S_{z}$ component of the NV electron spin under the influence of the pulse
sequence. The ticks for $S_{z}$ are displayed on the right axis of (c).

###  II.2 Quantum processes on larger spin registers

The operations on the quantum register are initiated solely by driving the
electron spin of the NV center. We follow protocols suggested by several
experiments in this direction Hegde  _et al._ (2020); Zhang  _et al._ (2020).
Particularly, we consider pulse sequences of the form shown in Fig. 1b. The
full unitary evolution operator of the register is encoded in the parameters
$\\{t_{i},\theta_{i},\phi_{i}\\}$. Here, the angles $\theta_{i}$ and
$\phi_{i}$ describe a particular rotation on the Bloch sphere of the NV spin
and $t_{i}$ are times for which we let the system interact. To find the
parameters of the particular process, we encode the sequence of pulses into
subsequent unitary operations

$\displaystyle\mathcal{S}(\boldsymbol{t},\boldsymbol{\theta},\boldsymbol{\phi})\\!=\\!...U_{\theta_{4},\phi_{4}}U_{t_{3}}U_{\theta_{3},\phi_{3}}U_{t_{2}}U_{\theta_{2},\phi_{2}}U_{t_{1}}U_{\theta_{1},\phi_{1}}U_{t_{0}},$
(14)

where $U_{t_{i}}$ denotes the unitary evolution with the Hamiltonian $H$
(generally containing the spin-bath) and

$\displaystyle
U_{\theta_{j},\phi_{j}}=e^{-i\frac{\theta_{j}}{2}\left(\sigma_{x}\otimes\mathbb{1}_{v}\cos(\phi_{j})+\sigma_{y}\otimes\mathbb{1}_{v}\sin(\phi_{j})\right)}$
(15)

is encoding the pulse. $\sigma_{x}$ and $\sigma_{y}$ are Pauli matrices acting
on the NV-electron spin and $\mathbb{1}_{v}$ is the identity matrix for the
rest of the system. We optimize the infidelity
$\mathcal{I}_{p}=1-\mathcal{F}_{p}$, where $\mathcal{F}_{p}$ is the fidelity
of the state after the gate sequence without the dissipative couplings

$\displaystyle\rho_{p}=\mathcal{S}(\boldsymbol{t},\boldsymbol{\theta},\boldsymbol{\phi})\,\rho(0)\,\mathcal{S}^{\dagger}(\boldsymbol{t},\boldsymbol{\theta},\boldsymbol{\phi})$
(16)

with respect to the target state $\rho_{T}$ and defined as

$\displaystyle\mathcal{F}_{p}=\frac{\text{tr}\left(\rho_{T}^{\dagger}\rho_{p}\right)}{\text{tr}\left(\rho_{T}^{\dagger}\rho_{T}\right)}.$
(17)

Further, we assume the pulse unitaries $U_{\theta_{i},\phi_{i}}$ to be
instantaneous and perfect. For the optimization, we use the python build in
function scipy.basinhopping.

###  II.3 Generalized Cluster Correlation Expansion

To efficiently simulate the effect of the environment on the spin register, we
use the generalized Cluster Correlation Expansion (gCCE) taking the full
system Hamiltonian for the unitary evolution into account. Within this
methodology, higher order correlations within the environment are considered
by subsequently calculating higher order terms in the expansion. We start with
the well-established zeroth order (gCCE0) corresponding to the mean field
approximation and yielding the simplified bath Hamiltonian

$\displaystyle H^{(\text{E0})}_{B}=$
$\displaystyle\sum_{j}\left(-g_{e}B_{z}\tilde{\mathcal{E}}^{(j)}_{z}+S_{z}\bar{K}_{zz}^{(j)}\tilde{\mathcal{E}}_{z}^{(j)}+I_{N,z}\bar{A}_{zz}^{(j)}\tilde{\mathcal{E}}_{z}^{(j)}\right.$
$\displaystyle\left.+\sum_{k}I^{(k)}_{C,z}\bar{F}_{zz}^{(k,j)}\tilde{\mathcal{E}}_{z}^{(j)}\right)+\sum_{ij}\tilde{\mathcal{E}}_{z}^{(i)}\bar{G}_{zz}^{(ij)}\tilde{\mathcal{E}}_{z}^{(j)},$
(18)

where the bath spins $\tilde{\mathcal{E}}_{z}^{(j)}=\pm\frac{1}{2}$ are no
longer operators, but chosen by classical probability (static disorder) and E0
stands for gCCE0. From now on, a particular collective bath-state (e.g.
$\\{\tilde{\mathcal{E}}_{z}^{(1)},\tilde{\mathcal{E}}_{z}^{(2)},\tilde{\mathcal{E}}_{z}^{(3)},..,\tilde{\mathcal{E}}_{z}^{(N)}\\}=\\{\frac{1}{2},-\frac{1}{2},-\frac{1}{2},..,\frac{1}{2}\\}$)
is encoded by the superscript $n=1,\ldots 2^{N}$ counting all possible
reservoir configurations. Considering room temperature, we assume the states
of the bath spins to be random, particularly, in this work, we will use a
completely mixed thermal bath-state $\rho_{B}=\mathbb{1}_{B}/2^{N}$. The so
called gCCE0 contribution to the density matrix reads therefore

$\displaystyle\rho^{(n)}_{\text{E0}}(t_{e})=\mathcal{S}_{n}(\boldsymbol{t},\boldsymbol{\theta},\boldsymbol{\phi})\rho(0)\mathcal{S}_{n}^{\dagger}(\boldsymbol{t},\boldsymbol{\theta},\boldsymbol{\phi}).$
(19)

The unitary operator of the systems evolution within the mean field bath is

$\displaystyle U_{t,n}=\text{exp}\left(-iH^{(E0,n)}t\right),$ (20)

where $n$ labels a specific bath configuration and
$H^{(E0,n)}=H_{R}+H_{B}^{(E0,n)}$. For simulating a completely mixed density
matrix, we need to average over $M$ different random bath-states yielding

$\displaystyle\rho_{\text{E0}}(t)=\frac{1}{M}\sum_{n=1}^{M}\rho^{(n)}_{\text{E0}}(t).$
(21)

The next order (gCCE1) is then calculated by additionally taking the full
interaction between one specific bath spin $l$ and the system into account,
see (13). For this, we start with the initial state

$\displaystyle\tilde{\rho}^{(n,l)}_{\text{E1}}(0)=\rho(0)\otimes\rho^{(n)}_{l}(0),$
(22)

where $n$ denotes the bath configuration and $l$ the spin we add to the
density matrix. The time evolved density matrix of the system for the specific
spin $l$ reads

$\displaystyle\rho^{(n,l)}_{\text{E1}}(t_{e})=\text{tr}_{B,l}\left(\mathcal{S}_{n,l}\,\tilde{\rho}^{(n,l)}_{\text{E1}}(0)\,\mathcal{S}_{n,l}^{\dagger}\right).$
(23)

Note that the pulse sequence operators $\mathcal{S}_{n,l}$ now also contain
the full interaction of the spin $l$ in the environment (encoded in
$H_{B}^{(E1,n,l)}$) with the register. With the definition of these
quantities, a density matrix element of the full gCCE1 contribution is defined
as

Figure 3: (a) and (c) spin-bath configurations for the plots (b)-(d),
respectively. The orange sphere represents the register while the black dots
represent the bath spin. The green lines are guides for the eye. (b) Long time
behavior of the logarithmic negativity for the spin-bath of (a). The gCCE0
contribution sufficiently describes the influence of bath. (d) Long time
behavior of the logarithmic negativity for the spin-bath of (c). The revivals
of the dissipative evolution are correctly captured by the gCCE2
approximation. Here, only the correlation between the encircled pairs of (c)
are taken into account. In both cases the short time behavior of the evolution
is sufficiently described by gCCE0. The gCCE1 contribution coincides with the
gCCE0, since longitudinal relaxation processes are irrelevant on the studied
time scales. Each data point was averaged over 450 different bath-states
(M=450).

$\displaystyle\langle a|\rho_{\text{E1}}^{(n)}(t)|b\rangle=\langle
a|\rho^{(n)}_{\text{E0}}(t)|b\rangle\,\prod_{l}\frac{\langle
a|\rho^{(n,l)}_{\text{E1}}(t)|b\rangle}{\langle
a|\rho^{(n)}_{\text{E0}}(t)|b\rangle},$ (24)

where each numerator in the product contains the full interaction with a
different environmental spin. Since $\rho^{(n,l)}_{\text{E1}}(t)$ contains
also the mean field contribution for all other spins apart from $l$, we have
to divide by it for each term to only multiply corrections to the mean field
approximation. The final gCCE1 density matrix is then calculated via

$\displaystyle\rho_{\text{E1}}(t)=\frac{1}{M}\sum_{n=1}^{M}\rho^{(n)}_{\text{E1}}(t).$
(25)

With this framework also higher order correlation within the spin-bath can be
taken into account, as we show in App. B. In the following, we use this method
to calculate the density matrix of the full spin register influenced by the
environment. To analyze the performance in presence of dissipation, we define
a second fidelity

$\displaystyle\langle\mathcal{F}_{f}\rangle=\frac{\text{tr}\left(\rho_{p}^{\dagger}\rho(t_{e})\right)}{\text{tr}\left(\rho_{p}^{\dagger}\rho_{p}\right)},$
(26)

encoding the deviation from unitary evolution $\rho_{p}$ due to the coupling
with the bath. Note that $\rho(t_{e})$ denotes the exact solution for the
reduced density matrix of the register at the end of the pulse sequence in
presence of the environment and will be approximated by different gCCE orders
in the next sections. Although the calculation of the full fidelity
$\langle\mathcal{F}\rangle=\mathcal{F}_{p}\langle\mathcal{F}_{f}\rangle$ is
straight forward, we will mainly focus on the distinct effect of the spin-bath
and therefore on $\mathcal{F}_{f}$ in the results section. We find that for
short times the gCCE0 contribution, i.e. treating the spin-bath in the mean
field approximation, is sufficient (e.g. $\rho(t)\approx\rho_{\text{E0}}(t)$).
However, we also analyze under which circumstances higher order terms become
necessary in the next section (e.g. $\rho(t)\approx\rho_{\text{E2}}(t)$) (see
also App. B).

##  III Results

###  III.1 Entanglement gate between the NV and one ${}^{13}C$ in presence of
a small spin-bath

Before we discuss the full system, we benchmark the gCCE via exact
diagonalization of a two qubit spin register within a spin-bath consisting of
nine impurities. Thereby, we assume only one ${}^{13}C$ to be close to the
$\text{NV}^{-}$ comparable to the experimental system of Ref. Hegde  _et al._
(2020) and described by the effective system Hamiltonian via

$\\!\\!H_{R}\\!=\\!D_{gs}S_{z}^{2}\\!-\\!\left(g_{e}B_{z}\\!-\\!N_{zz}\right)S_{z}\\!+\\!\boldsymbol{S}\bar{{M}}^{(1)}\boldsymbol{I}^{(1)}_{C}\\!-\\!g_{C}B_{z}I^{(1)}_{C,z},$
(27)

where we assume that the nuclear nitrogen spin remains in the state $m_{N}=1$
and $N_{zz}=-1.76$MHz is its diagonal dipole coupling to the NV electron spin
Körner  _et al._ (2021).

In the following, we simulate an entanglement gate between the NV and
${}^{13}C$ in presence of a small spin-bath. We show the results for the time
evolution of the entanglement process in Fig. 2. In (a) and (b), we display
the evolution of the $x$\- and $z$-component of the nuclear spin due to the
pulse sequence on the NV electron spin (the grey line in each plot shows the
$S_{z}$ component of the NV under influence of the pulse sequence). In (c), we
show the logarithmic negativity ($E_{N}$) as a measure for the entanglement in
the system. This quantity is suitable to quantify the bipartite entanglement
also for mixed density matrices and is defined as

$\displaystyle
E_{N}[{\rho}]=\log_{2}\left(\left\|{\rho}^{T_{A}}\right\|_{1}\right),$ (28)

where $\left\|...\right\|_{1}$ is the trace norm and ${\rho}^{T_{A}}$ the
partially transposed density matrix with respect to one subsystem (here the
${}^{13}C$ nuclear spin)Vidal and Werner (2002); Plenio (2005). We find that
at the end of the pulse sequence ($t_{e}=7.1\,\mu s$ vertical dotted line)
$E_{N}\approx 1$ signaling the success of the process in the case without an
environment (grey dashed line). The unitary evolution shows coherent revivals
of the gate-induced entanglement in the system for $t>t_{e}$. Considering the
bath, we distribute nine impurity spins equally within a spherical shell with
$r_{min}=30$nm and $r_{max}=60\,$nm. This corresponds on average to a volume
density of $\eta_{P1}\approx 65\,$ppb. The red dashed line corresponds to the
result obtained via exact diagonalization. We find that the $I_{C,x}$
component deviates significantly from the unitary evolution. This can be
explained by the fact that the $S_{x}$ component of the NV spin shrinks due to
dephasing making the process on the full register erroneous. Hence, in the
presence of the spin-bath the fidelity of the gate is reduced as can be seen
by the deviation to the unitary evolution at $t=t_{e}$ in Fig.2a-c. For longer
times this deviation increases significantly. To benchmark the above
introduced gCCE, we also show the results obtained by the mean field
contribution (gCCE0) and find that it perfectly describes the exact result
within this time range (blackline). For this simulation, we averaged over
$M=450$ bath-states.

However, on longer time scales higher order terms become important. Especially
for a spin-bath consisting of P1-impurities correlations can build up fast
with respect to the timescales of the process on the register, because the
electron spins are strongly interacting. In Fig. 3, we show the long time
evolution of the logarithmic negativity for two different spin-baths.
Particularly, the configuration in Fig. 3a and b are two examples for a random
distribution of nine P1-impurities within a spherical shell with
$r_{min}=30\,$nm and $r_{max}=90\,$nm corresponding on average to an abundance
$\eta_{P1}\sim 26\,$ppb.

Figure 4: $S_{z}$ component of the NV electron spin under the influence of the
pulse sequences for the different quantum processes on the full register.
Figure 5: Histograms showing the fidelity of different quantum processes as a
function of the abundance of P1-impurities. We simulated 300 different spin-
baths for each volume density. Each spin-bath is averaged over 200 bath-states
($M=200$). (a) GHZ-gate between the $\text{NV}^{-}$ and the two ${}^{13}C$
nuclear spins. (b) Entanglement gate between the ${}^{13}C$ nuclear spins
mediated by the $\text{NV}^{-}$ electron spin.

Fig. 3b shows the evolution of the logarithmic negativity for the spin-bath we
display in (a). Here, we find that the mean field contribution describes the
evolution sufficient also for longer times. This is not the case though for
the spin-bath in Fig. 3c. At a certain point in time the gCCE2 contribution -
also taking correlations within the spin-bath into account - describes the
evolution better then gCCE0 and gCCE1. Particularly, we obtain that the
revivals of the entanglement, which are due to resonances with the spin-bath,
are only correctly simulated when we account for the growth of correlations in
the environment. For this specific case it is enough to include the
correlations between the closest pairs of bath spins (highlighted in Fig. 3c
by the circles), since correlation between spins that are further away builds
up on longer time scales. We performed this selection by only taking the gCCE2
contribution of spins that are closer than $70\,$nm and have a second spin
within $60\,$nm in their surrounding. Reversing the discussion - for a small
number of spins in the environment - the gCCE can also be used as a tool to
analyze a given spin-bath and to find which specific correlations become
important for the register on different timescales. Generally, for for longer
times and the mentioned volume densities also higher order corrections (gCCE3
and higher) become important. We leave this analysis for future work.

For larger spin-baths the exact diagonalization is not possible any more,
since the Hilbert space grows exponentially. Without the possibility to
exactly diagonalize the full system for comparison, the gCCE has to be
performed until convergence. In accordance with previous findings, we observe
the mean field term to be a good approximation only on short timescales
$t<25\,\mu s$, densities of $\eta<65\,\text{ppb}$, and $r_{\text{min}}=30$nm.
Specifically, its performance compared to the gCCE2 accounting for spin-spin
interactions is shown for 100 different spin-baths in App. B.

###  III.2 Four-qubit quantum processes in presence of dephasing

In the following, we study a four-spin-register (one NV electron spin + two
$\times$ ${}^{13}C$ \+ one Nitrogen nuclear spin, see Fig. 1) in a spin-bath
containing P1-impurities. We randomly distribute the bath spins within a
spherical shell with $r_{\text{min}}=30\,$nm and $r_{\text{max}}=250\,$nm and
study different volume densities (for a volume density of $45$ppb this
corresponds to 518 bath spins). We assume that the register is initialized in
the state
$|\psi_{0}\rangle=|0\rangle_{NV}\\!\otimes\\!|\frac{1}{2}\rangle_{{}^{13}C}\\!\otimes\\!|\frac{1}{2}\rangle_{{}^{13}C}\\!\otimes\\!|1\rangle_{N}$
and study two different quantum processes that are again solely initiated by
pulsing the NV spin. The first pulse sequence yields a GHZ state between the
NV electron spin and both ${}^{13}C$ nuclear spins. The nitrogen nuclear spin
is supposed to stay in the $|1\rangle$ state. The process can be expressed as

$\displaystyle|0,\frac{1}{2},\frac{1}{2},1\rangle\xrightarrow{\mathcal{S}_{\text{GHZ}}(\boldsymbol{t},\boldsymbol{\phi},\boldsymbol{\theta})}|\text{GHZ}\rangle\\!\otimes\\!|1\rangle_{N},$
(29)

where
$|\text{GHZ}\rangle=\frac{1}{\sqrt{2}}\left(|0,\frac{1}{2},\frac{1}{2}\rangle+|-1,-\frac{1}{2},-\frac{1}{2}\rangle\right)$.
The optimization algorithm yields the pulse structure displayed via the
$S_{z}$ component in Fig. 4. Since the register is fully connected (general
spin-spin interactions), even in the non dissipative case we find a process
fidelity $\mathcal{F}_{p}<1$. For the GHZ-gate entangling both ${}^{13}C$
nuclear spins and the NV electron spin, we find $\mathcal{F}_{p}\approx
0.9954$ and a duration of $t_{e}\approx 21\,\mu s$. This control sequence
might be further optimized using more sophisticated algorithms. We leave the
optimization of the processes for future work and solely analyze the influence
of the spin-bath. The above analysis shows, that for the spin-baths we
consider the pure dephasing approximation is valid within these timescales and
the mean field term of the gCCE gives satisfactory results. To understand the
influence of the volume density as well as the different configurations of
impurities better, we simulated 300 different baths for each volume density
and plot a histogram of the corresponding fidelities
$\langle\mathcal{F}_{f}\rangle$ in Fig. 5a. There, the tallest bar shows the
probability to have a fidelity $\langle\mathcal{F}_{f}\rangle>0.997$ to be
above $80\%$ for a impurity density of $\eta_{P1}=5\,\text{ppb}$. On the
x-axis, we increase the impurity density subsequently by $5\text{ppb}$ up to
$50\,\text{ppb}$ and find that the distribution broadens as expected. This
broadening yields insight into the most probable magnitude of errors induced
by P1-impurities. The second process we discuss is an entanglement gate
between the two ${}^{13}C$s mediated by the only driving the NV electron spin,
hence

$\displaystyle|0,\frac{1}{2},\frac{1}{2},1\rangle\xrightarrow{\mathcal{S}_{\text{Bell}}(\boldsymbol{t},\boldsymbol{\phi},\boldsymbol{\theta})}|0\rangle_{NV}\\!\otimes\\!|\text{Bell}\rangle_{{}^{13}C}\\!\otimes\\!|1\rangle_{N},$
(30)

where
$|\text{Bell}\rangle_{{}^{13}C}=\frac{1}{\sqrt{2}}\left(|\frac{1}{2},\frac{1}{2}\rangle+|-\frac{1}{2},-\frac{1}{2}\rangle\right)$.
The corresponding gate sequence is also displayed in Fig. 4 yielding a process
fideltiy of $\mathcal{F}_{p}\approx 0.9968$. The errors originating from the
spin-bath are summarized in the histogram Fig. 5b. We find that the overall
fidelity is increased for this process compared to the GHZ-gate. This can be
understood by the fact that the electron spin of the NV - which is more
strongly coupled to the environment - is not in a coherent state at the end of
the gate and the state therefore less affected by decoherence. This has
naturally a also crucial influence on the protection of the target state.
While the evolution after the GHZ-gate is strongly influenced by the spin-bath
the evolution after the entanglement gate between the ${}^{13}C$ nuclear spins
remains almost unitary for a long time making the the ${}^{13}C$’s a suitable
quantum memory. However, note that in both processes the histogram falls off
quite rapidly with the enhancement of the volume density of the P1-impurities.
The probability to have a fidelity $\langle\mathcal{F}_{f}\rangle>0.997$ is
already below $80\,\%$ for sparse spin environments with volume density of
$\eta=10\,$ppb. We stress that in a real diamond crystal the P1-impurities are
only one of many sources of errors. In this study, we are still neglecting
errors in the initialization, in the gate sequences and in the readout. Also
the ${}^{13}C$ spin-bath is not considered. Our study nevertheless describes a
’best case’ scenario and yields an upper limit for the fidelity for processes
without sophisticated dynamical decoupling sequences and without error
correction.

##  IV Discussion and Outlook

In summary, we theoretically studied the performance of a quantum register in
diamond in presence of a spin-bath. The register is build of one
$\text{NV}^{-}$ and one or two close by ${}^{13}C$ spins. We initiate the
processes by only pulsing the NV electron spin and studied the corresponding
dissipative time evolution of the full system. To model the spin-bath, we used
the generalized cluster correlation expansion (Onizhuk  _et al._ (2021)) for
the full register. We benchmarked this technique via an exact diagonalization
of a small spin register within a spin-bath with nine impurities and found
that for short times $<25\,\mu$s and P1 abundances $\eta_{P1}$ of several tens
of ppb the mean field contribution (gCCE0) is dominating the dissipative
evolution. Further, we showed that in an environment consisting of P1-impurity
spins the correlations within the spin-bath can become important on timescales
of several tens of microseconds.

For practical purposes, we found that only taking into account correlations
between closest pairs of spins can already give the desired correction to the
gCCE1 contribution. While keeping the accuracy of a full gCCE2-simulation, it
substantially speeds up the computation. For a four spin register, we analyzed
two different quantum gates and their fidelities with experimentally realistic
parameters and depending on the abundance of P1-impurities. This yields
quantitative predictions how gate-performances deteriorate with increasing
P1-densities and thus provides thresholds to be met by future material design
and pulse shaping.

The methodology presented here, now allows to easily implement sophisticated
decoupling sequences or error correction schemes for consecutive gate
operations Zhang  _et al._ (2014); Zhang and Suter (2015); Kern and Alber
(2005); Ng  _et al._ (2011); Choi  _et al._ (2020); Bradley  _et al._ (2021);
Waldherr  _et al._ (2014). Our findings thus contribute to the rising
activities in the less explored field of how to protect multi NV-based quantum
gates from surrounding degrees of freedom. Since they are based on a framework
treating the complete density matrix of the entire register, they offer a
realistic test-bed for future concepts with quantitative predictions.

###### Acknowledgements.

We are grateful to R. Ghassemizadeh, W. Hahn, D. Hähnel, J. Stockburger and B.
Donvil for fruitful discussions. This work has been supported by the State of
Baden - Württemberg under KQCBW/SiQuRe.

Figure 6:  Histograms for the numerical error of the mean field (gCCE0) and the gCCE2 approximation with respect to the exact diagonalization for the fidelity after a dissipative evolution of $t=20\,\mu$s in (a) and for $t=65\,\mu$s in (b). The histograms are calculated by sampling over 100 different environments for $\eta_{\text{P1}}\approx 65\,$ppb in (a) and $\eta_{\text{P1}}\approx 26\,$ppb in (b). The error is less than one percent on average in both approximation after averaging over $M=200$ states in (a), making the results in Fig. 5 accurate. The situation in (b) is the same as in Fig. 3 ($\eta_{\text{P1}}\approx 26\,$ppb). The spin baths for which correlation become important are highlighted and have a significantly reduced error in the gCCE2 approximation.  | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---|---|---|---
Bell state: two spin register |  |  |  |  |  |  |  |  |
$\boldsymbol{t}[\mu\text{s}]$ | $0$ | $4.061$ | $1.763$ | $1.276$ | $0$ | $-$ | $-$ | $-$ | $-$
$\boldsymbol{\theta}[\text{rad}]$ | $-$ | $1.059$ | $3.566$ | $1.627$ | $3.360$ | $-$ | $-$ | $-$ | $-$
$\boldsymbol{\phi}[\text{rad}]$ | $-$ | $0.669$ | $1.952$ | $0.255$ | $0.503$ | $-$ | $-$ | $-$ | $-$
GHZ-state |  |  |  |  |  |  |  |  |
$\boldsymbol{t}[\mu\text{s}]$ | $0.367$ | $3.648$ | $2.482$ | $2.115$ | $2.675$ | $1.296$ | $2.653$ | $3.928$ | $2.059$
$\boldsymbol{\theta}[\text{rad}]$ | $-$ | $5.052$ | $0.146$ | $3.980$ | $1.118$ | $0.415$ | $5.398$ | $1.622$ | $1.990$
$\boldsymbol{\phi}[\text{rad}]$ | $-$ | $6.326$ | $5.465$ | $1.965$ | $3.141$ | $5.772$ | $3.330$ | $3.212$ | $2.983$
Bell state ${}^{13}C$ |  |  |  |  |  |  |  |  |
$\boldsymbol{t}[\mu\text{s}]$ | $0.4028$ | $3.492$ | $2.634$ | $2.695$ | $2.531$ | $0.854$ | $2.620$ | $3.863$ | $2.960$
$\boldsymbol{\theta}[\text{rad}]$ | $-$ | $5.369$ | $0.515$ | $5.422$ | $1.750$ | $2.332$ | $4.820$ | $1.695$ | $1.527$
$\boldsymbol{\phi}[\text{rad}]$ | $-$ | $6.716$ | $4.788$ | $2.692$ | $2.274$ | $5.522$ | $2.632$ | $3.376$ | $2.554$
Table 1: Pulse sequences for the processes simulated in the main text as
described in Fig. 1b.

##  Appendix A Interaction tensors

To make the specific nature of the simulated register more compareble to
experiments, we show the dipole interaction tensors of the fully connected
register in this Appendix. We start with $\overline{M}$ encoding the
interaction between the $\text{NV}^{-}$ and the ${}^{13}C$ nuclear spins

$\displaystyle\\!\\!\boldsymbol{S}\bar{M}^{(k)}\boldsymbol{I}^{(k)}_{C}\\!=\\!-\frac{\hbar\mu_{0}}{4\pi}\frac{\gamma_{e}\gamma_{C}}{r^{3}}[3(\boldsymbol{S}\\!\cdot\\!\hat{n})(\boldsymbol{I}^{(k)}_{C}\\!\cdot\\!\hat{n})\\!-\\!\boldsymbol{S}\\!\cdot\\!\boldsymbol{I}^{(k)}_{C}],$
(31)

where $r$ is the distance of the ${}^{13}C$ to the $\text{NV}^{-}$ electron
spin, $\hat{n}$ the unit vector pointing to the ${}^{13}C$,
$\gamma_{e}=-1.761\times 10^{11}\,T^{-1}s^{-1}$ and $\gamma_{C}=6.728\times
10^{7}\,T^{-1}s^{-1}$ are the gyromagnetic ratios of the electron and the
nuclear spin, respectively. The first ${}^{13}C$ (k=1) has the spherical
coordinates
$\boldsymbol{r}=(r=0.89\,\text{nm},\varphi=0^{\circ},\vartheta=78^{\circ})$
yielding

$\displaystyle\bar{M}^{(1)}_{xx}$ $\displaystyle\approx 310.8\text{kHz},$ (32)
$\displaystyle\bar{M}^{(1)}_{yy}$ $\displaystyle\approx 166.2\text{kHz},$
$\displaystyle\bar{M}^{(1)}_{zz}$ $\displaystyle\approx-144.6\text{kHz},$
$\displaystyle\bar{M}^{(1)}_{xy}$
$\displaystyle=\bar{M}^{(1)}_{yx}=\bar{M}^{(1)}_{yz}=\bar{M}^{(1)}_{zy}=0,$
$\displaystyle\bar{M}^{(1)}_{zx}$ $\displaystyle=\bar{M}^{(1)}_{xz}\approx
101.4\text{kHz}.$

Analogously, we find for the second ${}^{13}C$ (k=2) having the spherical
coordinates
$\boldsymbol{r}=(r=1.0\text{nm},\varphi=53.97^{\circ},\vartheta=107.94^{\circ})$

$\displaystyle\bar{M}^{(2)}_{xx}$ $\displaystyle\approx-7.381\text{kHz},$ (33)
$\displaystyle\bar{M}^{(2)}_{yy}$ $\displaystyle\approx 90.16\text{kHz},$
$\displaystyle\bar{M}^{(2)}_{zz}$ $\displaystyle\approx-84.26\text{kHz},$
$\displaystyle\bar{M}^{(2)}_{xy}$ $\displaystyle=\bar{M}^{(2)}_{yx}\approx
152.4\text{kHz},$ $\displaystyle\bar{M}^{(2)}_{yz}$
$\displaystyle=\bar{M}^{(2)}_{zy}\approx-84.23\text{kHz},$
$\displaystyle\bar{M}^{(2)}_{zx}$
$\displaystyle=\bar{M}^{(2)}_{xz}\approx-61.2\text{kHz}.$

The interaction tensors between the two ${}^{13}C$ between the ${}^{13}C$s and
the nitrogen nuclear spin do not play a role in this specific register since
the interactions are of the order of several Hz. However, we included it into
our simulation for future analysis. Since the interaction between the nuclear
spin nitrogen atom and the electron spin of the NV is not solely defined by
the dipolar interaction but also by the Fermi-Contact contribution, we use
$N_{\parallel}=-1.73\,$MHz and $N_{\perp}=-2.16\,$MHz as given in Körner  _et
al._ (2021).

##  Appendix B gCCE: Accuracy and higher orders

###  B.1 Accuracy of the mean field approximation (gCCE0)

As mentioned in the main text, the gCCE0 contribution coincides with the mean
field approximation. We only take this contribution for the histograms in Fig.
5. As a justification of this procedure, we benchmark the gCCE0 and gCCE2
against the exact diagonalization for 100 different spin-baths. Each of the
environments consists out of nine impurity spins. First, we distribute them
within a spherical shell with $r_{\text{min}}=30\,$nm and
$r_{\text{max}}=60\,$nm. For each bath configuration, we calculate the
fidelity ($\mathcal{F}_{f}$) for the two qubit register - with respect to the
unitary evolution - after the entangling gate at $t=20\,\mu$s. We calculate
the relative numerical error with respect to the exact solution. We show the
distribution of relative errors in Fig. 6a. Note that, on average, the here
simulated spin-baths are more dense as the ones simulated in Fig. 5, making
the results of Sec. III.2 very accurate. The average error is $\sim 0.4\%$ for
gCCE0 and gCCE2. In Fig. 6b, we analyze the effect of the gCCE2 contributions
for different spin-baths with abundances comparable to Fig. 3
($\eta_{\text{P1}}\approx 26\,$ppb). Here, we simulate longer timescales to
make the effect of correlation in environment visible. We find cases with a
large numerical error in the gCCE0 approximation which are cured by the gCCE2
contribution. The remaining error in the gCCE2 case (on average $\sim 1\,\%$),
can be attributed to the averaging over different bath-states. Note that we
included correlations between only for spins that are closer then $60\,$nm and
have a second spin within $55\,$nm in their surrounding.

###  B.2 Higher orders

Here we briefly dicuss the calculation of higher orders in the gCCE. For the
gCCE2 contribution, we start with the density matrix

$\displaystyle\tilde{\rho}^{(n,l,q)}_{\text{E2}}(0)=\rho(0)\otimes\rho^{(n)}_{l,q}(0),$
(34)

where we included two spins ($l$ and $q$) of the environment into the full
density matrix. The reduced density matrix then becomes

$\displaystyle\rho^{(n,l,q)}_{\text{E2}}(t_{e})=\text{tr}_{B,l,q}\left(\mathcal{S}_{n,l,q}\,\tilde{\rho}^{(n,l,q)}_{\text{E2}}(0)\,\mathcal{S}_{n,l,q}^{\dagger}\right).$
(35)

To obtain the correction to the gCCE1, we use

$\displaystyle\langle a|\rho_{\text{E2}}^{(n)}(t)|b\rangle$
$\displaystyle=\langle
a|\rho^{(n)}_{\text{E0}}(t)|b\rangle\,\left[\prod_{l}\frac{\langle
a|\rho^{(n,l)}_{\text{E1}}(t)|b\rangle}{\langle
a|\rho^{(n)}_{\text{E0}}(t)|b\rangle}\right]$ (36)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\times\prod_{l,q}\frac{\langle
a|\rho^{(n,l,q)}_{\text{E2}}(t)|b\rangle}{\langle
a|\rho^{(n)}_{\text{E0}}(t)|b\rangle^{-1}\langle
a|\rho^{(n,l)}_{\text{E1}}(t)|b\rangle\langle
a|\rho^{(n,q)}_{\text{E1}}(t)|b\rangle}.$

Averaging over different bath-states $n$ yields the gCCE density matrix to
second order. Generally j-th order correction of the gCCE can be recursively
written as

$\displaystyle\langle a|\tilde{\rho}^{(n)}_{Ej}|b\rangle=\frac{\left\langle
a\left|\rho_{Ej}^{(C_{j})}\right|b\right\rangle}{\prod_{{C}^{\prime}}\langle
a|\tilde{\rho}_{\left\\{{C}^{\prime}\subset C_{j}\right\\}}|b\rangle},$ (37)

where $C_{j}$ contains a number of indices coinciding with number of included
bath spins $j$ (e.g. $C_{2}\rightarrow l,p$ and $C_{3}\rightarrow l,p,s$).
Using this definition a general order $j$ reads

$\displaystyle\langle a|\rho_{\text{Ej}}^{(n)}|b\rangle$
$\displaystyle=\langle
a|\rho^{(n)}_{\text{E0}}|b\rangle\,\left[\prod_{C_{1}}\langle
a|\tilde{\rho}^{(n)}_{E1}|b\rangle\right]$ (38)
$\displaystyle\times\left[\prod_{C_{2}}\langle
a|\tilde{\rho}^{(n)}_{E2}|b\rangle\right]...\times\left[\prod_{C_{j}}\langle
a|\tilde{\rho}^{(n)}_{Ej}|b\rangle\right].$

This result is again to be averaged over different states $n$. For more
information we refer to Onizhuk  _et al._ (2021).

##  Appendix C Tables of pulse sequences

We display the parameters for the pulse sequences (as explained in Fig. 1b)
used to initiate the quantum processes on the register in table 1.

## References

  * Doherty  _et al._ (2013) M. W. Doherty, N. B. Manson, P. Delaney, F. Jelezko, J. Wrachtrup, and L. C. Hollenberg,  _The nitrogen-vacancy colour centre in diamond,_ Physics Reports 528, 1 (2013) The nitrogen-vacancy colour centre in diamond.
  * Pezzagna and Meijer (2021) S. Pezzagna and J. Meijer,  _Quantum computer based on color centers in diamond,_ Applied Physics Reviews 8, 011308 (2021) .
  * Wolfowicz  _et al._ (2021) G. Wolfowicz, F. J. Heremans, C. P. Anderson, S. Kanai, H. Seo, A. Gali, G. Galli, and D. D. Awschalom,  _Quantum guidelines for solid-state spin defects,_ Nature Reviews Materials 6, 906 (2021) .
  * Jelezko  _et al._ (2004a) F. Jelezko, T. Gaebel, I. Popa, A. Gruber, and J. Wrachtrup,  _Observation of Coherent Oscillations in a Single Electron Spin,_ Phys. Rev. Lett. 92, 076401 (2004a) .
  * Dolde  _et al._ (2013) F. Dolde, I. Jakobi, B. Naydenov, N. Zhao, S. Pezzagna, C. Trautmann, J. Meijer, P. Neumann, F. Jelezko, and J. Wrachtrup,  _Room-temperature entanglement between single defect spins in diamond,_ Nature Physics 9, 139 (2013) .
  * Haruyama  _et al._ (2019) M. Haruyama, S. Onoda, T. Higuchi, W. Kada, A. Chiba, Y. Hirano, T. Teraji, R. Igarashi, S. Kawai, H. Kawarada, Y. Ishii, R. Fukuda, T. Tanii, J. Isoya, T. Ohshima, and O. Hanaizumi,  _Triple nitrogen-vacancy centre fabrication by C5N4Hn ion implantation,_ Nature Communications 10, 2664 (2019) .
  * Jakobi  _et al._ (2016) I. Jakobi, S. A. Momenzadeh, F. F. de Oliveira, J. Michl, F. Ziem, M. Schreck, P. Neumann, A. Denisenko, and J. Wrachtrup,  _Efficient creation of dipolar coupled nitrogen-vacancy spin qubits in diamond,_ Journal of Physics: Conference Series 752, 012001 (2016) .
  * Groot-Berning  _et al._ (2021) K. Groot-Berning, G. Jacob, C. Osterkamp, F. Jelezko, and F. Schmidt-Kaler, New Journal of Physics 23, 063067 (2021) .
  * Yamamoto  _et al._ (2013) T. Yamamoto, C. Müller, L. P. McGuinness, T. Teraji, B. Naydenov, S. Onoda, T. Ohshima, J. Wrachtrup, F. Jelezko, and J. Isoya,  _Strongly coupled diamond spin qubits by molecular nitrogen implantation,_ Phys. Rev. B 88, 201201 (2013) .
  * Bauch  _et al._ (2020) E. Bauch, S. Singh, J. Lee, C. A. Hart, J. M. Schloss, M. J. Turner, J. F. Barry, L. M. Pham, N. Bar-Gill, S. F. Yelin, and R. L. Walsworth,  _Decoherence of ensembles of nitrogen-vacancy centers in diamond,_ Phys. Rev. B 102, 134210 (2020) .
  * Jelezko  _et al._ (2004b) F. Jelezko, T. Gaebel, I. Popa, M. Domhan, A. Gruber, and J. Wrachtrup,  _Observation of Coherent Oscillation of a Single Nuclear Spin and Realization of a Two-Qubit Conditional Quantum Gate,_ Phys. Rev. Lett. 93, 130501 (2004b) .
  * Childress  _et al._ (2006) L. Childress, M. V. G. Dutt, J. M. Taylor, A. S. Zibrov, F. Jelezko, J. Wrachtrup, P. R. Hemmer, and M. D. Lukin,  _Coherent Dynamics of Coupled Electron and Nuclear Spin Qubits in Diamond,_ Science 314, 281 (2006) .
  * Neumann  _et al._ (2008) P. Neumann, N. Mizuochi, F. Rempp, P. Hemmer, H. Watanabe, S. Yamasaki, V. Jacques, T. Gaebel, F. Jelezko, and J. Wrachtrup,  _Multipartite Entanglement Among Single Spins in Diamond,_ Science 320, 1326 (2008) .
  * Jacques  _et al._ (2009) V. Jacques, P. Neumann, J. Beck, M. Markham, D. Twitchen, J. Meijer, F. Kaiser, G. Balasubramanian, F. Jelezko, and J. Wrachtrup,  _Dynamic Polarization of Single Nuclear Spins by Optical Pumping of Nitrogen-Vacancy Color Centers in Diamond at Room Temperature,_ Phys. Rev. Lett. 102, 057403 (2009) .
  * Bernien  _et al._ (2013) H. Bernien, B. Hensen, W. Pfaff, G. Koolstra, M. S. Blok, L. Robledo, T. H. Taminiau, M. Markham, D. J. Twitchen, L. Childress, and R. Hanson,  _Heralded entanglement between solid-state qubits separated by three metres,_ Nature 497, 86 (2013) .
  * Hensen  _et al._ (2015) B. Hensen, H. Bernien, A. E. Dréau, A. Reiserer, N. Kalb, M. S. Blok, J. Ruitenberg, R. F. L. Vermeulen, R. N. Schouten, C. Abellán, W. Amaya, V. Pruneri, M. W. Mitchell, M. Markham, D. J. Twitchen, D. Elkouss, S. Wehner, T. H. Taminiau, and R. Hanson,  _Loophole-free Bell inequality violation using electron spins separated by 1.3 kilometres,_ Nature 526, 682 (2015) .
  * Bradley  _et al._ (2021) C. E. Bradley, S. W. de Bone, P. F. W. Moller, S. Baier, M. J. Degen, S. J. H. Loenen, H. P. Bartling, M. Markham, D. J. Twitchen, R. Hanson, D. Elkouss, and T. H. Taminiau,  _Robust quantum-network memory based on spin qubits in isotopically engineered diamond,_ (2021), 10.48550/ARXIV.2111.09772 .
  * Pompili  _et al._ (2021) M. Pompili, S. L. N. Hermans, S. Baier, H. K. C. Beukers, P. C. Humphreys, R. N. Schouten, R. F. L. Vermeulen, M. J. Tiggelman, L. dos Santos Martins, B. Dirkse, S. Wehner, and R. Hanson,  _Realization of a multinode quantum network of remote solid-state qubits,_ Science 372, 259 (2021) .
  * Hegde  _et al._ (2020) S. S. Hegde, J. Zhang, and D. Suter,  _Efficient Quantum Gates for Individual Nuclear Spin Qubits by Indirect Control,_ Phys. Rev. Lett. 124, 220501 (2020) .
  * Zhang  _et al._ (2020) J. Zhang, S. S. Hegde, and D. Suter,  _Efficient Implementation of a Quantum Algorithm in a Single Nitrogen-Vacancy Center of Diamond,_ Phys. Rev. Lett. 125, 030501 (2020) .
  * Bradley  _et al._ (2019) C. E. Bradley, J. Randall, M. H. Abobeih, R. C. Berrevoets, M. J. Degen, M. A. Bakker, M. Markham, D. J. Twitchen, and T. H. Taminiau,  _A Ten-Qubit Solid-State Spin Register with Quantum Memory up to One Minute,_ Phys. Rev. X 9, 031045 (2019) .
  * Bar-Gill  _et al._ (2012) N. Bar-Gill, L. M. Pham, C. Belthangady, D. Le Sage, P. Cappellaro, J. R. Maze, M. D. Lukin, A. Yacoby, and R. Walsworth,  _Suppression of spin-bath dynamics for improved coherence of multi-spin-qubit systems,_ Nature Communications 3, 858 (2012) .
  * Bar-Gill  _et al._ (2013) N. Bar-Gill, L. M. Pham, A. Jarmola, D. Budker, and R. L. Walsworth,  _Solid-state electronic spin coherence time approaching one second,_ Nature Communications 4, 1743 (2013) .
  * de Sousa and Das Sarma (2003) R. de Sousa and S. Das Sarma,  _Theory of nuclear-induced spectral diffusion: Spin decoherence of phosphorus donors in Si and GaAs quantum dots,_ Phys. Rev. B 68, 115322 (2003) .
  * Yao  _et al._ (2006) W. Yao, R.-B. Liu, and L. J. Sham,  _Theory of electron spin decoherence by interacting nuclear spins in a quantum dot,_ Phys. Rev. B 74, 195301 (2006) .
  * Witzel and Das Sarma (2006) W. M. Witzel and S. Das Sarma,  _Quantum theory for electron spin decoherence induced by nuclear spin dynamics in semiconductor quantum computer architectures: Spectral diffusion of localized electron spins in the nuclear solid-state environment,_ Phys. Rev. B 74, 035322 (2006) .
  * Breuer (2004) H.-P. Breuer,  _Exact quantum jump approach to open systems in bosonic and spin baths,_ Phys. Rev. A 69, 022115 (2004) .
  * Dobrovitski  _et al._ (2008) V. V. Dobrovitski, A. E. Feiguin, D. D. Awschalom, and R. Hanson,  _Decoherence dynamics of a single spin versus spin ensemble,_ Phys. Rev. B 77, 245212 (2008) .
  * Zhao  _et al._ (2012) N. Zhao, S.-W. Ho, and R.-B. Liu,  _Decoherence and dynamical decoupling control of nitrogen vacancy center electron spins in nuclear spin baths,_ Phys. Rev. B 85, 115303 (2012) .
  * Yang and Liu (2008) W. Yang and R.-B. Liu,  _Quantum many-body theory of qubit decoherence in a finite-size spin bath,_ Phys. Rev. B 78, 085315 (2008) .
  * Hall  _et al._ (2014) L. T. Hall, J. H. Cole, and L. C. L. Hollenberg,  _Analytic solutions to the central-spin problem for nitrogen-vacancy centers in diamond,_ Phys. Rev. B 90, 075201 (2014) .
  * Myers  _et al._ (2017) B. A. Myers, A. Ariyaratne, and A. C. B. Jayich,  _Double-Quantum Spin-Relaxation Limits to Coherence of Near-Surface Nitrogen-Vacancy Centers,_ Phys. Rev. Lett. 118, 197201 (2017) .
  * Onizhuk  _et al._ (2021) M. Onizhuk, K. C. Miao, J. P. Blanton, H. Ma, C. P. Anderson, A. Bourassa, D. D. Awschalom, and G. Galli,  _Probing the Coherence of Solid-State Qubits at Avoided Crossings,_ PRX Quantum 2, 010311 (2021) .
  * Kwiatkowski and Cywiński (2018) D. Kwiatkowski and L. Cywiński,  _Decoherence of two entangled spin qubits coupled to an interacting sparse nuclear spin bath: Application to nitrogen vacancy centers,_ Phys. Rev. B 98, 155202 (2018) .
  * Bragar and Cywiński (2015) I. Bragar and L. Cywiński,  _Dynamics of entanglement of two electron spins interacting with nuclear spin baths in quantum dots,_ Phys. Rev. B 91, 155310 (2015) .
  * Ivády (2020) V. Ivády,  _Longitudinal spin relaxation model applied to point-defect qubit systems,_ Phys. Rev. B 101, 155203 (2020) .
  * Norambuena  _et al._ (2018) A. Norambuena, E. Muñoz, H. T. Dinani, A. Jarmola, P. Maletinsky, D. Budker, and J. R. Maze,  _Spin-lattice relaxation of individual solid-state spins,_ Phys. Rev. B 97, 094304 (2018) .
  * Cambria  _et al._ (2021) M. C. Cambria, A. Gardill, Y. Li, A. Norambuena, J. R. Maze, and S. Kolkowitz,  _State-dependent phonon-limited spin relaxation of nitrogen-vacancy centers,_ Phys. Rev. Research 3, 013123 (2021) .
  * Note (1) In natural diamond the abundance $\eta_{{}^{13}C}=1.1\%$ yielding to a reasonable probability to find strongly coupled ${}^{13}C$’s ($\sim 100\text{MHz}$). Purified diamond (e.g. $\eta_{{}^{13}C}\approx 0.01\%$) leads to a sparse ${}^{13}C$ bath with better coherence properties for the NV spin.
  * Körner  _et al._ (2021) W. Körner, D. F. Urban, and C. Elsässer,  _Influence of extended defects on the formation energy, hyperfine structure, and zero-field splitting of NV centers in diamond,_ Phys. Rev. B 103, 085305 (2021) .
  * Vidal and Werner (2002) G. Vidal and R. F. Werner,  _Computable measure of entanglement,_ Phys. Rev. A 65, 032314 (2002) .
  * Plenio (2005) M. B. Plenio,  _Logarithmic Negativity: A Full Entanglement Monotone That is not Convex,_ Phys. Rev. Lett. 95, 090503 (2005) .
  * Zhang  _et al._ (2014) J. Zhang, A. M. Souza, F. D. Brandao, and D. Suter,  _Protected Quantum Computing: Interleaving Gate Operations with Dynamical Decoupling Sequences,_ Phys. Rev. Lett. 112, 050502 (2014) .
  * Zhang and Suter (2015) J. Zhang and D. Suter,  _Experimental Protection of Two-Qubit Quantum Gates against Environmental Noise by Dynamical Decoupling,_ Phys. Rev. Lett. 115, 110502 (2015) .
  * Kern and Alber (2005) O. Kern and G. Alber,  _Controlling Quantum Systems by Embedded Dynamical Decoupling Schemes,_ Phys. Rev. Lett. 95, 250501 (2005) .
  * Ng  _et al._ (2011) H. K. Ng, D. A. Lidar, and J. Preskill,  _Combining dynamical decoupling with fault-tolerant quantum computation,_ Phys. Rev. A 84, 012305 (2011) .
  * Choi  _et al._ (2020) J. Choi, H. Zhou, H. S. Knowles, R. Landig, S. Choi, and M. D. Lukin,  _Robust Dynamic Hamiltonian Engineering of Many-Body Spin Systems,_ Phys. Rev. X 10, 031002 (2020) .
  * Waldherr  _et al._ (2014) G. Waldherr, Y. Wang, S. Zaiser, M. Jamali, T. Schulte-Herbrüggen, H. Abe, T. Ohshima, J. Isoya, J. F. Du, P. Neumann, and J. Wrachtrup,  _Quantum error correction in a solid-state hybrid spin register,_ Nature 506, 204 (2014) .


# Feynman diagram approach to dynamical Casimir effect in optimechanical
cavity

Yu-Song Cao111caoyusong15@mails.ucas.ac.cn$~{}^{1,2}$,YanXia Liu$~{}^{2}$

###### Abstract

In this paper we study an optomechnical system enclosed by an optical cavity
with one mirror attached to a spring as a closed quantum system. We provide a
different angle of studying the phenomenons related to the dynamical Casimir
effect via Feynman diagram technique. Dressing effects of phonon, photon and
coupling strength are discussed. The energy shift of ground state is obtained.
The dynamical Casimir effect is modeled by the scattering processes converting
phonons to photons and the corresponding scattering amplitudes are computed.
The force-force correlation function of the radiation pressure is derived,
whose non-Gaussian probability distribution is revealed.

$~{}^{1}$School of Physical Sciences, University of Chinese Academy of
Sciences, Beijing, 100049, China

$~{}^{2}$Beijing National Laboratory for Condensed Matter Physics, Institute
of Physics, Chinese Academy of Sciences, Beijing 100190, China

##  1 Introduction

The magic lies in the center of quantum mechanics is the uncertainty
principle, as a result of which, the quantum systems exhibit non-vanishing
zero-point energy. When the degrees of freedom of the system goes infinite, we
encounter the first intrinsic divergence in quantum field theory: the absolute
value of vacuum energy [1]. This crisis is solved by introducing the normal
ordering technique, which indicates the absolute value of vacuum energy has no
physical meaning because it is unobservable. However, it is not the case for
the real world. In gravitational systems, the Einstein tensor is coupled to
the energy momentum tensor thus the vacuum energy gives rise to the
accelerating expansion of the universe [2]. The vacuum energy shows its
presence also in the flat timespace, it is discovered that if the quantum
vacuum is disturbed non-adiabatically, particles can emerge from the vacuum,
which gives rise to the astonishing phenomenons such as dynamical Casimir
effect [3, 4, 5, 6], Hawking radiation [7, 8, 9] and spontaneous radiation
[1]. So far, the dynamical Casimir effect has been verified in many
experimental systems [10, 11]. It is notable that the dynamical Casimir effect
serves as an analog of studying the Unruh effect [12] and Hawking radiation
[3, 13, 14, 15].

Many theoretical treatments of the dynamical Casimir effect are developed
semiclassically, where the motion of the boundary is given classically prior
and the field dynamics is described quantized field equation [16, 17, 18, 19,
20, 21]. Although the semiclassical theory proved its success, however, people
still want a fully quantized description of dynamical Casimir effect. The
first reason is philosophical consistence while the second reason is more
realistic. In systems requiring the high accuracy of the measurement of the
mirror’s position, the quantum fluctuation must be taken into consideration
[22]. Moreover, the semiclassical theory has an Achillies’ heel of failing to
give an integrated description of dynamical properties of the boundary. In
semiclassical theory, the backaction force quantum field exerts on the mirror
is given in an average level [17, 18, 19, 20, 21, 23].

The fully quantum description of the dynamical Casimir effect is made by Law
[24]. After which various new effects as a result of the quantum nature of the
mirror have been discovered [25, 26, 27, 28, 29, 30, 31, 32]. The backaction
effect, is of course studied. In [25, 26], the dynamical Casimir effect is
modeled by the evolution of the mean value of the field operator and
excitation number of the harmonic oscillator and the dispassion to the
enviorment is considered crucial. In this paper, we study the behavior of
Hamiltonian in [24] with Feynman diagram approach. The Feynman rules with
single mode approximation is obtained, which will show its power in studying
the dynamical effects of the system. With Feynman diagrams we study the
dressing effects of the photon, phonon and the coupling strength. We give a
picturise description on how the virtual particles presents in physical photon
as well as physical phonon. We also verify the dressing effect of vacuum,
where the non-vanishing mean value of photon number in ground state will
modify the ground state energy [30]. Beyond the single mode approximation, we
briefly discuss the effective interaction between photon modes induced by the
fluctuation of phonon vacuum. The dynamical Casimir effect is modeled by the
scattering processes converting phonon to photons, and the corresponding
transition amplitudes are obtained. We also derive the correlation function of
the backaction force, which is non-Gaussian distributed in our model. And we
also point out that the correlation function of backaction force is very
small, meaning the Langivan-type equation serves as a good approximmation
describing the mirror’s motion. Throughout this paper the systems is
considered in weak coupling regime and the dispassion is disregarded.

Figure 1: The sketch of our system: An optical cavity with one fixed mirror
and another attached to a spring. The cavity field acts radiation pressure on
the mirror and the motion of the mirror squeezes the cavity field.

This paper is organized as follows. In Sec.2, we give a brief introduction of
the model Hamiltonian and present the Feynman rules. In Sec.3, we study the
dressing effects of phonon, photons and the coupling strength. We also
demonstrate the existence of the virtual particles in vacuum state, which
shifts the ground state energy of the system. The vacuum induced photon-photon
interaction is also briefly discussed. In Sec.4, we compute the transition
amplitudes of the processes modeling dynamical Casimir effect and derive the
correlation function of the backaction force, whose non-Gaussian possibility
distribution is revealed. In Sec.5, we give a conclusion of our results and
give some discussion.

Throughout this paper the units are chosen as $\hbar=c=1$.

##  2 Model Hamiltonian and Feynamn rules

The schematic of the system is demonstrated in Fig.(1). The mass of the
movable mirror is $m$ and the strength of the spring is $k$. $L$ is the length
of the cavity when the spring is at relaxation. Under the action of the
spring, the mirror behaves like a simple harmonic oscillator and the radiation
pressure on the mirror can be regarded as driving force. In the meantime, the
motion of the mirror will also squeeze the cavity field non-adiabatically. In
our discussion, the radiation pressure is small compared to the elastic force
from the spring, thus the mirror’s sphere of activity is much smaller compared
to $L$.

The Hamiltonian of this system is given as $H=H_{0}+V$, with the free term

$H_{0}=\omega_{c}a^{\dagger}a+\omega_{m}b^{\dagger}b,$ (1)

where $a$ is the cavity field operator whose frequency is
$\omega_{c}=n{\pi}/L$, with $n$ for the arbitrarily chosen mode number, and
$b$ is the annihilation operator of the harmonic oscillator of form
$b=\sqrt{\frac{m\omega_{m}}{2}}(x+\frac{ip}{m\omega_{m}})$, with $x$ and $p$
being the position and momentum operators of the harmonic oscillator, thus the
frequency of harmonic oscillator is given by $\omega_{m}=\sqrt{\frac{k}{m}}$.
Note that here we use the single mode simplification which will be sufficient
for making our point. Meanwhile, there are two polarizations of the
electromagnetic field inside the cavity, which gives only a redundant freedom
in the effects we will discuss later [31]. So we will consider only one
polarization in this case. Throughout this paper, the exciton number of the
harmonic oscillator is referred as phonon number [3].

The radiation pressure on the mirror is given by [25, 26]

$F=\frac{\omega_{c}}{2L}(a+a^{\dagger})^{2},$ (2)

from which we obtain the interaction term as

$\begin{split}V&=Fx\\\
&=ga^{\dagger}a(b+b^{\dagger})+\frac{g}{2}[a^{2}+(a^{\dagger})^{2}](b+b^{\dagger}),\end{split}$
(3)

where the normal ordering technique is applied [30] and

$g=-\sqrt{\frac{1}{2m\omega_{m}}}\frac{\omega_{c}}{L}.$ (4)

The interaction term is split into two parts: $V=V_{om}+V_{DCE}$, the first
part

$V_{om}=ga^{\dagger}a(b+b^{\dagger})$ (5)

is called the standard optomechanical coupling term which conserves the number
of photon [32] and the second term

$V_{DCE}=\frac{g}{2}[a^{2}+(a^{\dagger})^{2}](b+b^{\dagger})$ (6)

is called dynamical Casimir effect term being responsible for the dynamical
Casmir effect. In this paper, the cavity length $L$ is a macroscopic quantity
and mode number $n$ is taken not very large. Thus we obtain the weak coupling
condition $\omega_{c},\omega_{m}{\gg}g$. In experimentation, the
optomechanical coupling term $V_{om}$ plays the dominant role and the effects
from dynamical Casimir effect term is hard to observe when the coupling $g$ is
too small [32].

Figure 2: Feynman rules in momentum space. (a) External lines. (b)
Propagators. (c) Optomechanical interaction vertices, the numbers denote the
number of the lines attached to the vertex. (d) Dynamical Casimir effect
interaction vertices, the numbers denote the number of the lines attached to
the vertex. Note that the cross mark in (a) and (b) can be any type of
vertices. And the boson symmetry factor is absorbed in the interaction
strength in (c) and (d). In (a) and (b), the legs are drawn short to save
space thus we name them ”Corgidiagrams”.

Like the case of Rabi model [33], this optomechanical system can also be
viewed as an anology to quantum electrodynamics. Thus we are encouraged to
apply the Feynman diagram technique. Since we are working in non-relativistic
regime, the field propagators goes only from the backward towards the future,
making us unable to use Wick’s theorem. After some algebra, the Feynman rules
is obtained and is shown in Fig.(2), whose detailed derivation is given in
Appendix A. Note that the Bose symmetry factor is absorbed into the coupling
strength in Fig.(2). From Fig.(2) we can see the two kinds of interaction in
Eq.(3) gives two types of the vertices, the optomechanical ones are marked by
squares while the dynamical Casimir effect ones marked by circles. Note that
the different kinds of interaction vertices give rise to different value of
coupling strength, as shown in Fig.(2). The time arrow in Fig.(2) and
throughout this paper is taken from left to right. Here are three more
constraints in building Feynman diagrams with Fig.(2): One, the energy
conserves at each vertices. Two, take the integral of all the undermined
energy appear in loop diagrams as ${\int}\frac{dE}{2\pi}$. And three, the
energy of the incoming and outgoing lines is conserved. Whose derivation can
also be found in Appendix A.

##  3 Dressing effects

In classical field theory, say, electrodynamics, the dynamics of the field is
given by Maxwell equation when the motion of source is specified. However, the
classical electrodynamics fails give the dynamics of the source under the
action of the field [34, 35]. This difficulty is inherited by the
semiclassical theory where the field are introduced as quantum variables while
the source are described by classical variables [1]. This implies a full
quantum theory is needed to give a self-consistent description of the
backaction effects. In quantum electrodynamics, the dynamics of electrons are
also quantized. As a result of which, bare electrons travelling through the
electromagnetic vacuum will get dressed with a cloud of virtual photons, which
gives self-consistence modifications to the physical quantities of the
electrons such as mass, charge and spin [1, 33, 36]. In this section, we will
discuss the dressing effects of the photons, phonon and the coupling strength
in our system, among which the dressing effect of phonon is the first
backaction effect we encounter.

###  3.1 Frequency shift of harmonic oscillator

In quantum electrodynamics, the loop diagrams of the fermion lines gives the
electromagnetic mass. In this part, we will mimic such procedure. To the
lowest order of $g$, there exists two one particle irreducible diagrams shown
in Fig.(3). Under the chain approximation [36], the frequency shift
$\delta\omega_{m}$ is given as

Figure 3: One particle irreducible diagrams of phonon propagator. (a) Diagram
of two virtual photons. (b) Diagram of four virtual photons.

$\delta\omega_{m}=\Sigma_{2}^{(2)}(\omega_{m})+\Sigma_{2}^{(4)}(\omega_{m}),$
(7)

where $-i\Sigma_{2}^{(2)}(\omega_{m})$ and $-i\Sigma_{2}^{(4)}(\omega_{m})$
are the values of Fig.(3a) and Fig.(3b), respectively. The superscripts denote
the order of diagrams while the subscripts denotes the number of intermediate
particles.

With the Feynman rules, the value of Fig.(3a) is

$\begin{split}-i\Sigma_{2}^{(2)}(\omega_{m})&=-\frac{g^{2}}{2}\int_{-\infty}^{+\infty}\frac{dE_{1}dE_{2}}{2\pi}\frac{i}{E_{1}-\omega_{c}+i\epsilon}\frac{i}{E_{2}-\omega_{c}+i\epsilon}\delta(E_{1}+E_{2}-\omega_{m})\\\
&=i\frac{g^{2}}{2}\frac{1}{2\omega_{c}-\omega_{m}}.\end{split}$ (8)

The integral in Eq.(8) is convergent. And the value of Fig.(3b) is

$\begin{split}-i\Sigma_{4}^{(2)}(\omega_{m})&=-g^{2}\int_{-\infty}^{+\infty}\frac{dE_{1}dE_{2}dE_{3}dE_{4}}{(2\pi)^{3}}\frac{i}{E_{1}-\omega_{c}+i\epsilon}\frac{i}{E_{2}-\omega_{c}+i\epsilon}\frac{i}{E_{3}-\omega_{m}+i\epsilon}\\\
&\times\frac{i}{E_{4}-\omega_{b}+i\epsilon}\delta(E_{1}+E_{2}+E_{3}+E_{4}-\omega_{m})\\\
&=ig^{2}\frac{1}{2\omega_{c}+\omega_{m}}.\end{split}$ (9)

This integral is also convergent. When the coupling $g$ is sufficient small,
the frequency shift of harmonic oscillator Eq.(7) vanishes due to the
$V_{DCE}$ being negligible. In the rotation wave approximation applied in [25,
26], only diagram Fig.(3a) survives.

Taking Eqs.(8) and (9) into Eq.(7) we obtain the frequency shift of harmonic
oscillator

$\delta\omega_{m}=-g^{2}(\frac{1}{2}\frac{1}{2\omega_{c}-\omega_{m}}+\frac{1}{2\omega_{c}+\omega_{m}}),$
(10)

which is identical what one can get from perturbation method to the state
$\ket{0,1}$, where the first number in the ket denotes photon number while the
second number denotes phonon number. From Fig.(3) we can see that the cloud of
virtual particles, both photon and phonon shifts the phonon frequency. As a
result, the observed frequency $\Omega_{m}=\omega_{m}+\delta\omega_{m}$ is no
longer only determined by the strength of the spring and the mass of mirror,
but also by the length of the cavity and the cavity field frequency.

Besides the correction to the physical quantities, another consequence of the
dressing effect is the renormalization of the field strength. In [25, 26], the
mean value of which is used to model the dynamical Casimir effect.

The field equation of the phonon field obtained by Heisenberg motion equation
is

$\frac{db}{dt}=-i[\omega_{m}b+ga^{\dagger}a+\frac{g}{2}[a^{2}+(a^{\dagger})^{2}],$
(11)

the solution of which can be formally written as

$b(t)=b_{0}(t)+\int dt^{\prime}G_{r}(t-t^{\prime})f(t^{\prime}),$ (12)

where $b_{0}(t)$ is the solution of the homogenous motion equation Eq.(28) and
the $G_{r}(t-t^{\prime})$ is the retarded Green function of the phonon field,
the source term $f(t)$ is

$f=-i[ga^{\dagger}a+\frac{g}{2}[a^{2}+(a^{\dagger})^{2}].$

Eq.(12) gives us a picturized description of a physical phonon being composed
of a bare phonon with a cloud of virtual particles. To keep in consistent with
canonical quantization condition, the interacting field equation Eq.(12) must
satisfy the boson statistics $[b,b^{\dagger}]=1$ which was imposed on the free
field $b_{0}$ originally. As a result of which, the free field $b_{0}$ takes a
percentage of $Z_{b}$ of the physical field $b$ and is no longer normalized.
That is why the coefficient $Z_{b}$ called the field strength normalization
factor. The field strength renormalization factor $Z_{b}$ can be computed via

$\begin{split}Z_{b}^{-1}&=1-\frac{d\Sigma(\Omega_{m})}{dE}\\\
&=1+g^{2}[\frac{1}{2(2\omega_{c}-\omega_{m})^{2}}-\frac{1}{(2\omega_{c}+\omega_{m})^{2}}].\end{split}$
(13)

From Eqs.(10) and (13) we can see $\Omega_{m}$ and $Z_{b}$ go divergent at
resonance $\omega_{m}=2\omega_{c}$. However, this does not cause any
inconsistency. The reason is, when tuned in resonance, the dynamical Casimir
effect term $V_{DCE}$ will introduce process $\text{b}{\rightarrow}2\text{a}$,
which is forbidden by energy conservation if there is detuning. This means the
phonon is no longer stable in resonance. As a result of which, the energy of
the phonon will no longer be determined but will experience a distribution
width due to the time-energy uncertainty [37]. The detailed computation of the
decay width and the life time of the phonon in resonance will be given in the
next section.

###  3.2 Frequency shift of cavity field

To the lowest order of $g$, there are two one particle irreducible diagrams of
photon lines demonstrated in Fig.(4), which gives the frequency shift as

$\delta\omega_{c}=\Pi^{(2)}_{om}(\omega_{c})+\Pi^{(2)}_{DCE}(\omega_{c}),$
(14)

where $-i\Pi^{(2)}_{om}(\omega_{c})$ and $-i\Pi^{(2)}_{DCE}(\omega_{c})$ are
the values of Fig.(4a) and Fig.(4b), respectively. The superscript denotes the
order of the diagrams while the subscripts $om$ and $DCE$ denote the type of
vertices. Here we note that Fig.(4a) survives when $V_{DCE}$ becomes
negligible and the photon frequency shift Eq.(14) vanishes in the rotation
wave approximation [25, 26].

The values of the diagrams in Fig.(4(a)) and Fig.(4(b)) are

$\begin{split}-i\Pi^{(2)}_{om}(\omega_{c})&=-g^{2}\int_{-\infty}^{\infty}\frac{dE_{1}dE_{2}}{2\pi}\frac{i}{E_{1}-\omega_{c}+i\epsilon}\frac{i}{E_{2}-\omega_{m}+i\epsilon}\delta(E_{1}+E_{2}-\omega_{c})\\\
&=ig^{2}\frac{1}{\omega_{m}},\\\
-i\Pi^{(2)}_{DCE}(\omega_{c})&=-\frac{3g^{2}}{2}\int_{-\infty}^{+\infty}\frac{dE_{1}dE_{2}dE_{3}dE_{4}}{(2\pi)^{3}}\frac{i}{E_{1}-\omega_{c}+i\epsilon}\frac{i}{E_{2}-\omega_{c}+i\epsilon}\frac{i}{E_{3}-\omega_{c}+i\epsilon}\frac{i}{E_{4}-\omega_{m}+i\epsilon}\\\
&\times\delta(E_{1}+E_{2}+E_{3}+E_{4}-\omega_{c})\\\
&=i\frac{3g^{2}}{2}\frac{1}{2\omega_{c}+\omega_{m}},\end{split}$ (15)

respectively. The integrals in the equations are convergent. Thus we obtain
the frequency shift of photons

$\delta\omega_{c}=-g^{2}(\frac{1}{\omega_{m}}+\frac{3}{2}\frac{1}{2\omega_{c}+\omega_{m}}).$
(16)

Figure 4: One particle irreducible diagrams of photon propagator. (a) Diagram
with two virtual particles. (b) Diagram with four virtual particles.

This result is identical to that from perturbation method to state
$\ket{1,0}$. Like the picture we introduced in the dress effect of the phonon,
a physical photon also contains a bare photon dressed in a cloud of virtual
particles, whose field strength normalization factor is given by

$\begin{split}Z_{a}^{-1}&=1-\frac{d\Pi(\Omega_{c})}{dE}\\\
&=1-\frac{3g^{2}}{(2\omega_{c}+\omega_{m})^{2}}.\end{split}$ (17)

###  3.3 Vertex functions

The irreducible vertex loop diagrams to the lowest order of $g$ is given in
Figs.(5a) and (5b), from which we can see that the optomechanical interaction
plays central part in this effect.

Figure 5: Irreducible vertex loop diagrams. (a) Loop diagram of dynamical
Casimir effect interaction vertex. (b) Loop diagram of optomechanical
interaction vertex. (c) Loop diagram of dynamical Casimir effect vertices with
more than two virtual photon lines. Diagram (c) gives correction to
Figs.(3(b)) and (4(b)).

The corresponding vertex functions are

$\begin{split}\Gamma_{DCE}(E_{1})&=\frac{ig^{3}}{4}\int\frac{dE}{2\pi}\frac{i}{E-\omega_{m}+i\varepsilon}\frac{i}{(E_{1}-E)-\omega_{c}+i\varepsilon}\frac{i}{(\omega_{m}+E-E_{1})-\omega_{c}+i\varepsilon}\\\
&=-\frac{ig^{3}}{4}\frac{1}{(E_{1}-\omega_{c}-\omega_{m})(\omega_{m}-2E_{1})},\\\
\Gamma_{opto}(E_{1})&=-\frac{ig^{3}}{8}\int\frac{dE}{2\pi}\frac{i}{E-\omega_{m}+i\varepsilon}\frac{i}{(E_{1}-E)-\omega_{c}+i\varepsilon}\frac{i}{(\omega_{m}+E-E_{1})-\omega_{c}+i\varepsilon}\\\
&=\frac{ig^{3}}{8}[\frac{1}{(E_{1}-\omega_{c}-\omega_{m})(2E_{1}-\omega_{m})}+\frac{1}{(\omega_{c}+E_{1})(\omega_{m}-2E_{1})}],\end{split}$
(18)

respectively, where the variable in the $\Gamma$ function denotes the energy
of the outgoing lines in the diagrams with undetermined energy. Thus the
correction to the coupling strength are given by $\Gamma_{DCE}(E_{1})$ and
$\Gamma_{opto}(E_{1})$, respectively.

In diagrams with many virtual particles like Fig.(3b) and Fig.(4b) we will
encounter the type of vertex graph in Fig.(5c)

$\begin{split}\Gamma_{x}(E_{1},E_{2})&=\frac{ig^{3}}{4}\int\frac{dE}{2\pi}\frac{i}{E-\omega_{m}+i\varepsilon}\frac{i}{E_{1}+E-\omega_{c}+i\varepsilon}\frac{i}{E_{2}-E-\omega_{c}+i\varepsilon}\\\
&=\frac{ig^{3}}{4}\frac{1}{(E_{1}-\omega_{c}-\omega_{m})(E_{1}+E_{2}-2\omega_{c})},\end{split}$
(19)

where there are two outgoing lines with undetermined energy marked by $E_{1}$
and $E_{2}$.

###  3.4 Vacuum effect

In ground state, the quantum harmonic oscillator is not at still but exhibits
zero-point motion, which squeezes the vacuum and may also generate radiation.
Due to the energy conservation, the radiated particles cannot be real. In this
part we will make use of the vacuum diagrams among which the one of lowest
order of $g$ is shown in Fig.(6). Note that this effect is purely induced by
dynamical Casimir effect interaction.

The value of Fig.(6) is

$\begin{split}-i\Lambda&=-\frac{g^{2}}{2}\int_{-\infty}^{+\infty}\frac{dE_{1}dE_{2}dE_{3}}{(2\pi)^{2}}\frac{i}{E_{1}-\omega_{c}+i\epsilon}\frac{i}{E_{2}-\omega_{c}+i\epsilon}\frac{i}{E_{3}-\omega_{m}+i\epsilon}\delta(E_{1}+E_{2}+E_{3})\\\
&=i\frac{g^{2}}{2}\frac{1}{2\omega_{c}+\omega_{m}},\end{split}$ (20)

which gives rise to the shift of the system’s ground state energy by

${\delta}E_{g}=\Lambda.$ (21)

There is a sketchy way of understanding this, we pretend there is a propagator
for the vacuum whose one particle irreducible loop diagram is Fig.(6), then
the shift of the vacuum energy can be obtained by chain approximation.

Eq.(21) is an exact match to that from perturbation theory to $\ket{0,0}$ and
verifies the result in [30]. And the integral in Eq.(21) is convergent. This
can be interpreted as the zero-point movement of the mirror will generate
”virtual radiation”.

Figure 6: The vacuum-vacuum diagram.

When more than one mode of the cavity field is taken into consideration, the
photons of the different modes can interact through the vacuum fluctuation of
the phonon field. In this case, the interaction term becomes

$V=\sum_{ij}g_{ij}(b+b^{\dagger})a_{i}^{\dagger}a_{j}+\frac{g_{ij}}{2}(b+b^{\dagger})(a_{i}a_{j}+a_{i}^{\dagger}a_{j}^{\dagger}),$
(22)

where $a_{i}$ and $a_{i}^{\dagger}$ denote the annihilation and creation
operator of the cavity field mode $i$ and
$g_{ij}=(\frac{1}{2})^{\frac{3}{2}}(-1)^{i+j}\frac{1}{L\sqrt{m}}\sqrt{\frac{\omega_{i}\omega_{j}}{\omega_{m}}}$
denotes the coupling strength when the photon lines of mode $i$ and $j$ attach
to the vertex. In deriving Eq.(22), the normal ordering technique is applied
to get rid of the infinite vacuum energy [28, 30].

To the lowest order of $g$, the effective photon-photon interaction diagram is
shown in Fig.(7). We can see this effect is induced solely by optomechanical
interaction. The value of Fig.(7) is thus given by

Figure 7: The photon-photon interaction diagram, where the photon propagator
can be of any mode.

$\begin{split}\sigma&=\sum_{ab}(-ig_{ab})^{4}\int\frac{dE}{2\pi}\frac{i}{E-\omega_{i}-\omega_{m}+i\varepsilon}\frac{i}{E-\omega_{k}-\omega_{m}+i\varepsilon}\frac{i}{E-\omega_{z}+i\varepsilon}\frac{i}{E-\omega_{k}-\omega_{c}-\omega_{b}+i\varepsilon}\\\
&=\sum_{ab}ig_{ab}^{4}[\frac{1}{(\omega_{i}-\omega_{k})(\omega_{i}+\omega_{m}-\omega_{a})(\omega_{i}+\omega_{m}-\omega_{k}-\omega_{l}-\omega_{b})}\\\
&+\frac{1}{(\omega_{k}-\omega_{i})(\omega_{k}+\omega_{m}-\omega_{a})(\omega_{m}-\omega_{l}-\omega_{b})}+\frac{1}{(\omega_{a}-\omega_{i}-\omega_{m})(\omega_{a}-\omega_{k}-\omega_{m})(\omega_{a}-\omega_{k}-\omega_{c}-\omega_{b})}\\\
&+\frac{1}{(\omega_{j}+\omega_{b}-\omega_{m})(\omega_{i}+\omega_{j}+\omega_{b}-\omega_{k}-\omega_{m})(\omega_{i}+\omega_{j}+\omega_{b}-\omega_{a})}].\end{split}$
(23)

Note here we ignored the Boson statistics of the outgoing photons and this is
a process of order $g^{4}$.

##  4 Dynamical Casimir effect and backaction

The dynamical Casimir effect is energy transferring from the mirror to the
field in nature. In scattering theory language, this can be modeled by the
processes $r\text{b}{\rightarrow}s\text{a}$ when parameters are properly tuned
as $r\omega_{m}=s\omega_{a}$, where $r$ and $s$ are the number of the phonon
and photon involved. Throughout this section we will be working in tree-level
diagrams.

In Fig.(8) we presented some processes of dynamical Casimir effect to the
lowest order of $g$. From which we can see that the there are always two
phonon taking part in the dynamical Casimir effect, i.e. $r=2$. First we
consider the process $\text{b}\rightarrow 2\text{a}$ shown in Fig.(8a). The
decay width of the phonon can be directly seen from Fig.(8a) as

$\Gamma=g^{2}.$

Thus the life time of phonon is obtained

$\tau=\frac{1}{\Gamma}=\frac{1}{g^{2}}.$ (24)

The diagrams of $k\text{b}{\rightarrow}2\text{a}$ with $k=2,3,4$ are shown in
Figs.(8b), (8c) and (8d), respectively. From which we can deduce that for
processes $k\text{b}{\rightarrow}2\text{a}$ with arbitrary $k$, there will be
$k-1$ optomechanical vertices and one dynamical Casimir effect vertices in the
diagrams, thus the corresponding transition amplitude

$A_{k{\rightarrow}2}\sim\frac{1}{2^{k-1}}g^{k}.$ (25)

Figure 8: Diagrams of the lowest order dynamical Casimir effect process: (a)
Diagram of decay of a phonon, (b) Diagram of $2\text{b}{\rightarrow}2\text{a}$
process, (c) Diagram of $3\text{b}{\rightarrow}2\text{a}$ process, (d) Diagram
of $4\text{b}{\rightarrow}2\text{a}$ process. Note that we dropped the
markings for energies because they can be easily be deduced by the energy
conservation rules at each vertices and of the whole process.

Since the coupling constant $g$ is small, the transition amplitude Eq.(25)
decrease rapidly with growing $k$. This is qualitatively in agreement with the
result obtained in the previous work [32].

When photons are generated out of the vacuum, they take the energy away from
the mirror. In the meantime, they exert radiation pressure on the mirror. This
is the most commonly known backaction effect. The statistical nature of the
backaction force has been studied persistently [17, 18, 19, 20, 21, 23, 38]
and the debate on whether the mirror will exhibit Brownian motion has also
lasted for a very long time [20, 23, 38]. In the following part, we will use
the Feynman diagram approach to explore the statistic nature of the backaction
force of our system.

The basic tool will be used is the correlation function of two operators $A$
and $B$

$\bra{\Omega}T[A(t_{1})B(t_{2})]\ket{\Omega}=\frac{\bra{0,0}T[U_{I}(+\infty,t_{1})A_{I}(t_{1})U_{I}(t_{1},t_{2})B_{I}(t_{2})U_{I}(t_{2},-\infty)]\ket{0,0}}{\bra{0,0}U_{I}(+\infty,-\infty)\ket{0,0}},$
(26)

where $T$ is the time-ordering symbol, $U_{I}(t)$ is the evolution operator in
interaction picture and $A$ and $B$ are arbitrary time-dependent operators,
$\ket{\Omega}$ is the ground state of the full Hamiltonian.

Figure 9: Diagrams of the force-force correlation function to the lowest
order: (a) connected diagram, (b) unconnected diagram. The empty circles at
the end of the photon lines denotes they are not of interaction origin.

Taking $A$ and $B$ as the backaction force operator $F$ in Eq.(2) we obtain
the its correlation function. Power expand the evolution operator gives the
force-force correlation function

$\bra{\Omega}T[F(t_{1})F(t_{2})]\ket{\Omega}=\frac{\omega_{c}^{2}}{4L^{2}}e^{-2i\omega_{c}(t_{1}-t_{2})}$
(27)

to the lowest order of $g$, the corresponding diagrams are given in
Fig.(9),where $G_{0}(t)$ is the free Green function of Eq.(1). Note that the
unconnected diagram Fig.(9b) is ruled out by the denominator in Eq.(26).
Eq.(27) demonstrates that the backaction force at different time have non-
vanishing correlations and is of order $L^{-3}$, which is a small number. This
means the Langevin-type equations serves as a approximation of modeling the
motion of the mirror to the order of $L^{-3}$. The techniques used in deriving
Eq.(27) can be found in Appendix A.

##  5 Conclusion and Discussion

In this paper we have developed the Feynman diagram technique to study the
fully quantized optomechanical system. We studied the dressing effects of
phonon, photon and strength coupling, and demonstrated how the vacuum
fluctuation modifies the physical quantities and physical particles of the
system. The phenemenon of non-zero photon population will modify the ground
state energy is verified and modeled as virtual radiation process. The
dynamical Casimir effect is modeled elegantly by the scattering process
converting the phonon pair into photons and the corresponding transition rates
are derived. We also revealed the non-Gaussian nature of the possibility
distribution of the backaction force and showed the Langivan-type equation
serves as a good approximation of modeling the mirror’s motion.

The optimechanical system studied in this paper is conceptual a model of
interaction between the mechanical degree of freedom and quantized field, thus
implying many similarity to the black hole with Hawking radiation.
Unfortunately, due to the lack of a manifestable quantum gravity theory, we
can now modeling the black hole evaporation only on a semiclassical level [8].
The techniques developed and new effects revealed in our system may be of help
in the further investigation of the quantum effects in black hole evaporation,
through which the road towards the theory of quantum gravity may be lighted
[39].

### Acknowledgements

We thank W. H. Sang for comments on early drafts of this paper. Y. S. Cao
thanks K. F. Lyu for discussion about the loop diagram computation. We thank
R. Zhang for useful discussion.

##  Appendix A Green function and Feynman rules

The motion equation of free field $a$ dictated by Eq.(1) is

$(\frac{d}{dt}+i\omega_{c})a=0,$ (28)

the corresponding Green function is

$G_{a}(t-t^{\prime})=\int\frac{dE}{2\pi}\frac{i}{E-\omega_{c}+i\varepsilon}e^{-iE(t-t^{\prime})},$
(29)

which is in consistence with the one obtained from Huygens principle

$G_{a}(t-t^{\prime})=\theta(t-t^{\prime})\bra{0}a(t)a^{\dagger}(t^{\prime})\ket{0}$

. In energy representation

$G_{a}[E]=\frac{i}{E-\omega_{c}+i\varepsilon},$ (30)

where we used parenthesis and square brackets to distinguish the primitive and
Fourier transformed functions.

The free Green operator of the system satisfies

$(\frac{d}{dt}-iH_{0})G_{0}(t-t^{\prime})=\delta(t-t^{\prime}),$ (31)

whose solution is obtained as [37]

$G_{0}(t-t^{\prime})=\theta(t-t^{\prime})e^{-iH_{0}(t-t^{\prime})},$ (32)

where $\theta(x)$ is the step function. This gives us the relation between
Green operator and time free evolution operator $U_{0}(t)$ as

$G(t-t^{\prime})=\theta(t-t^{\prime})U_{0}(t-t^{\prime}).$ (33)

The free Green function Eq.(29) can thus be written as

$G_{a}(t-t^{\prime})=\theta(t-t^{\prime})\bra{1}G_{0}(t-t^{\prime})\ket{1}.$
(34)

For multi-particle propagation we have

$G_{a}^{n}(t-t^{\prime})=\theta(t-t^{\prime})\bra{n}G_{0}(t-t^{\prime})\ket{n}.$
(35)

The $S$-matrix

$\begin{split}S&=U_{I}(+\infty,-\infty)\\\
&=\sum_{n=0}^{\infty}\frac{(-i)^{n}}{n!}T\int dt_{1}dt_{2}\cdots
dt_{n-1}dt_{n}[V(t_{1})V(t_{2})\cdots V(t_{n-1})V(t_{n})],\end{split}$ (36)

where $T$ is time-ordering symbol and $U_{I}(t,t^{\prime})$ is the time
evolution operator in interaction picture.

Express the arbitrary $S$-matrix element in Eq.(36) in Schrodinger picture,
its general term formula reads

$\begin{split}&\frac{(-i)^{n}}{n!}\bra{\alpha}T\int dt_{1}dt_{2}\cdots
dt_{n-1}dt_{n}[V(t_{1})V(t_{2})\cdots V(t_{n-1})V(t_{n})]\ket{\beta}\\\
&=\frac{(-i)^{n}}{n!}\bra{\alpha}\int dt_{1}dt_{2}\cdots
dt_{n-1}dt_{n}[G_{0}(+\infty,t_{1})VG_{0}(t_{1},t_{2})\cdots
G_{0}(t_{n-1},t_{n})VG_{0}(t_{n},-\infty)]\ket{\beta},\end{split}$ (37)

where two facts used in deriving the formula above. First, state
$\ket{\alpha}$ obeys the following rule evolving in time

$\ket{\alpha,t}=U_{0}^{\dagger}(t,t^{\prime})\ket{\alpha,t^{\prime}},$

and so is for $\ket{\beta}$. Second, the explicit form the Green operator in
Eq.(33) allows us to eliminate the time-ordering symbol in the second row due
to the step function in Green operator.

For a general term in Eq.(37), the Fourier transformation

$\begin{split}&\int dt_{1}dt_{2}\cdots
dt_{n-1}dt_{n}[G_{0}(+\infty,t_{1})VG_{0}(t_{1},t_{2})\cdots
G_{0}(t_{n-1},t_{n})VG_{0}(t_{n},-\infty)]\\\ &=\frac{1}{2\pi}\int
dE_{1}dE_{2}\cdots dE_{n-1}dE_{n}G_{0}[E_{1}]VG_{0}[E_{2}]\cdots
G_{0}[E_{n-1}]VG_{0}[E_{n}]\delta(E_{1}-E_{2})\delta(E_{2}-E_{3})\cdots\delta(E_{n-1}-E_{n})\end{split}$
(38)

gives its formula in energy representation. The reason why the interaction $V$
does not contain energy variable $E$ can be understood by the whole system
being closed, which means the energy is conserved in each interaction
vertices.

From Eq.(35), we have

$\bra{n}G_{0}[E]\ket{n}=\int\frac{dE_{1}}{2\pi}\frac{dE_{2}}{2\pi}\cdots\frac{dE_{n-1}}{2\pi}\frac{dE_{n}}{2\pi}G_{a}[E_{1}]G_{a}[E_{2}]\cdots
G_{a}[E_{n-1}]G_{a}[E_{n}]\delta(\sum_{i=1}^{n}E_{i}-E),$ (39)

whose diagrammatic representation is the Feynman rules introduced in Sec. 2.

## References

  * [1] C. Itykson and J-B. Zuber, Quantum Field Theory, Dover Publications (2006).
  * [2] S. Weinberg, Rev. Mod. Phys. 61, 1 (1989).
  * [3] P. D. Nation, J. R. Johansson J.R., M. P. Blencowe and F. Nori, Rev. Mod. Phys. 1, 84 (2012).
  * [4] P. W. Milonni, The Quantum Vacuum: An Introduction to Quantum Electrodynamics, Academic Press (New York) (1994).
  * [5] G. T. Moore, J. Math. Phys. 11, 2679 (1970).
  * [6] S. A. Fulling and P. C. W. Davies, Proc. R. Soc. Lond. A: Math. Phys. Sci. 348, 393 (1976).
  * [7] S. W. Hawking, Commun. Math. Phys. 43, 199 (1975).
  * [8] A. Fabbri and J. Navarro-Salas, Modeling Black Hole Evaporation, Imperial College Press (2005).
  * [9] L. Parker, Phys. Rev. 183, 1057 (1969).
  * [10] W. Kim, J. H. Brownell and R. Onofrio, Phys. Rev. Lett. 96, 200402 (2006).
  * [11] W. Qin, V. Macrì, A. Miranowicz, S. Savasta and F. Nori, Phys. Rev. A 100, 062501 (2019).
  * [12] W. G. Unruh, Phys. Rev. D 14, 870 (1976).
  * [13] R. D. Carlitz and R. S. Willey, Phys. Rev. D 36, 2327 (1987).
  * [14] R. D. Carlitz and R. S. Willey, Phys. Rev. D 36, 2336 (1987).
  * [15] J. Haro and E. Elizalde, Phys. Rev. D 77, 045011 (2008).
  * [16] P. A. M. Neto and L. A. S. Machado, Phys. Rev. A 54, 3420 (1996).
  * [17] M. F. Maghrebi, R. Golestanian and M. Kardar, Phys. Rev. D 87, 025016 (2013).
  * [18] L. H. Ford and A. Vilenkin, Phys. Rev. D 25, 2569 (1982).
  * [19] P. A. M. Neto, J. Phys. A: Math. Gen. 27, 2167 (1994).
  * [20] Q. D. Wang and W. G. Unruh, Phys. Rev. D 89, 085009 (2014).
  * [21] Q. D. Wang and W. G. Unruh, Phys. Rev. D 92, 063520 (2015).
  * [22] T. J. Kippenberg and K. J. Vahala, Science 321, 1172 (2008).
  * [23] G. Gour and L. Sriramkumar, Found. Phy. 29, 1917 (1999).
  * [24] C. K. Law, Phys. Rev. A 3, 2537 (1995).
  * [25] S. Butera and I. Carusotto, EPL 128, 24002 (2019).
  * [26] S. Butera and I. Carusotto, Phys. Rev. A 99, 053815 (2019).
  * [27] R. Golestanian and M. Kardar, Phys. Rev. A 58, 1713 (1998).
  * [28] F. Armata and R. Passante, Phys. Rev. D 91, 025012 (2015).
  * [29] F. Armata, M. S. Kim, S. Butera, L. Rizzuto, and R. Passante, Phys. Rev. D 96, 045007 (2017).
  * [30] S. Butera and R. Passante, Phys. Rev. Lett. 111, 060403 (2013).
  * [31] N. F. Del Grosso, F. C. Lombardo and P. I. Villar, arXiv:1910.10218.
  * [32] V. Marci, A. Ridolfo, O. D. Stefano, A. F. Kockum, F. Nori and S. Savasta, Phys. Rev. X 8, 011031 (2018).
  * [33] O. D. Stefano, R. Stassi, L. Garziano, A. F. Kockum, S. Savasta and F. Nori, New J. Phys. 19 053010 (2017).
  * [34] J. P. Liu, Electrodynamics, Higher Education Press (2004).
  * [35] Z. S. Zhang, Electrodynamics and Special Relativity, Peking University Press (2013).
  * [36] M. E. Peskin and D. V. Schroeder, An Introduction to Quantum Field Theory, Perseus Books (1995).
  * [37] J. P. Liu, Quantum Mechanics, Higher Education Press (2012).
  * [38] M. -T. Jaekel and S. Reynaud, J. Phys. I (France) 3, 1 (1993).
  * [39] A. Gomberoff and D. Malrof, Lectures on Quantum Gravity, Springer (2005).