# Spontaneous decay of artificial atoms in a three-qubit system

Ya. S. Greenberg<EMAIL_ADDRESS>Novosibirsk State Technical
University, Novosibirsk, Russia  A. A. Shtygashev  Novosibirsk State Technical
University, Novosibirsk, Russia  A. G. Moiseev  Novosibirsk State Technical
University, Novosibirsk, Russia

###### Abstract

We study the evolution of qubits amplitudes in a one-dimensional chain
consisting of three equidistantly spaced noninteracting qubits embedded in an
open waveguide. The study is performed in the frame of single-excitation
subspace, where the only qubit in the chain is initially excited. We show that
the dynamics of qubits amplitudes crucially depend on the value of $kd$, where
$k$ is the wave vector, $d$ is a distance between neighbor qubits. If $kd$ is
equal to an integer multiple of $\pi$, then the qubits are excited to a
stationary level. In this case, it is the dark states which prevent qubits
from decaying to zero even though they do not contribute to the output
spectrum of photon emission. For other values of $kd$ the excitations of
qubits exhibit the damping oscillations which represent the vacuum Rabi
oscillations in a three-qubit system. In this case, the output spectrum of
photon radiation is determined by a subradiant state which has the lowest
decay rate. We also investigated the case with the frequency of a central
qubit being different from that of the edge qubits. In this case, the qibits
decay rates can be controlled by the frequency detuning between the central
and the edge qubits.

qubits, microwave circuits, waveguide, transmission line, quantum measurements

###### pacs:

84.40.Az, 84.40.Dc, 85.25.Hv, 42.50.Dv, 42.50.Pq

##  I Introduction

Superconducting qubits coupled to photons propagating in an open waveguide
Asta2010 ; Hoi2011 ; Loo2013  allow for the investigation of the fascinating
world of quantum light-matter interactions in one dimensionRoy2017 ; Gu2017 ;
Shev2019 . Even though the properties of multi-qubit 1D systems have been
extensively studied both theoretically Alb2019 ; Zhang2019 ; Ruos2017 ;
Lalum2013 ; Chang2012  and experimentally Mirho2019 ; Brehm2021 ; Loo2014 ,
less attention has been paid to a detail investigation of dynamic properties
of the few-qubit systems which are the building blocks for quantum gates
Barn2017 ; Kenfa2018 . Moreover, the scaling laws for decay rates that have
been found for multi-qubit systems Tsoi2008  cannot obviously be applied for
systems containing few qubits.

As is known, the superconducting qubits can be technologically addressed and
controlled individually Shev2019 . Therefore, it is important to know the
evolution of the probability amplitude of any qubit in a superconducting
circuitry.

Here we explore the dynamic properties of a quantum circuit consisting of a
three-qubit linear chain, which is strongly coupled to a common waveguide. The
motivation for this choice is that it is a simplest system with non-trivial
properties for which a full analytical treatment can be obtained. We
investigate a dynamic behavior of the qubits amplitudes with the only qubit in
the chain being initially excited. For three-qubit system, we find the
analytic expressions for the complex energies and for the collective states
which define a temporal behavior of qubits amplitudes. We show that the
dynamics of qubits amplitudes crucially depend on the value of $kd$, where $k$
is the wave vector, $d$ is a distance between neighbor qubits. If $kd$ is
equal to integer multiple of $\pi$, the qubits are excited to a stationary
level. In this case, these are the dark states which prevent qubits from
decaying to zero even though they do not contribute to the output photon
spectrum. For other values of $kd$, the excitations of qubits have oscillatory
behavior and are gradually damped out to zero. We also investigate the case
with the frequency of a central qubit being different from that of the edge
qubits. In this case, the qibits decay rates can be controlled by the
frequency detuning between the central and the edge qubits. As the frequency
detuning between central and edge qubits increases, the decay rates of qubits
are also increases. This property is very important for the implementation of
the efficient control and readout protocols where a fast reset of the excited
qubits to their ground state is essential Zhou2021 ; Maq2018 .

The paper is structured as follows.

In Sec. II, we begin by introducing a Jaynes-Cummings Hamiltonian for atom-
light interactions. We truncate the Hilbert space to a single-excitation
subspace and obtain a set of linear integro-differential equations for the
qubits amplitudes.

In Sec. III, in the frame of Wigner-Weisskopf approximation we derive a set of
linear differential equations for the qubits amplitudes $\beta_{n}(t)$, which
allow for a direct numerical simulations.

A comprehensive analysis of the dynamics of the qubits amplitudes is given in
Sec. IV. We show that for $kd=n\pi$ the qubits amplitudes become ”frozen” at
the constant level. The reason for this is the dark states which prevent
qubits from decaying to zero. For the values of $kd$ which are not integer
multiple of $\pi$ the qubits amplitudes gradually damp out to zero. In this
section we also calculated the probability amplitude of the photon emission.
We find, that for $kd=(2n+1)\pi$ the evolution of the photon amplitude
consists of clearly seen steps. These steps can be attributed to the
interrelation between temporal behaviors of the different qubits amplitudes.

In Sec. V we calculate a spectral density of photon radiation from a three-
qubit chain. We show that for $kd=n\pi$ the dark states do not contribute to
output radiation near resonance. In this case, the output spectral density has
a Lorentzian lineshape. For $kd=(2n+1)\pi/2$ a spectral density exhibits two
peaks which are a signature of the vacuum Rabi oscillations of qubits
amplitudes. In general, if $kd$ is not equal to an integer multiple of $\pi$,
the width of a spectral line is defined by the deep subradiant states.

In Sec. VI we formulate the dynamics of a three-qubit chain with the aid of a
non-Hermitian Hamiltonian which is obtained after the elimination of the
photon variables. We find the collective states which are eigenvectors of non-
Hermitian Hamiltonian and show how the qubits amplitudes $\beta_{n}(t)$ can be
expressed in terms of these collective states.

In Sec. VII we consider a three-qubit system in which the frequency
$\Omega_{0}$ of the central qubit is different from that of the edge qubits.
We show that using the detuning $\delta\Omega=\Omega-\Omega_{0}$ as external
parameter we can control the decay rates of the qubits’ amplitudes. As the
frequency detuning between central and edge qubits increases, the decay rates
of qubits are also increases. For this case, a spectral density of photon
radiation is also calculated for different values of $\delta\Omega$.

The main results of the paper are summarized in the concluding Section VIII.

##  II Formulation of the problem

We consider a linear chain of three equally spaced qubits which are coupled to
photon field in an open waveguide (see Fig. 1).

Figure 1: Schematic illustration of a single-excitation subspace for a three-
qubit chain in an open waveguide. (a) A single qubit is excited, two qubits
are in the ground state. (b) three qubits are in the ground state and a single
photon propagates in the waveguide.

A distance between neighbor qubits is equal to $d$. The Hilbert space of every
qubit consists of the excited state $|e\rangle$ and the ground state
$|g\rangle$. The Hamiltonian which accounts for the interaction between qubits
and the electromagnetic field is as follows (we use $\hbar=1$ throughout the
paper):

$H=H_{0}+\sum\limits_{k}{\omega_{k}a_{k}^{+}a_{k}}+H_{int}$ (1)

where $H_{0}$ \- is Hamiltonian of bare qubits.

$H_{0}=\frac{1}{2}\sum\limits_{n=1}^{3}{\left({1+\sigma_{z}^{(n)}}\right)\Omega_{n}}$
(2)

$H_{\operatorname{int}}=\sum\limits_{n=1}^{3}{}\sum\limits_{k}{}g_{k}^{(n)}e^{-ikx_{n}}\sigma_{-}^{(n)}a_{k}^{+}+h.c.$
(3)

The quantity $g_{k}^{(n)}$ in (3) is the coupling between $n$-th qubit and the
photon field in a waveguide. Below we consider a single-excitation subspace
with either a single photon is in a waveguide and all qubits are in the ground
state, Fig. 1b, or there are no photons in a waveguide with the only n-th
qubit in the chain being excited, Fig. 1a. Therefore, we truncate Hilbert
space to the following states:

$\begin{array}[]{l}\left|{G,1_{k}}\right\rangle=\left|{g_{1},g_{2},g_{3}}\right\rangle\otimes\left|{1_{k}}\right\rangle;\\\
\left|{1,0_{k}}\right\rangle=\left|{e_{1},g_{2},g_{3}}\right\rangle\otimes\left|{0_{k}}\right\rangle;\\\
\left|{2,0_{k}}\right\rangle=\left|{g_{1},e_{2},g_{3}}\right\rangle\otimes\left|{0_{k}}\right\rangle;\\\
\left|{3,0_{k}}\right\rangle=\left|{g_{1},g_{2},e_{3}}\right\rangle\otimes\left|{0_{k}}\right\rangle\\\
\end{array}$ (4)

The Hamiltonian (3) preserves the number of excitations (number of excited
qubits + number of photons). In our case the number of excitations is equal to
one (see Fig. 1). Therefore, at any instant of time the system will remain
within a single-excitation subspace. The wave function of an arbitrary single-
excitation state can then be written in the form:

$\left|\Psi\right\rangle=\sum\limits_{n=1}^{3}{\beta_{n}(t)e^{-i\Omega_{n}t}}\left|{n,0_{k}}\right\rangle+\sum\limits_{k}{\gamma_{k}(t)e^{-i\omega_{k}t}}\left|{G,1_{k}}\right\rangle$
(5)

where $\beta_{n}(t)$ is the amplitude of $n$-th qubit, $\gamma_{k}(t)$ is a
single-photon amplitude which is related to a spectral density of spontaneous
emission.

$S(\omega_{k},t)=|\gamma_{k}(t)|^{2}.$ (6)

The function (5) is normalized to unity:

$\sum_{n=1}^{3}|\beta_{n}(t)|^{2}+\sum_{k}|\gamma_{k}(t)|^{2}=1.$ (7)

From (7) we can find the full probability of photon emission from the three-
qubit system.

$P_{ph}(t)=\sum_{k}|\gamma_{k}(t)|^{2}=1-\sum_{n=1}^{3}|\beta_{n}(t)|^{2}.$
(8)

In fact, the quantity $P_{ph}(t)$ is the probability to find the emitted
photon at the moment $t$.

The equations for qubits amplitudes $\beta_{n}(t)$ in (5) can be found from
time-dependent Schrodinger equation $id|\Psi\rangle/dt=H|\Psi\rangle$. For the
amplitudes $\gamma_{k}(t),\beta_{n}(t)$ we obtain:

$\begin{array}[]{l}\frac{{d\beta_{n}}}{{dt}}=-\sum\limits_{k}{|g_{k}^{(n)}|^{2}}\int\limits_{0}^{t}{\beta_{n}(t^{\prime})e^{-i(\omega_{k}-\Omega_{n})(t-t^{\prime})}dt^{\prime}}\\\
-\sum\limits_{m\neq
n}^{3}{}\sum\limits_{k}{g_{k}^{*(n)}g_{k}^{(m)}e^{-ik(x_{m}-x_{n})}e^{i(\Omega_{n}-\Omega_{m})t}}\\\
\times\int\limits_{0}^{t}{\beta_{m}(t^{\prime})e^{-i(\omega_{k}-\Omega_{m})(t-t^{\prime})}dt^{\prime}}\\\
\end{array}$ (9)

$\gamma_{k}(t)=-i\sum\limits_{n=1}^{3}{}g_{k}^{(n)}e^{-ikx_{n}}\int\limits_{0}^{t}{\beta_{n}(t^{\prime})e^{i(\omega_{k}-\Omega_{n})t^{\prime}}dt^{\prime}}$
(10)

According to (10) there are no photons in the system at $t=0$. Our goal is to
find the evolution of the amplitudes $\beta_{n}(t)$ for any qubit in the chain
when the only $n_{0}$-th qubit ($n_{0}=1,2$, or $3$) is initially excited:

$\begin{cases}&\beta_{n_{0}}(0)=1,\\\ &\beta_{n}(0)=0,\;\;n\neq
n_{0}.\end{cases}$ (11)

##  III Equations for the qubits amplitudes

We assume that the first and the third qubits are identical
($\Omega_{1}=\Omega_{3}\equiv\Omega$, $g_{k}^{(1)}=g_{k}^{(3)}\equiv g_{k}$).
The frequency and the coupling of the second qubit are different
($\Omega_{2}\equiv\Omega_{0}$, $g_{k}^{(2)}\equiv g_{k}^{(0)})$. A distance
between central qubit and the edge qubits is equal to $d$. We take the origin
in the location of the second qubit: $x_{1}=-d,\;x_{2}=0,\;x_{3}=+d$. In the
frame of Wigner-Weisskopf approximation, the equations (9) for the qubits
amplitudes can be reduced to the following set of linear differential
equations (see Appendix A for the derivation).

$\begin{array}[]{l}\frac{{d\bar{\beta}_{1}}}{{dt}}=-\frac{\Gamma}{2}\bar{\beta}_{1}(t)-i\frac{{\Omega-\Omega_{0}}}{2}\bar{\beta}_{1}(t)\\\
\\\
-\bar{\beta}_{2}(t)\frac{1}{2}\left({\frac{{\Omega_{0}}}{\Omega}}\right)^{1/2}\sqrt{\Gamma\Gamma_{0}}e^{ik_{0}d}-\frac{\Gamma}{2}\bar{\beta}_{3}(t)e^{2ikd}\\\
\end{array}$ (12)

$\begin{array}[]{l}\frac{{d\bar{\beta}_{2}}}{{dt}}=-\frac{{\Gamma_{0}}}{2}\bar{\beta}_{2}(t)+i\frac{{\Omega-\Omega_{0}}}{2}\bar{\beta}_{2}(t)\\\
\\\
-\frac{1}{2}\left({\frac{\Omega}{{\Omega_{0}}}}\right)^{1/2}\sqrt{\Gamma\Gamma_{0}}e^{ikd}\left({\bar{\beta}_{1}(t)+\bar{\beta}_{3}(t)}\right)\\\
\end{array}$ (13)

$\begin{array}[]{l}\frac{{d\bar{\beta}_{3}}}{{dt}}=-\frac{\Gamma}{2}\bar{\beta}_{3}(t)-i\frac{{\Omega-\Omega_{0}}}{2}\bar{\beta}_{3}(t)-\bar{\beta}_{1}(t)\frac{\Gamma}{2}e^{2ikd}\\\
\\\
-\bar{\beta}_{2}(t)\frac{1}{2}\left({\frac{{\Omega_{0}}}{\Omega}}\right)^{1/2}\sqrt{\Gamma\Gamma_{0}}e^{ik_{0}d}\\\
\end{array}$ (14)

where

$\begin{cases}\bar{\beta}_{1,3}(t)&=e^{-i(\Omega-\Omega_{0})t/2}{\beta}_{1,3}(t),\\\
\bar{\beta}_{2}(t)&=e^{i(\Omega-\Omega_{0})t/2}{\beta}_{2}(t).\end{cases}$
(15)

$k=\Omega/v_{g}$, $k_{0}=\Omega_{0}/v_{g}$, $\Gamma$, $\Gamma_{0}$ are the
rates of spontaneous emission into the waveguide mode from edge qubits and
from the central qubit, respectively.

##  IV The dynamics of three identical qubits

For three identical qubits we obtain from (12)-(14) the following equations:

$\frac{{d\beta_{n}}}{{dt}}=-\frac{\Gamma}{2}\sum\limits_{m=1}^{3}{}\beta_{m}(t)e^{ikd\left|{m-n}\right|},\quad(n=1,2,3)$
(16)

which can be expanded explicitly as:

$\begin{array}[]{l}\frac{{d\beta_{1}}}{{dt}}=-\frac{\Gamma}{2}\beta_{1}(t)-\frac{\Gamma}{2}\beta_{2}(t)e^{ikd}-\frac{\Gamma}{2}\beta_{3}(t)e^{2ikd}\\\
\\\
\frac{{d\beta_{2}}}{{dt}}=-\frac{\Gamma}{2}\beta_{2}(t)-\frac{\Gamma}{2}e^{ikd}\left({\beta_{1}(t)+\beta_{3}(t)}\right)\\\
\\\
\frac{{d\beta_{3}}}{{dt}}=-\frac{\Gamma}{2}\beta_{3}(t)-\frac{\Gamma}{2}\beta_{1}(t)e^{2ikd}-\frac{\Gamma}{2}\beta_{2}(t)e^{ikd}\\\
\end{array}$ (17)

where $k=\Omega/v_{g}$, and for simplicity we remove the bar over
$\beta_{n}(t)$.

The general solution of equations (17) can be written as follows:

${\rm{\beta}}_{i}(t)=\sum\limits_{m=1}^{3}{}a_{m}^{(i)}e^{\lambda_{m}t};\;i=1,2,3$
(18)

The quantities $\lambda_{m}$ are characteristic roots, which can be found by
equating to zero the determinant of equations (17)

$\left({\lambda_{i}+\frac{\Gamma}{2}}\right)\delta_{mn}+\frac{\Gamma}{2}e^{ikd\left|{m-n}\right|}\left({1-\delta_{mn}}\right)$
(19)

The quantities $a_{m}^{(n)}$ in (18) are defined by initial conditions:

$\begin{cases}&\sum_{m=1}^{3}a_{m}^{(n_{0})}=1,\\\
&\sum_{m=1}^{3}a_{m}^{(i)}=0,\;\;i\neq n_{0}.\end{cases}$ (20)

where $n_{0}$ is the number of initially excited qubit.

From determinant of (19) we obtain the exact expressions for $\lambda_{m}$:

$\begin{array}[]{*{20}c}{\lambda_{1}=-\frac{\Gamma}{4}e^{i2kd}-\frac{\Gamma}{4}e^{ikd}\sqrt{e^{i2kd}+8}-\frac{\Gamma}{2}}\\\
\\\
{\lambda_{2}=-\frac{\Gamma}{4}e^{i2kd}+\frac{\Gamma}{4}e^{ikd}\sqrt{e^{i2kd}+8}-\frac{\Gamma}{2}}\\\
\\\ {\lambda_{3}=\frac{\Gamma}{2}e^{i2kd}-\frac{\Gamma}{2}}\\\ \end{array}$
(21)

From (21) we see that the decay rates $\Gamma_{i}=-\mathrm{Re}(\lambda_{i})$
depend on $kd$. This dependence is shown in Fig. 2.

Figure 2: Dependence of $\mathrm{Re}(\lambda_{i})$ on $kd$.
$\mathrm{Re}(\lambda_{1})$-red line, $\mathrm{Re}(\lambda_{2})$-blue line,
$\mathrm{Re}(\lambda_{3})$-green line. All values of
$\mathrm{Re}(\lambda_{i})$ are negative.

In general, as is seen from (21), all $\Gamma_{i}$’s contribute to the decay
rate of a concrete qubit. More important is that a sum

$\sum_{i=1}^{3}\lambda_{i}=-\frac{3}{2}\Gamma$ (22)

does not depend on $kd$. The expression (22) is a special case of the more
general sum rule for $N$ qubit system.

$\sum_{i=1}^{N}\Gamma_{i}=N\Gamma$ (23)

where $\Gamma_{i}=2\mathrm{Re}(\lambda_{i})$. The sum rule (23) states that
there are no other losses in the system other than the coherent spontaneous
emission into a waveguide.

Consider now the solution of equations (17). By subtracting the third equation
in (17) from the first one we obtain:

$\frac{d}{dt}(\beta_{1}-\beta_{3})=-\frac{\Gamma}{2}\left(1-e^{i2kd}\right)(\beta_{1}-\beta_{3})$
(24)

It follows from (24) that if the edge qubits are initially not excited,
$\beta_{1}(0)-\beta_{3}(0)=0$, then this difference remains zero for all
times, $\beta_{1}(t)=\beta_{3}(t)\equiv\beta(t)$. This is quite reasonable
from symmetry consideration: if a central qubit is initially excited, the
temporal behavior of the amplitudes of the edge qubits must be the same.
Therefore, for this case, three equations (17) can be reduced to two
equations:

$\begin{array}[]{l}\frac{{d\beta_{2}}}{{dt}}=-\frac{\Gamma}{2}\beta_{2}(t)-\Gamma\beta(t)e^{ikd}\\\
\\\
\frac{{d\beta}}{{dt}}=-\frac{\Gamma}{2}\beta(t)\left({1+e^{2ikd}}\right)-\frac{\Gamma}{2}\beta_{2}(t)e^{ikd}\\\
\end{array}$ (25)

where $\beta_{2}(0)=1,\;\;\beta(0)=0$.

The characteristic roots of (25) are equal to $\lambda_{1}$ and $\lambda_{2}$
which are given in (21). Therefore, the solution of equations (25) reads:

$\begin{cases}\beta_{2}(t)&=b_{1}e^{\lambda_{1}t}+b_{2}e^{\lambda_{2}t},\\\
\beta(t)&=a_{1}e^{\lambda_{1}t}+a_{2}e^{\lambda_{2}t},\end{cases}$ (26)

where from initial conditions

$\begin{cases}&b_{1}+b_{2}=1,\\\ &a_{1}+a_{2}=0.\end{cases}$ (27)

Another two conditions follow from (25) for time derivatives at $t=0$.

$\begin{cases}&b_{1}\lambda_{1}+b_{2}\lambda_{2}=-\cfrac{\Gamma}{2},\\\
&a_{1}\lambda_{1}+a_{2}\lambda_{2}=-\cfrac{\Gamma}{2}e^{ikd}.\end{cases}$ (28)

From (27) and (28) we obtain

$b_{1}=-\frac{{\frac{\Gamma}{2}+\lambda_{2}}}{{\lambda_{1}-\lambda_{2}}};\;b_{2}=\frac{{\frac{\Gamma}{2}+\lambda_{1}}}{{\lambda_{1}-\lambda_{2}}}$
(29)

$a_{1}=-\frac{\Gamma}{2}\frac{{e^{ikd}}}{{\lambda_{1}-\lambda_{2}}};\;a_{2}=\frac{\Gamma}{2}\frac{{e^{ikd}}}{{\lambda_{1}-\lambda_{2}}}$
(30)

Using the explicit expressions (21) we obtain

$\begin{array}[]{l}{{a}}_{\rm{1}}{\rm{=}}\frac{1}{{{R}}}{{;a}}_{\rm{2}}=-\frac{1}{{R}}\\\
\\\
{{b}}_{\rm{1}}{\rm{=}}\frac{{{{R}}-{\rm{e}}^{{\rm{ikd}}}}}{{{{2R}}}}{{;b}}_{\rm{2}}{\rm{=}}\frac{{{{R+e}}^{{\rm{ikd}}}}}{{{{2R}}}}\\\
\end{array}$ (31)

where $R=\sqrt{e^{i2kd}+8}$.

Now we assume that the first qubit in the chain is initially excited.

$\beta_{1}(0)=1,\;\beta_{2}(0)=0,\;\beta_{3}(0)=0.$ (32)

In this case, all amplitudes behave differently, so that the equations (17)
should be used from which the solution can be straightforwardly obtained

$\begin{array}[]{l}{\rm{\beta}}_{\rm{1}}{\rm{(t)=}}\frac{{{{b}}_{\rm{2}}}}{{\rm{2}}}e^{\lambda_{1}t}+\frac{{b_{1}}}{2}e^{\lambda_{2}t}+\frac{1}{2}e^{{}_{\lambda_{3}t}}\\\
\\\
{\rm{\beta}}_{\rm{2}}{\rm{(t)=}}\frac{1}{{{R}}}\left({{\rm{e}}^{{\rm{\lambda}}_{\rm{1}}{\rm{t}}}-{\rm{e}}^{{\rm{\lambda}}_{\rm{2}}{\rm{t}}}}\right)\\\
\\\
{\rm{\beta}}_{\rm{3}}{\rm{(t)=}}\frac{{{{b}}_{\rm{2}}}}{{\rm{2}}}e^{\lambda_{1}t}+\frac{{b_{1}}}{2}e^{\lambda_{2}t}-\frac{1}{2}e^{{}_{\lambda_{3}t}}\\\
\end{array}$ (33)

where $b_{1}$ and $b_{2}$ are given in (31). For special cases when $kd$ is
integer multiple of $\pi$ we may obtain from (26) and (33) very simple forms
for the temporal behavior of the qubits’ amplitudes.

###  IV.1 $kd=\pi n$

For this case, $\lambda_{1}=-\frac{3}{2}\Gamma$, $\lambda_{2}=\lambda_{3}=0$
if $n$ is even number, and $\lambda_{1}=0$, $\lambda_{2}=-\frac{3}{2}\Gamma$,
$\lambda_{3}=0$ if $n$ is odd number. The calculations show that no matter
which qibit in the chain is initially excited (central or edge qubit) its
evolution is the same.

$\beta_{exc}(t)=\frac{1}{3}e^{-3\Gamma t/2}+\frac{2}{3}$ (34)

The evolution of unexcited qubits are also the same within a phase factor of
$\pi$.

$\beta_{unexc}(t)=(-1)^{n}\left(\frac{1}{3}e^{-3\Gamma
t/2}-\frac{1}{3}\right)$ (35)

The expression (35) is valid for any $n$ if the central qubit is excited. If
the first qubit is excited and $n$ is the odd number, the amplitudes of
unexcited qubits evolve with the opposite phases:
$\beta_{2}(t)=\beta_{unexc}(t)$, $\beta_{3}(t)=-\beta_{unexc}(t)$.

It worth noting the equality of the amplitudes of the second and the third
qubit (if the first qubit is initially excited) even though they are located
at different distance from the excited qubit. This is the consequence of
Wigner-Weisskopf (or Markov) approximation: in this special case there are no
interference between qubits so that all unexcited qubits feel the photon field
simultaneously no matter how far they are from the excited qubit. From (34)
and (35) we find a full probability of photon emission (8)

$P_{ph}(t)=1-2|\beta_{unexc}(t)|^{2}-|\beta_{exc}(t)|^{2}=\frac{1}{3}\left(1-e^{-3\Gamma
t}\right)$ (36)

The probability of qubits amplitudes and a full probability of photon emission
are shown in Fig. 3

Figure 3: The probability amplitudes for $kd=n\pi$ of excited,
$|\beta_{exc}(t)|^{2}$ (the upper line) and unexcited $|\beta_{unexc}(t)|^{2}$
(the lower line) qubits. The probability of a full photon emission is shown
between the qubit lines.

As is seen from Fig. 3, as the time proceeds the qubits amplitudes become
”frozen” at some level. We see from (34) and (35) that this frozenness is
solely due to the dark states ($\lambda_{n}=0$), which prevent the qubits’
amplitudes from decaying to zero. The frozenness can be lifted if $kd$ is not
equal to an integer multiple of $\pi$. Then all qubits are damped to zero with
the rate being determined by the root of $\lambda_{n}$ which has the lowest
real part.

###  IV.2 $kd=\pi(2n+1)/2$

If n is even number then $\lambda_{1}=-\frac{\Gamma}{4}(1+i\sqrt{7})$,
$\lambda_{2}=-\frac{\Gamma}{4}(1-i\sqrt{7})$, $\lambda_{3}=-\Gamma$ , and if n
is odd number then $\lambda_{1}=-\frac{\Gamma}{4}(1-i\sqrt{7})$,
$\lambda_{2}=-\frac{\Gamma}{4}(1+i\sqrt{7})$, $\lambda_{3}=-\Gamma$

If the central qubit is initially excited we obtain:

$\begin{array}[]{l}\beta_{1}(t)=\beta_{3}(t)=(-1)^{n+1}i\frac{2}{{\sqrt{7}}}e^{-\frac{1}{4}\Gamma
t}\sin\left({\frac{{\sqrt{7}}}{4}\Gamma t}\right)\\\ \\\
\beta_{2}(t)=e^{-\frac{1}{4}\Gamma
t}\left({\cos\left({\frac{{\sqrt{7}}}{4}\Gamma
t}\right)+(-1)^{n}\frac{1}{{\sqrt{7}}}\sin\left({\frac{{\sqrt{7}}}{4}\Gamma
t}\right)}\right)\\\ \end{array}$ (37)

For initially excited edge qubit we obtain:

$\begin{array}[]{l}\beta_{1}(t)=\frac{{e^{-\frac{1}{4}\Gamma
t}}}{2}\left({\cos\left({\frac{{\sqrt{7}}}{4}\Gamma
t}\right)+\frac{1}{{\sqrt{7}}}\sin\left({\frac{{\sqrt{7}}}{4}\Gamma
t}\right)}\right)+\frac{1}{2}e^{-\Gamma t}\\\ \\\
\beta_{2}(t)=(-1)^{n+1}i\frac{2}{{\sqrt{7}}}e^{-\frac{1}{4}\Gamma
t}\sin\left({\frac{{\sqrt{7}}}{4}\Gamma t}\right)\\\ \\\
\beta_{3}(t)=\frac{{e^{-\frac{1}{4}\Gamma
t}}}{2}\left({\cos\left({\frac{{\sqrt{7}}}{4}\Gamma
t}\right)+\frac{1}{{\sqrt{7}}}\sin\left({\frac{{\sqrt{7}}}{4}\Gamma
t}\right)}\right)-\frac{1}{2}e^{-\Gamma t}\\\ \end{array}$ (38)

The oscillatory behavior of the qubits amplitudes in these expressions is a
signature of the vacuum Rabi oscillations with the frequency
$\sqrt{7}\Gamma/4$.

The evolution of qubits amplitudes for half-integer multiple of $\pi$ with the
central qubit being initially excited is shown in Fig. 4.

Figure 4: The evolution of qubits amplitudes for $kd$ equals the half-integer
multiple of $\pi$ with the central qubit being initially excited.

The evolution of $P_{ph}(t)$ reveals two clear visible steps where,
$dP_{ph}(t)/dt=0$. These steps can be attributed to the interrelation between
the temporal dynamics of the different amplitudes. As is seen from this
figure, the steps on $P_{ph}(t)$ curve are in the vicinity of extremum points
of qubits amplitudes. Physically, these steps are the signature of the
trapping of the photon radiation. As is seen in Fig. 3, near the first step
the radiation emitted by excited qubit is absorbed by unexcited qubits, so
that the rate of the output radiation is not changed. The evolution of qubits’
amplitudes for half-integer multiple of $\pi$ with the first qubit being
initially excited is shown in Fig. 5.

Figure 5: The evolution of qubits amplitudes for $kd$ equals the half-integer
multiple of $\pi$ with the first qubit being initially excited

Here, we see that the evolution of qubits amplitude depends on their distance
from the excited qubit. The greater is the distance of a qubit from the
excited one, the less it is affected by the excitation.

##  V Spectral density of photon radiation

The quantity $\gamma_{k}(t)$ in (5) allows for the calculation of a spectral
density (6) of spontaneous emission into a waveguide, $S(\omega,t)$. The
equation (10) for identical qubits reads:

$\gamma(\omega,t)=-ig_{k}\sum\limits_{n=1}^{3}{}e^{-ikx_{n}}\int\limits_{0}^{t}{\beta_{n}(t^{\prime})e^{i(\omega-\Omega)t^{\prime}}dt^{\prime}}$
(39)

Using the expression (18) for $\beta_{n}(t)$ we obtain

$\gamma(\omega,t)=g_{k}\sum\limits_{n,m=1}^{3}{}a_{j}^{(n)}e^{-ikx_{n}}\frac{{1-e^{i(\omega-\Omega-i\lambda_{m})t}}}{{\omega-\Omega-i\lambda_{m}}}$
(40)

As $g_{k}$ in (39) and (40) depends on the waveguide length, $L$ and the light
velocity $v_{g}$ (see (108)), we factorize $\gamma_{k}(t)$ as follows:

$\gamma_{k}(t)\equiv\left({\frac{{v_{g}}}{{2L\Omega}}}\right)^{1/2}\left({\frac{\Gamma}{\Omega}}\right)^{1/2}f_{k}(t)$

Therefore, we define spectral density in the following form:

$S(\omega,t)=\frac{{\left|{\gamma_{k}(t)}\right|^{2}}}{{\left({\frac{{v_{g}}}{{2L\Omega}}}\right)}}=\frac{\Gamma}{\Omega}\left|{f_{k}(t)}\right|^{2}$
(41)

The equation (40) is a rather general expression. As the time proceeds only
the term with the lowest real part of $\lambda_{m}$ survives. The concise
analytical results can be obtained only for several simple cases. Below we
calculate the spectral density for $kd=n\pi$, using the amplitudes (34) and
(35) in (39).

When a central qubit is initially excited we obtain for $\gamma_{k}(t)$ the
following result:

$\begin{array}[]{l}\gamma_{k}(t)=-g_{k}\frac{2}{3}\left({1-(-1)^{n}\cos\left({\frac{\omega}{\Omega}n\pi}\right)}\right)\frac{{e^{i(\omega-\Omega)t}-1}}{{\omega-\Omega}}\\\
\\\
-g_{k}\frac{1}{3}\left({1+(-1)^{n}2\cos\left({\frac{\omega}{\Omega}n\pi}\right)}\right)\frac{{e^{i\left({\omega-\Omega+i\frac{{3\Gamma}}{2}}\right)t}-1}}{{\left({\omega-\Omega+i\frac{{3\Gamma}}{2}}\right)}}\\\
\end{array}$ (42)

If the first qubit is initially excited we obtain:

$\begin{array}[]{l}\gamma_{k}(t)=-g_{k}\frac{1}{3}\left({2e^{i\left({\frac{\omega}{\Omega}n\pi}\right)}-e^{-i\left({\frac{\omega}{\Omega}n\pi}\right)}-(-1)^{n}}\right)\frac{{e^{i(\omega-\Omega)t}-1}}{{\omega-\Omega}}\\\
\\\
-g_{k}\frac{1}{3}\left({(-1)^{n}+2\cos\left({\frac{\omega}{\Omega}n\pi}\right)}\right)\frac{{e^{i\left({\omega-\Omega+i\frac{{3\Gamma}}{2}}\right)t}-1}}{{\left({\omega-\Omega+i\frac{{3\Gamma}}{2}}\right)}}\\\
\end{array}$ (43)

It worth noting that in expressions (42),(43) a wave vector $k$ is related to
a running frequency $\omega$: $k=\omega/v_{g}$. The first lines in (42) and
(43) are related to the dark states ($\lambda_{n}=0$). They do not lead to a
singularity near resonance $\omega\approx\Omega$, since $kd$ \- dependent
prefactors in this terms tend to zero more rapidly than the denominator does.
At the point of a resonance these terms exactly equal to zero. A physical
reason for this is that the dark states do not interact with a photon field
and, therefore, cannot contribute to the photon emission. Therefore, near the
resonance, a spectral density of spontaneous emission can be approximated by a
Lorentzian form with a full width at half the height of the resonance line
being equal to $3\Gamma$.

$\gamma_{k}(t)\approx
g_{k}\frac{{e^{i\left({\omega-\Omega+i\frac{{3\Gamma}}{2}}\right)t}-1}}{{\left({\omega-\Omega+i\frac{{3\Gamma}}{2}}\right)}}\xrightarrow[t\to\infty]{}\frac{{g_{k}}}{{\left({\omega-\Omega+i\frac{{3\Gamma}}{2}}\right)}}$
(44)

$S(\omega,t\to\infty)=\frac{\Gamma\Omega}{{(\omega-\Omega)^{2}+\left({\frac{{3\Gamma}}{2}}\right)^{2}}}$
(45)

For $kd=(2n+1)\pi/2$, the qubits amplitudes are given in equations (37) and
(38). The calculation of (39) with these amplitudes results in the following
expressions for $\gamma_{k}(t)$ with the central qubit being initially
excited:

$\begin{gathered}\gamma_{k}(t)=-g_{k}\left[{\frac{1}{2}\left({1-\frac{{(-1)^{n}i}}{{\sqrt{7}}}}\right)-\frac{{2(-1)^{n}}}{{\sqrt{7}}}\cos
kd}\right]\hfill\\\
\times\frac{{e^{i(\omega-\Omega+\frac{{\sqrt{7}}}{4}\Gamma+i\frac{1}{4}\Gamma)t}-1}}{{\omega-\Omega+\frac{{\sqrt{7}}}{4}\Gamma+i\frac{1}{4}\Gamma}}\hfill\\\
-g_{k}\left[{\frac{1}{2}\left({1+\frac{{(-1)^{n}i}}{{\sqrt{7}}}}\right)+\frac{{2(-1)^{n}}}{{\sqrt{7}}}\cos
kd}\right]\hfill\\\
\times\frac{{e^{i(\omega-\Omega-\frac{{\sqrt{7}}}{4}\Gamma+i\frac{1}{4}\Gamma)t}-1}}{{\omega-\Omega-\frac{{\sqrt{7}}}{4}\Gamma+i\frac{1}{4}\Gamma}}\hfill\\\
\end{gathered}$ (46)

If the edge qubit is initially excited we obtain:

$\begin{gathered}\gamma_{k}(t)=-ig_{k}\left[{\frac{1}{2}\left({1-\frac{i}{{\sqrt{7}}}}\right)\cos
kd+\frac{{(-1)^{n}}}{{\sqrt{7}}}}\right]\hfill\\\
\times\frac{{e^{i(\omega-\Omega+\frac{{\sqrt{7}}}{4}\Gamma+i\frac{1}{4}\Gamma)t}-1}}{{\omega-\Omega+\frac{{\sqrt{7}}}{4}\Gamma+i\frac{1}{4}\Gamma}}\hfill\\\
-ig_{k}\left[{\frac{1}{2}\left({1+\frac{i}{{\sqrt{7}}}}\right)\cos
kd-\frac{{(-1)^{n}}}{{\sqrt{7}}}}\right]\hfill\\\
\times\frac{{e^{i(\omega-\Omega-\frac{{\sqrt{7}}}{4}\Gamma+i\frac{1}{4}\Gamma)t}-1}}{{\omega-\Omega-\frac{{\sqrt{7}}}{4}\Gamma+i\frac{1}{4}\Gamma}}\hfill\\\
-ig_{k}\sin
kd\frac{{e^{i(\omega-\Omega+i\frac{1}{2}\Gamma)t}-1}}{{\omega-\Omega+i\frac{1}{2}\Gamma}}\hfill\\\
\end{gathered}$ (47)

In expressions (46), (47) $kd=\frac{\omega}{\Omega}(2n+1)\pi/2$ .

The spectral density of photon radiation (41) calculated from (46) and (47)
for $n=1,kd=1.5\pi,t=20/\Gamma,\Gamma/\Omega=10^{-3}$ is shon in Fig. 6.

Figure 6: Spectral density for three-qubit system calculated from (46) and
(47) for $n=1,kd=1.5\pi,t=20/\Gamma,\Gamma/\Omega=10^{-3}$. (a) second qubit
is initially excited, (b) first qubit is initially excited.

Two peaks at these plots are a clear signature of vacuum Rabi oscillations..

##  VI Effective non-Hermitian Hamiltonian and collective states

The elimination of photon variables allows us to express the photon-mediated
interaction between identical qubits in terms of a non-Hermitian effective
Hamiltonian, which in the Markovian approximation reads Chang2012 ; Lehm1970 :

$H_{eff}=-i\frac{\Gamma}{2}\sum\limits_{m,n=1}^{3}{}e^{ik\left|{x_{m}-x_{n}}\right|}\sigma_{m}^{+}\sigma_{n}$
(48)

where $k=\Omega/v_{g}$, $v_{g}$ is the velocity of electromagnetic wave in a
waveguide, $x_{n}$ is the position of $n$-th qubit,
$\sigma^{+}_{n},\sigma_{n}$ are raising and lowering spin operators for $n$-th
qubit.

The rate of spontaneous emission $\Gamma$ of an individual qubit is defined by
the Fermi golden rule:

$\Gamma=2\pi\sum\limits_{k}{\left|{g_{k}}\right|^{2}\delta(\omega_{k}-\Omega)}$
(49)

It follows from (48) that the photon-mediated interaction between qubits in
such a system results in the coherent
$J_{mn}=\Gamma\sin\left({k\left|{x_{m}-x_{n}}\right|}\right)/2$ and
dissipative $\Gamma_{mn}=\Gamma\cos\left({k\left|{x_{m}-x_{n}}\right|}\right)$
rates. The coherent rate shifts the positions of the qubits resonances, while
the dissipative rate gives rise to the additional spontaneous emission into
the waveguide mode. Unlike the real atoms with short-range dipole-dipole
interaction, here a coherent interaction $J_{mn}$ is a long-range one: every
qubit is sensitive to its distant neighbor.

The wave function for Hamiltonian (48) can be expressed in terms of a
superposition of the single excited states

$\Psi(t)=\sum\limits_{n=1}^{3}\beta_{n}(t)\left|n\right\rangle$ (50)

Even though the wave functions (5) and (50) are different, the qubits
amplitudes $\beta_{n}(t)$ in these expressions are the same quantities‘.
Indeed, the equations for $\beta_{n}(t)$ in (50) can be derived from the time-
dependent Schrodinger equation ${{id\Psi/dt=H}}_{{{eff}}}{{\Psi}}$ :

$\frac{{d\beta_{n}}}{{dt}}=-\frac{\Gamma}{2}\sum\limits_{m=1}^{3}{}\beta_{m}(t)e^{ik\left|{x_{m}-x_{n}}\right|},\quad(n=1,2,3)$
(51)

The equations (51) are nothing but a set of equations (17). Therefore, within
a single excitation subspace two Hamiltonians, (1) and (48) are equivalent in
that they provide the same equations for the qubits amplitudes.

In a single-excitation subspace, the Hamiltonian (48) has three collective
eigenfunctions:

$\left|{\Psi_{i}(t)}\right\rangle=e^{-i\bar{E}_{i}t}\sum\limits_{n=1}^{3}{}\alpha_{n}^{(i)}\left|n\right\rangle;\;(i=1,2,3)$
(52)

where $\bar{E}_{i}$ is a complex energy

$\bar{E}_{i}=E_{i}-i\frac{{\Gamma_{i}}}{2}$ (53)

The quantities $E_{i}$ and $\Gamma_{i}$ depend on the system parameters
$\Gamma,k,x_{n}$. The wave vectors with $\Gamma_{i}<\Gamma$ are called
subradiant states, those with $\Gamma_{i}>\Gamma$ are called superradiant
states. The coefficients in (52) can be obtained from the Schrodinger equation
$H_{eff}\Psi=E\Psi$.

$E_{j}\sum\limits_{n=1}^{3}{}\alpha_{n}^{(j)}\left|n\right\rangle=-i\frac{\Gamma}{2}\sum\limits_{m,n=1}^{3}{}e^{ikd\left|{m-n}\right|}\alpha_{n}^{(j)}\left|m\right\rangle;j=1,2,3$
(54)

Complex energies can be found by equating to zero the determinant of the
matrix:

$\left({E+i\frac{\Gamma}{2}}\right)\delta_{mn}+i\frac{\Gamma}{2}e^{ikd\left|{m-n}\right|}\left({1-\delta_{mn}}\right);\quad(m,n=1,2,3)$
(55)

A comparison between the determinants of (19) and (55) shows that the
quantities $\lambda_{n}$ are related to those of complex energies,
$\bar{E}_{n}$ : $\bar{E}_{n}=i\lambda_{n}$ where $\lambda_{n}$ are given in
(21).

From (54) we obtain three equations for the coefficients, .

$\begin{gathered}\left({E_{j}+i\frac{\Gamma}{2}}\right)\alpha_{1}^{(j)}+i\frac{\Gamma}{2}e^{ikd}\alpha_{2}^{(j)}+i\frac{\Gamma}{2}e^{2ikd}\alpha_{3}^{(j)}=0\hfill\\\
\left({E_{j}+i\frac{\Gamma}{2}}\right)\alpha_{2}^{(j)}+i\frac{\Gamma}{2}e^{ikd}\alpha_{1}^{(j)}+i\frac{\Gamma}{2}e^{ikd}\alpha_{3}^{(j)}=0\hfill\\\
\left({E_{j}+i\frac{\Gamma}{2}}\right)\alpha_{3}^{(j)}+i\frac{\Gamma}{2}e^{ikd}\alpha_{2}^{(j)}+i\frac{\Gamma}{2}e^{2ikd}\alpha_{1}^{(j)}=0\hfill\\\
\end{gathered}$ (56)

where for simplicity we remove the bar over $E_{i}$.

Because Hamiltonian (48) is non-Hermitian, the eigenfunctions (52) are neither
normalized nor orthonormal. It is known that a correct calculation of the
coefficients $\alpha_{n}^{(i)}$ in (52) requires a bi-orthogonal set of
eigenfunctions $\left|{\bar{\Psi}_{i}(t)}\right\rangle$ which are a solution
of the Schrodinger equation for $H_{eff}^{\dagger}$. In our case
$H_{eff}^{\dagger}=H_{eff}^{*}$ with the consequence that the complex
conjugate of an eigenstate $\left|{\Psi_{i}(t)}\right\rangle$ of $H_{eff}$ is
an eigenstate of $H_{eff}^{\dagger}$. Therefore, the conditions for
normalization and orthonormality between eigenfunctions of these two sets lead
to the following equations for the coefficients Keck2003 ; Brody2014 .

$\sum\limits_{n=1}^{3}{}\left({\alpha_{n}^{(i)}}\right)^{2}=1,\quad i=1,2,3$
(57)

$\sum\limits_{n=1}^{3}{}\alpha_{n}^{(i)}\alpha_{n}^{(j)}=0,\quad i\neq j,\quad
i,j=1,2,3$ (58)

It worth noting that the coefficients $\alpha_{n}^{(i)}$ in (57) and (58) are
in general complex quantities. In what follows we use equations (56) and the
conditions (57) and (58) for the calculation of the coefficients
$\alpha_{n}^{(i)}$.

As the energies $E_{j}$ are obtained from the determinant of equations (56),
these equations are not independent. Therefore, for the calculations of the
coefficients $\alpha_{n}^{(i)}$ we may take any two of them. For subsequent
calculations we take first and second equations in (56). First, we calculate
the coefficients $\alpha_{n}^{(1)}$ . With $E_{1}=i\lambda_{1}$ determined in
(21) two first equations in (56) read as follows:

$\begin{gathered}\left({e^{ikd}+R}\right)\alpha_{1}^{(1)}-2\alpha_{2}^{(1)}-2e^{ikd}\alpha_{3}^{(1)}=0\hfill\\\
\left({e^{ikd}+R}\right)\alpha_{2}^{(1)}-2\alpha_{1}^{(1)}-2\alpha_{3}^{(1)}=0\hfill\\\
\end{gathered}$ (59)

where $R$ is given in (31).

Additional third equation is given by the normalizing condition (57).
Therefore, for the wave-function $\left|{\Psi_{1}(t)}\right\rangle$ we find
the following set of the coefficients $\alpha_{n}^{(1)}$ .

$\alpha_{1}^{(1)}=\alpha_{3}^{(1)}=\pm\frac{{e^{ikd}\left({e^{ikd}+R}\right)+2}}{D(R)}$
(60)

$\alpha_{2}^{(1)}=\pm\frac{{3e^{ikd}+R}}{D(R)}$ (61)

where

$D(R)=\left({4e^{4ikd}+34e^{2ikd}+14e^{ikd}R+4e^{3ikd}R+16}\right)^{1/2}$ (62)

Similar calculations for $E_{2}=i\lambda_{2}$ provide the result that differs
from (60), (61) only by the sign of $R$. Therefore, for the wave-function
$\left|{\Psi_{2}(t)}\right\rangle$ we find the following set of the
coefficients $\alpha_{n}^{(2)}$ .

$\alpha_{1}^{(2)}=\alpha_{3}^{(2)}=\pm\frac{{e^{ikd}\left({e^{ikd}-R}\right)+2}}{D(-R)}$
(63)

$\alpha_{2}^{(2)}=\pm\frac{{3e^{ikd}-R}}{D(-R)}$ (64)

The signs in right hand side of equations (60), (61), and (63), (64) must be
the same (plus or minus) for all four expressions. This requirement follows
from the orthonormality condition (58):

$\alpha_{1}^{(1)}\alpha_{1}^{(2)}+\alpha_{2}^{(1)}\alpha_{2}^{(2)}+\alpha_{3}^{(1)}\alpha_{3}^{(2)}=0$
(65)

Finally, for $E_{3}=i\lambda_{3}$ we obtain two equations

$\begin{gathered}\alpha_{1}^{(3)}+e^{-ikd}\alpha_{2}^{(3)}+\alpha_{3}^{(3)}=0\hfill\\\
\alpha_{1}^{(3)}+e^{ikd}\alpha_{2}^{(3)}+\alpha_{3}^{(3)}=0\hfill\\\
\end{gathered}$ (66)

which provide with account for the normalizing condition (57) the following
result:

$\alpha_{1}^{(3)}=\pm\frac{1}{{\sqrt{2}}};\alpha_{3}^{(3)}=\mp\frac{1}{{\sqrt{2}}};\alpha_{2}^{(3)}=0;\;$
(67)

Therefore, the coefficients of collective states
$\left|{\Psi_{1}(t)}\right\rangle$ and $\left|{\Psi_{2}(t)}\right\rangle$
depend on $kd$, while those of $\left|{\Psi_{3}(t)}\right\rangle$ are $kd$\-
independent.

Below, we consider special cases. For $kd=2\pi n$, where $n$ is integer, we
obtain from (60), (61)
$\alpha_{1}^{(1)}=\alpha_{2}^{(1)}=\alpha_{3}^{(1)}=\pm\frac{1}{{\sqrt{3}}}$.
However, the calculation of the coefficients $\alpha_{n}^{(2)}$ from (63),
(64) is not straightforward: at this point both the numerator and the
denominator in these equations are equal to zero. In order to resolve this
uncertainty we put in these equations $kd=2\pi n+\epsilon$, where $\epsilon$
tends to zero. In this case, both the numerator and the denominator tend to
zero as $\epsilon$, and their ratio is finite. The calculations show that as
$kd$ tends to $2\pi n$ the coefficients $\alpha_{n}^{(2)}$ tend to their
finite values:
$\alpha_{1}^{(2)}=\alpha_{3}^{(2)}=\frac{1}{{\sqrt{6}}};\;\alpha_{2}^{(2)}=-\frac{2}{{\sqrt{6}}}$.

For $kd=(2n+1)\pi$ the picture is vise versa:
$\alpha_{1}^{(1)}=\alpha_{3}^{(1)}=\frac{1}{{\sqrt{6}}};\;\alpha_{2}^{(1)}=-\frac{2}{{\sqrt{6}}}$,
while
$\alpha_{1}^{(2)}=\alpha_{2}^{(2)}=\alpha_{3}^{(2)}=\pm\frac{1}{{\sqrt{3}}}$.

For $kd=(2n+1)\pi/2$ we obtain from (60), (61), and (63), (64):

$\begin{gathered}\alpha_{1}^{(1)}=\alpha_{3}^{(1)}=\frac{{-{\text{i}}+(-1)^{n}\sqrt{7}}}{{\sqrt{2}\left({7-{\text{i}}(-1)^{n}5\sqrt{7}}\right)^{1/2}}}\hfill\\\
\alpha_{2}^{(1)}=\frac{{(-1)^{n}3-{\text{i}}\sqrt{7}}}{{\sqrt{2}\left({7-{\text{i}}(-1)^{n}5\sqrt{7}}\right)^{1/2}}}\hfill\\\
\end{gathered}$ (68)

$\begin{gathered}\alpha_{1}^{(2)}=\alpha_{3}^{(2)}=\frac{{-{\text{i}}-(-1)^{n}\sqrt{7}}}{{\sqrt{2}\left({7+i(-1)^{n}5\sqrt{7}}\right)^{1/2}}}\hfill\\\
\alpha_{2}^{(2)}=\frac{{(-1)^{n}3+i\sqrt{7}}}{{\sqrt{2}\left({7+i(-1)^{n}5\sqrt{7}}\right)^{1/2}}}\hfill\\\
\end{gathered}$ (69)

Below we summarize these results in the explicit forms of collective wave
functions for above special cases. If $kd$ is equal to integer multiple of
$\pi$ we obtain:

$\begin{gathered}\left|{\Psi_{1`}(t)}\right\rangle=e^{-\frac{3}{2}\Gamma
t}\frac{1}{{\sqrt{3}}}\left({\left|1\right\rangle+\left|2\right\rangle+\left|3\right\rangle}\right)\hfill\\\
\left|{\Psi_{2`}(t)}\right\rangle=\frac{1}{{\sqrt{6}}}\left({\left|1\right\rangle-2\left|2\right\rangle+\left|3\right\rangle}\right)\hfill\\\
\left|{\Psi_{3`}(t)}\right\rangle=\frac{1}{{\sqrt{2}}}\left({\left|1\right\rangle-\left|3\right\rangle}\right)\hfill\\\
\end{gathered}$ (70)

Therefore, two of these collective states $\left|{\Psi_{2}(t)}\right\rangle$
and $\left|{\Psi_{3}(t)}\right\rangle$ are dark as their decay widths are
zero.

For $kd=(2n+1)\pi/2$ the collective wave functions are as follows:

$\begin{gathered}\left|{\Psi_{1`}(t)}\right\rangle=\frac{{e^{-\frac{1}{4}\Gamma
t}e^{-i\frac{{\sqrt{7}}}{4}(-1)^{n}\Gamma
t}}}{{\sqrt{2}\left({7-i(-1)^{n}5\sqrt{7}}\right)^{1/2}}}\left({(-i+(-1)^{n}\sqrt{7})\left({\left|1\right\rangle+\left|3\right\rangle}\right)+((-1)^{n}3-i\sqrt{7})\left|2\right\rangle}\right)\hfill\\\
\left|{\Psi_{2`}(t)}\right\rangle=-\frac{{e^{-\frac{1}{4}\Gamma
t}e^{i\frac{{\sqrt{7}}}{4}(-1)^{n}\Gamma
t}}}{{\sqrt{2}\left({7+i(-1)^{n}5\sqrt{7}}\right)^{1/2}}}\left({(i+(-1)^{n}\sqrt{7})\left({\left|1\right\rangle+\left|3\right\rangle}\right)-((-1)^{n}3+i\sqrt{7})\left|2\right\rangle}\right)\hfill\\\
\left|{\Psi_{3`}(t)}\right\rangle=e^{-\Gamma
t}\frac{1}{{\sqrt{2}}}\left({\left|1\right\rangle-\left|3\right\rangle}\right)\hfill\\\
\end{gathered}$ (71)

where $n$ is any integer.

It is instructive to rewrite (71) in terms of the dark, $|D\rangle$ and
bright, $|B\rangle$ states for a _two qubits_ (the first qubit and the third
one) separated by $\lambda/2$. For $n=0$ we obtain from (71):

$\begin{gathered}\left|{\Psi_{1`}(t)}\right\rangle=\frac{{e^{-\frac{1}{4}\Gamma
t}e^{-i\frac{{\sqrt{7}}}{4}\Gamma
t}}}{{\sqrt{2}\left({7-i5\sqrt{7}}\right)^{1/2}}}\left({(-i+\sqrt{7})\sqrt{2}\left|D\right\rangle\otimes\left|{g_{2}}\right\rangle+(3-i\sqrt{7})\left|G\right\rangle\otimes\left|{e_{2}}\right\rangle}\right)\hfill\\\
\left|{\Psi_{2`}(t)}\right\rangle=-\frac{{e^{-\frac{1}{4}\Gamma
t}e^{i\frac{{\sqrt{7}}}{4}\Gamma
t}}}{{\sqrt{2}\left({7+i5\sqrt{7}}\right)^{1/2}}}\left({(i+\sqrt{7})\sqrt{2}\left|D\right\rangle\otimes\left|{g_{2}}\right\rangle-(3+i\sqrt{7})\left|G\right\rangle\otimes\left|{e_{2}}\right\rangle}\right)\hfill\\\
\left|{\Psi_{3`}(t)}\right\rangle=e^{-\Gamma
t}\left|B\right\rangle\otimes\left|{g_{2}}\right\rangle\hfill\\\
\end{gathered}$ (72)

where $|G\rangle=|g_{1}g_{3}\rangle$ and

$\left|D\right\rangle=\frac{{\left|{e_{1}g_{3}}\right\rangle+\left|{g_{1}e_{3}}\right\rangle}}{{\sqrt{2}}};\;\left|B\right\rangle=\frac{{\left|{e_{1}g_{3}}\right\rangle-\left|{g_{1}e_{3}}\right\rangle}}{{\sqrt{2}}}$
(73)

As is seen from the third equation in (72) the bright state of the two-qubit
system decays independently on the presence of the second (central) qubit.
However, the decay of the dark state can be revealed only through its
entanglement with the second qubit Mirho2019 . In fact, the interaction
between second qubit and the dark state formed by two edge qubits gives rise
to the vacuum Rabi oscillations which are shown in (37), (38), and in Fig.4,
Fig.5.

###  VI.1 Relation between qubits amplitudes and collective states

Here, we show how the qubits amplitudes $\beta_{n}(t)$ are related to
collective wave functions $\left|{\Psi_{i}(t)}\right\rangle$ (52). We write
the dynamic wave function $\left|{\Psi(t)}\right\rangle$ (50) as a
decomposition over the collective states (52).

$\left|{\Psi(t)}\right\rangle=\sum\limits_{i=1}^{3}{A_{i}\left|{\Psi_{i}(t)}\right\rangle}=\sum\limits_{i,n=1}^{3}{A_{i}e^{-i\bar{E}_{i}t}\alpha_{n}^{(i)}\left|n\right\rangle}$
(74)

From (50) we obtain:

$\beta_{n}(t)=\sum\limits_{i=1}^{3}{e^{-i\bar{E}_{i}t}A_{i}\alpha_{n}^{(i)}}$
(75)

with the initial conditions

$\begin{gathered}\beta_{n_{0}}(0)=\sum\limits_{i=1}^{3}{A_{i}\alpha_{n_{0}}^{(i)}}=1\\\
\beta_{n}(0)=\sum\limits_{i=1}^{3}{A_{i}\alpha_{n}^{(i)}}=0;\;n\neq
n_{0}\end{gathered}$ (76)

where $n_{0}$ is the sequence number of excited qubit.

The probability amplitude $|\beta_{n}(t)|^{2}$ can directly be expressed in
terms of collective state wave functions $|\Psi_{i}(t)\rangle$.

$\left|{\beta_{n}(t)}\right|^{2}=\left\langle
n\right|\left({\sum\limits_{i,j=1}^{3}{}A_{i}A_{j}^{*}\left|{\Psi_{i}}\right\rangle\left\langle{\Psi_{j}}\right|}\right)\left|n\right\rangle$
(77)

From linear algebraic equations (76) we can find coefficients $A_{i}$ and,
therefore, restore the qubits’ amplitudes $\beta_{n}(t)$. However, if the
number of qubits is large, this procedure is not convenient for computer
simulations. The main reason is that it requires first, the calculation of $N$
complex energies from determinant of $N\times N$ matrix analogous to (55),
second, the calculations of $\alpha_{n}^{(i)}$ from non linear conditions (57)
and (58), and third, the solution of a system of $N$ linear algebraic
equations analogous to (76). Every of these three steps is not simple from a
mathematical point of view. It is more convenient to directly compute the
qubits amplitudes $\beta_{n}(t)$ from a set of the linear differential
equations (16), which allow us to completely avoid all three steps we
mentioned above.

Nevertheless, we should like to mention some interesting consequences that
follow from equations (75), (76), (77). First, from (75) we see that the dark
states ($Im\bar{E}_{i}=0$) contribute to $\beta_{n}(t)$, even though they do
not contribute to the spectrum of the photon emission. Second, it follows from
(77) that if specific collective state $|\Psi_{i}(t)\rangle$ does not contain
qubit state $|n\rangle$, then this collective state does not take part in the
formation of the dynamics of the qubits amplitudes $\beta_{n}(t)$. As was
shown above (see (67)) the state $|\Psi_{3}(t)\rangle$ does not contain the
qubit state $|2\rangle$. Therefore, independently on the value of $kd$, only
two collective states $\left|{\Psi_{1}(t)}\right\rangle$ and
$\left|{\Psi_{2}(t)}\right\rangle$ take part in the formation of the dynamics
of the qubits amplitudes $\beta_{2}(t)$.

##  VII Three non-identical qubits

Here, we consider a system in which all three qubits are identical except for
the frequency of a second qubit which has a different value $\Omega_{0}$. For
this case, we rewrite the equations (12), (13), and (14).

$\begin{gathered}\frac{{d\beta_{1}}}{{dt}}=-\frac{\Gamma}{2}\beta_{1}(t)-i\frac{{\delta\Omega}}{2}\beta_{1}(t)\hfill\\\
-\beta_{2}(t)\frac{1}{2}\left({\frac{{\Omega_{0}}}{\Omega}}\right)^{1/2}\Gamma
e^{ik_{0}d}-\frac{\Gamma}{2}\beta_{3}(t)e^{2ikd}\hfill\\\ \end{gathered}$ (78)

$\begin{gathered}\frac{{d\beta_{2}}}{{dt}}=-\frac{\Gamma}{2}\beta_{2}(t)+i\frac{{\delta\Omega}}{2}\beta_{2}(t)\hfill\\\
-\frac{1}{2}\left({\frac{\Omega}{{\Omega_{0}}}}\right)^{1/2}\Gamma
e^{ikd}\left({\beta_{1}(t)+\beta_{3}(t)}\right)\hfill\\\ \end{gathered}$ (79)

$\begin{gathered}\frac{{d\beta_{3}}}{{dt}}=-\frac{\Gamma}{2}\beta_{3}(t)-i\frac{{\delta\Omega}}{2}\beta_{3}(t)-\beta_{1}(t)\frac{\Gamma}{2}e^{2ikd}\hfill\\\
-\beta_{2}(t)\frac{1}{2}\left({\frac{{\Omega_{0}}}{\Omega}}\right)^{1/2}\Gamma
e^{ik_{0}d}\hfill\\\ \end{gathered}$ (80)

where $\delta\Omega=\Omega-\Omega_{0}$,
$k=\frac{\Omega}{{v_{g}}},\;k_{0}=\frac{{\Omega_{0}}}{{v_{g}}}$

Characteristic roots can be found from the determinant of equations (78),
(79), and (80).

$\begin{gathered}\lambda_{1}=-\frac{\Gamma}{2}\left({1+\frac{1}{2}e^{2ikd}}\right)\hfill\\\
-\frac{\Gamma}{4}e^{ikd}\sqrt{e^{2ikd}+8e^{i(k_{0}-k)d}+4i\frac{{\delta\Omega}}{\Gamma}-4\left({\frac{{\delta\Omega}}{\Gamma}}\right)^{2}e^{-2ikd}}\hfill\\\
\end{gathered}$ (81)

$\begin{gathered}\lambda_{2}=-\frac{\Gamma}{2}\left({1+\frac{1}{2}e^{2ikd}}\right)\hfill\\\
+\frac{\Gamma}{4}e^{ikd}\sqrt{e^{2ikd}+8e^{i(k_{0}-k)d}+4i\frac{{\delta\Omega}}{\Gamma}-4\left({\frac{{\delta\Omega}}{\Gamma}}\right)^{2}e^{-2ikd}}\hfill\\\
\end{gathered}$ (82)

$\lambda_{3}=\frac{\Gamma}{2}e^{i2kd}-\frac{\Gamma}{2}-i\frac{{\delta\Omega}}{2}$
(83)

The calculations show that equations (81) and (82) provide no dark states
$(\operatorname{Re}(\lambda_{1,2})=0)$ if $\delta\Omega$ is not equal to zero.
The quantities $\operatorname{Re}(\lambda_{1,2})$ are always negative at any
$kd$. The exception is the third root $\lambda_{3}$, with real part being
equal to zero for $kd=n\pi$. Therefore, if the root $\lambda_{3}$ does not
contribute to the dynamics of qubits, then at any $kd$, the qubits’ amplitudes
will decay to zero with the rate being dependent on $\delta\Omega$.

The dependence of real parts of $\lambda_{n}$ on $kd$ is shown in Fig.7 for
$\delta\Omega=\Gamma$.

Figure 7: Dependence of $\operatorname{Re}(\lambda_{n})/\Gamma$ on $kd$ for
$\delta\Omega=\Gamma$ . Solid (red) line relates to $\lambda_{1}$, dashed
(blue) line relates to $\lambda_{2}$, and dashed-dotted (green) line relates
to $\lambda_{3}$. The arrows show the points of deep subradiant states where
$\operatorname{Re}(\lambda_{1,2})/\Gamma\approx-0.001$

If the second (central) qubit is excited, the amplitudes of the first and the
third qubit are the same $\beta_{1}(t)=\beta_{3}(t)\equiv\beta(t)$ . Then
three equations (78), (79), and (80) reduce to two equations, for $\beta_{2}$
and $\beta$.

$\frac{{d\beta}}{{dt}}=-\frac{\Gamma}{2}\left({1+e^{2ikd}+i\frac{{\delta\Omega}}{\Gamma}}\right)\beta(t)-\beta_{2}(t)\frac{\Gamma}{2}\left({\frac{{\Omega_{0}}}{\Omega}}\right)^{1/2}e^{ik_{0}d}$
(84)

$\frac{{d\beta_{2}}}{{dt}}=-\frac{\Gamma}{2}\left({1-i\frac{{\delta\Omega}}{\Gamma}}\right)\beta_{2}(t)-\beta(t)\Gamma\left({\frac{\Omega}{{\Omega_{0}}}}\right)^{1/2}e^{ikd}$
(85)

The characteristic roots of these equations are equal to $\lambda_{1}$ and
$\lambda_{2}$ given above in equations (81) and (82).

The solution of equations (84) and (85). with account for initial conditions
for the amplitudes, $\beta_{2}(0)=1,\;\beta(0)=0$ and their time derivatives

$\left.{\frac{{d\beta}}{{dt}}}\right|_{t=0}=-\frac{\Gamma}{2}\left({\frac{{\Omega_{0}}}{\Omega}}\right)^{1/2}e^{ik_{0}d}$
(86)

$\left.{\frac{{d\beta_{2}}}{{dt}}}\right|_{t=0}=-\frac{\Gamma}{2}\left({1-i\frac{{\delta\Omega}}{\Gamma}}\right)$
(87)

is similar to (26)

$\begin{gathered}\beta_{2}(t)=b_{1}e^{\lambda_{1}t}+b_{2}e^{\lambda_{2}t}\hfill\\\
\beta(t)=a_{1}e^{\lambda_{1}t}+a_{2}e^{\lambda_{2}t}\hfill\\\ \end{gathered}$
(88)

with

$\begin{gathered}b_{1}=-\frac{{\frac{\Gamma}{2}\left({1-i\frac{{\delta\Omega}}{\Gamma}}\right)+\lambda_{2}}}{{\lambda_{1}-\lambda_{2}}}\hfill\\\
b_{2}=\frac{{\frac{\Gamma}{2}\left({1-i\frac{{\delta\Omega}}{\Gamma}}\right)+\lambda_{1}}}{{\lambda_{1}-\lambda_{2}}}\hfill\\\
\end{gathered}$ (89)

$\begin{gathered}a_{1}=-\frac{\Gamma}{2}\left({\frac{{\Omega_{0}}}{\Omega}}\right)^{1/2}\frac{{e^{ik_{0}d}}}{{\lambda_{1}-\lambda_{2}}}\hfill\\\
a_{2}=\frac{\Gamma}{2}\left({\frac{{\Omega_{0}}}{\Omega}}\right)^{1/2}\frac{{e^{ik_{0}d}}}{{\lambda_{1}-\lambda_{2}}}\hfill\\\
\end{gathered}$ (90)

Using the detuning $\delta\Omega$ as the external parameter we can control the
decay rates of the qubits amplitudes. This is shown in Fig.8 for several
values of the detuning $\delta\Omega$ with the second qubit being initially
excited. In the case shown in Fig.8, the decay rates are governed by real
parts of two roots, $\lambda_{1}$ and $\lambda_{2}$ with
$|Re\lambda_{1}|\ll|Re\lambda_{2}|$. As $\delta\Omega/\Gamma$ increases from
$0$ to $1$, $|Re\lambda_{1}|$ decreases from $1.5$ to $1.346$ while
$|Re\lambda_{2}|$ increases from $0$ to $0.154$. Even though the first root
tends to slow down the decay rate as $\delta\Omega/\Gamma$ increases, its
influence for $t\gg\Gamma$ becomes negligible, so that the main contribution
to the decay rate for large times comes from the second root which speeds up
the decay rate as $\delta\Omega/\Gamma$ increases.

Therefore, the greater is $\delta\Omega$, the more is its influence on the
qubit decay rates.

Figure 8: Dependence of the decay rates on the detuning $\delta\Omega$ for
initially excited second qubit. $kd=2\pi$ . (a)-the decay rates of the
first(third) qubit; (b) the decay rates of the second qubit. The numbers in
the plots correspond to different detunings:
$1-\delta\Omega=0.2\Gamma;\,2-\delta\Omega=0.4\Gamma;\,3-\delta\Omega=0.6\Gamma;\,4-\delta\Omega=0.8\Gamma;\,5-\delta\Omega=\Gamma$.

If the first qubit is initially excited the solution is similar to (33):

$\begin{gathered}\beta_{1}(t)=\frac{b_{2}}{2}e^{\lambda_{1}t}+\frac{b_{1}}{2}e^{\lambda_{2}t}+\frac{1}{2}e^{{}_{\lambda_{3}t}}\hfill\\\
\beta_{2}(t)=a_{1}\left(e^{\lambda_{1}t}-e^{\lambda_{2}t}\right)\hfill\\\
\beta_{3}(t)=\frac{b_{2}}{2}e^{\lambda_{1}t}+\frac{b_{1}}{2}e^{\lambda_{2}t}-\frac{1}{2}e^{{}_{\lambda_{3}t}}\end{gathered}$
(91)

where $\lambda_{1}$, $\lambda_{2}$, $\lambda_{3}$, and $a_{1}$, $b_{1}$,
$b_{2}$ are given in (81), (82), (83), and (89), (90).

For this case, we show in Fig. 9 the temporal behavior of qubits amplitudes
for several values of the detuning $\delta\Omega$. In the absence of detuning
(the panel _a_) we obtain the result shown in Fig.3. However, for non zero
detuning the amplitudes of the first and the third qubits decay to 0.25, while
the second qubit decays to zero. This difference can be explained by the
influence of dark states (the last terms in the expressions for $\beta_{(}t)$
and $\beta_{3}(t)$ in (91)). In general, as the detuning increases the decay
rates increase as well.

Figure 9: Dependence of the decay rates on the detuning $\delta\Omega$ for
initially excited the first qubit. $kd=2\pi$, $\Gamma/\Omega=0.001$. Every
panel from _a_ to _f_ shows the decay rates of the the first (black, solid
line), second (red, solid line), and third (blue, dashed line) qubit,
respectively. From top (_a_) to bottom (_f_) the quantity
$\delta\Omega/\Gamma$ is varied from $0$ (_a_) to $1$ (_f_) with the increment
being equal to $0.2$

The subsystem consisting of a central qubit and the symmetric qubit array can
be described as an analogue to a cavity QED system. In this description, the
central qubit acts as a two-level atom and the symmetric qubit array mimics a
high-finesse cavity, with the qubits array acting as an atomic mirror
Chang2012 . It seems that this analogy is supported by Fig.3 where two dark
states (70) of the _whole_ system prevent the qubit amplitudes from the
damping. If this analogy worked we would expect that the detuning of a central
qubit would not lead to any damping at all. In real cavity the qubit which is
large detuned from the cavity resonance cannot exchange its energy with the
cavity via real photons. However, as is seen from Fig.8, a full analogy
between resonance cavity and the qubit array does not exist. Even a small
detuning of the central qubit from the frequency of the qubit array results in
the disruption of the dark states that, in turn, leads to the damping of the
qubits amplitudes.

###  VII.1 Spectral density of spontaneous photon emission

For the calculation of the spectral density of photon emission we use the
equation (10) together with the substitution of $\beta_{n}(t)$ from equations
(15).

$\begin{gathered}\gamma_{k}(t)=-ig_{k}e^{ikd}\int\limits_{0}^{t}{}\bar{\beta}_{1}(t^{\prime})e^{i(\Omega-\Omega_{0})t^{\prime}/2}e^{i(\omega-\Omega)t^{\prime}}dt^{\prime}\hfill\\\
-ig_{k}\int\limits_{0}^{t}{}\bar{\beta}_{2}(t^{\prime})e^{-i(\Omega-\Omega_{0})t^{\prime}/2}e^{i(\omega-\Omega_{0})t^{\prime}}dt^{\prime}\hfill\\\
-ig_{k}e^{-ikd}\int\limits_{0}^{t}{}\bar{\beta}_{3}(t^{\prime})e^{i(\Omega-\Omega_{0})t^{\prime}/2}e^{i(\omega-\Omega)t^{\prime}}dt^{\prime}\hfill\\\
\end{gathered}$ (92)

where $kd=\frac{\omega}{{v_{g}}}d$,
$g_{k}=\left({\frac{{v_{g}\Gamma}}{{2L}}}\right)^{1/2}$.

The amplitudes $\bar{\beta}_{n}(t)$ in (92) are just the amplitudes
$\beta_{n}(t)$ given above in (88) and (91).

Below we show two plots which demonstrate how the spectral density of photon
emission depends on the detuning $\delta\Omega$. The dependencies of spectral
density on detuning are shown in Fig.10 for $kd=2\pi$ and in Fig.11 for
$kd=3\pi/2$, respectively. In both cases the second qubit was initially
excited.

If $\delta\Omega=0$ and $kd=2\pi$ we see a single Lorentzian peak (Fig.10(a))
with the width equal to $3\Gamma$ (see also (45). If $\delta\Omega\neq 0$ a
Lorentzian peak splits into two narrow peaks. The distance between the peaks,
$\Delta\omega$ is approximately equal to $(|Im\lambda_{1}|+|Im\lambda_{2}|)$.
The width of the peaks is determined by the root with a lowest real part. In
the case shown in Fig.10(b-f) the width of the peaks and their height are
approximately equal to $|Re\lambda_{2}|$ and $\Omega/|Re\lambda_{2}|$,
respectively. To estimate the order of these quantities we have found
numerically the range of variation for $Re\lambda_{2}$, $Im\lambda_{1}$, and
$Im\lambda_{2}$. As $\delta\Omega/\Gamma$ increases from $0.2$ to $1$, the
quantities $Re\lambda_{2}$, $Im\lambda_{1}$, and $Im\lambda_{2}$ vary from
$-0.006\Gamma$ to $-0.154\Gamma$, from $+0.028\Gamma$ to $+0.202\Gamma$, and
from $-0.032\Gamma$ to $-0.205\Gamma$, respectively.

Figure 10: Dependence of spectral density on the detuning. Second qubit is
initially excited. $kd=2\pi,\Gamma/\Omega=0.001$, $t\gg 1/\Gamma$. From top
(a) to bottom (f) the quantity $\delta\Omega/\Gamma$ is varied from $0$ (a) to
$1$ (f) with the increment being equal to $0.2$.

If $\delta\Omega=0$ and $kd=3\pi/2$ there are two similar peaks shown in
Fig.11(a) (see also Fig.6(a)). The distance between the peaks and their widths
are determined by the roots $\lambda_{1}$ and $\lambda_{2}$:
$Im\lambda_{1}$=$-Im\lambda_{2}$=$(\sqrt{7}/4)\Gamma$,
$Re\lambda_{1}$=$Re\lambda_{2}$=$-0.25\Gamma$ (see Eq. (21)). As
$\delta\Omega/\Gamma$ increases from $0.2$ to $1$, the quantities
$Re\lambda_{1}$ and $Re\lambda_{2}$ are of the same order of magnitude:
$Re\lambda_{1}$ varies from $-0.287\Gamma$ to $-0.399\Gamma$, and
$Re\lambda_{2}$ varies from $-0.213\Gamma$ to $-0.100\Gamma$. This behavior is
confirmed in Fig.11 (b-f) where the width of the left peak increases, while
the width of the right peak decreases. The distance between the peaks are
determined by the quantities $Im\lambda_{1}$ and $Im\lambda_{2}$ which vary
from $-0.708\Gamma$ to $-0.890\Gamma$, and from $+0.638\Gamma$ to
$+0.821\Gamma$, respectively.

Figure 11: Dependence of spectral density on the detuning. Second qubit is
initially excited. $kd=3\pi/2,\Gamma/\Omega=0.001$, $t\gg 1/\Gamma$. From top
(a) to bottom (f) the quantity $\delta\Omega/\Gamma$ is varied from $0$ (a) to
$1$ (f) with the increment being equal to $0.2$

We also calculated the temporal behavior of the full probability of the photon
emission , $P_{ph}(t)$ as a function of the detuning. This dependence is shown
in Fig.12. It is evident from the definition (8) that $P_{ph}(0)=0$ and as the
time increases it approaches either $1$ if all qubits amplitudes damp out to
zero, or a constant value if not all qubits decay to zero. Because the
detuning influences the rate of the qubits damping, the output rate of the
photon emission, $dP_{ph}/dt$ also depends on the detuning. The greater is the
detuning, the greater is the rate of the photon emission.

Figure 12: The probability amplitude of the photon emission.
$kd=2\pi,\Gamma/\Omega=0.001$. (a) the second qubit is initially excited; (b)
the first qubit is initially excited. The numbers in the plots correspond to
different detunings: $1-\delta\Omega=0$,
$2-\delta\Omega=0.2\Gamma;\,3-\delta\Omega=0.4\Gamma;\,4-\delta\Omega=0.6\Gamma;\,5-\delta\Omega=0.8\Gamma;\,6-\delta\Omega=\Gamma$.

##  VIII Conclusion

In this paper we have thoroughly investigated the dynamics behavior of qubits
amplitudes in 1D chain consisting of three qubits embedded in an open
waveguide. Within a single-excitation subspace we have found the evolution of
qubits amplitudes if one of the qubits was initially excited. We have shown
that even though the dark states do not contribute to the output photon
emission, they influence the evolution of qubits amplitudes in that they
prevent the qubits amplitudes from decaying to zero. We have found the
collective eigenstates of a three-qubit system and have shown how the qubits
amplitudes can be expressed in terms of the amplitudes of the collective
states. We also calculated the spectral density of the output photon emission
from three-qubit system and studied its dependence on the $kd$ value. We
studied the case when the frequency of the second qubit was different from
that of the edge qubits. In this case, the rates of qubits decay crucially
depend on the frequency detuning between central and edge qubits. The greater
is the detuning, the greater is the rate of the qubits damping.

We hope that this research will prove useful for the development of the
efficient control and readout protocols for a few-qubit quantum processor.

###### Acknowledgements.

Ya. S. G. thanks A. Sultanov for fruitful discussions. The work is supported
by the Ministry of Education and Science of Russian Federation under the
project FSUN-2020-0004.

##  Appendix A Equations for qubits amplitudes in the Wigner-Weisskopf
approximation

We assume that the first and the third qubits are identical
$(\Omega_{1}=\Omega_{3}\equiv\Omega,\;g_{k}^{(1)}=g_{k}^{(3)}\equiv g_{k})$
The frequency and the coupling of the second qubit are different
($\Omega_{2}\equiv\Omega_{0},\,g_{k}^{(2)}\equiv g_{k}^{(0)}$). A distance
between central qubit and the edge qubits is equal to $d$. We take the origin
in the location of the second qubit: $x_{1}=-d,x_{2}=0,x_{3}=+d$. For this
case, we expand the equation (9) as a set of three equations.

$\begin{gathered}\frac{{d\beta_{1}}}{{dt}}=-\sum\limits_{k}{g_{k}^{2}}\int\limits_{0}^{t}{\beta_{1}(t^{\prime})e^{-i(\omega_{k}-\Omega)(t-t^{\prime})}dt^{\prime}}\hfill\\\
-\sum\limits_{k}{g_{k}g_{k}^{(0)}e^{-ikd}e^{i(\Omega-\Omega_{0})t}}\int\limits_{0}^{t}{\beta_{2}(t^{\prime})e^{-i(\omega_{k}-\Omega_{0})(t-t^{\prime})}dt^{\prime}}\hfill\\\
-\sum\limits_{k}{g_{k}^{2}e^{-i2kd}}\int\limits_{0}^{t}{\beta_{3}(t^{\prime})e^{-i(\omega_{k}-\Omega)(t-t^{\prime})}dt^{\prime}}\hfill\\\
\end{gathered}$ (93)

$\begin{gathered}\frac{{d\beta_{2}}}{{dt}}=-\sum\limits_{k}{\left({g_{k}^{(0)}}\right)^{2}}\int\limits_{0}^{t}{\beta_{2}(t^{\prime})e^{-i(\omega_{k}-\Omega_{0})(t-t^{\prime})}dt^{\prime}}\hfill\\\
-\sum\limits_{k}{g_{k}^{(0)}g_{k}e^{ikd}e^{i(\Omega_{0}-\Omega)t}}\int\limits_{0}^{t}{\beta_{1}(t^{\prime})e^{-i(\omega_{k}-\Omega)(t-t^{\prime})}dt^{\prime}}\hfill\\\
-\sum\limits_{k}{g_{k}^{(0)}g_{k}e^{-ikd}e^{i(\Omega_{0}-\Omega)t}}\int\limits_{0}^{t}{\beta_{3}(t^{\prime})e^{-i(\omega_{k}-\Omega)(t-t^{\prime})}dt^{\prime}}\hfill\\\
\end{gathered}$ (94)

$\begin{gathered}\frac{{d\beta_{3}}}{{dt}}=-\sum\limits_{k}{g_{k}^{2}}\int\limits_{0}^{t}{\beta_{3}(t^{\prime})e^{-i(\omega_{k}-\Omega)(t-t^{\prime})}dt^{\prime}}\hfill\\\
-\sum\limits_{k}{g_{k}^{2}e^{i2kd}}\int\limits_{0}^{t}{\beta_{1}(t^{\prime})e^{-i(\omega_{k}-\Omega)(t-t^{\prime})}dt^{\prime}}\hfill\\\
-\sum\limits_{k}{g_{k}g_{k}^{(0)}e^{ikd}e^{i(\Omega-\Omega_{0})t}}\int\limits_{0}^{t}{\beta_{2}(t^{\prime})e^{-i(\omega_{k}-\Omega_{0})(t-t^{\prime})}dt^{\prime}}\hfill\\\
\end{gathered}$ (95)

According to Wigner-Weisskopf approach we replace
$\beta_{1}(t^{\prime}),\beta_{2}(t^{\prime}),\beta_{3}(t^{\prime})$ in the
integrands with $\beta_{1}(t),\beta_{2}(t),\beta_{3}(t)$ and take them out of
the integrals.

$\begin{gathered}\frac{{d\beta_{1}}}{{dt}}=-\beta_{1}(t)\sum\limits_{k}{g_{k}^{2}}I_{k}(\Omega,t)\hfill\\\
-\beta_{2}(t)\sum\limits_{k}{g_{k}g_{k}^{(0)}e^{-ikd}e^{i(\Omega-\Omega_{0})t}}I_{k}(\Omega_{0},t)\hfill\\\
-\beta_{3}(t)\sum\limits_{k}{g_{k}^{2}e^{-i2kd}}I_{k}(\Omega,t)\hfill\\\
\end{gathered}$ (96)

$\begin{gathered}\frac{{d\beta_{2}}}{{dt}}=-\beta_{2}(t)\sum\limits_{k}{|g_{k}^{(0)}|^{2}}I_{k}(\Omega_{0},t)\hfill\\\
-\beta_{1}(t)\sum\limits_{k}{g_{k}^{(0)}g_{k}e^{ikd}e^{-i(\Omega-\Omega_{0})t}}I_{k}(\Omega,t)\hfill\\\
-\beta_{3}(t)\sum\limits_{k}{g_{k}^{(0)}g_{k}e^{-ikd}e^{-i(\Omega-\Omega_{0})t}}I_{k}(\Omega,t)\hfill\\\
\end{gathered}$ (97)

$\begin{gathered}\frac{{d\beta_{3}}}{{dt}}=-\beta_{3}(t)\sum\limits_{k}{g_{k}^{2}}I_{k}(\Omega,t)\hfill\\\
-\beta_{1}(t)\sum\limits_{k}{g_{k}^{2}e^{i2kd}}I_{k}(\Omega,t)\hfill\\\
-\beta_{2}(t)\sum\limits_{k}{g_{k}g_{k}^{(0)}e^{ikd}e^{i(\Omega-\Omega_{0})t}}I_{k}(\Omega_{0},t)\hfill\\\
\end{gathered}$ (98)

where

$\begin{gathered}I_{k}(\Omega,t)=\int\limits_{0}^{t}{e^{-i(\omega_{k}-\Omega)(t-t^{\prime})}dt^{\prime}}=\int\limits_{0}^{t}{e^{-i(\omega_{k}-\Omega)\tau}d\tau}\hfill\\\
\approx\int\limits_{0}^{\infty}{e^{-i(\omega_{k}-\Omega)\tau}d\tau}=\pi\delta(\omega_{k}-\Omega)-iP.v.\left({\frac{1}{{\omega_{k}-\Omega}}}\right)\hfill\\\
\end{gathered}$ (99)

where $P.v.$ is a Cauchy principal value integral.

We can remove oscillating exponents in (96)-(98) with the aid of the
substitution

$\begin{gathered}\beta_{1}(t)=e^{i(\Omega-\Omega_{0})t/2}\bar{\beta}_{1}(t)\hfill\\\
\beta_{2}(t)=e^{-i(\Omega-\Omega_{0})t/2}\bar{\beta}_{2}(t)\hfill\\\
\beta_{3}(t)=e^{i(\Omega-\Omega_{0})t/2}\bar{\beta}_{3}(t)\hfill\\\
\end{gathered}$ (100)

In addition, we assume the coupling constants $g_{k}$ are even functions of
$k$ ($g_{k}=g_{-k}$). Then from (96)-(98) we obtain

$\begin{gathered}\frac{{d\bar{\beta}_{1}}}{{dt}}=-\bar{\beta}_{1}(t)\sum\limits_{k}{g_{k}^{2}}I_{k}(\Omega,t)-i\frac{{\Omega-\Omega_{0}}}{2}\bar{\beta}_{1}(t)\hfill\\\
-\bar{\beta}_{2}(t)2\sum\limits_{k>0}{g_{k}g_{k}^{(0)}\cos(kd)}I_{k}(\Omega_{0},t)\hfill\\\
-\bar{\beta}_{3}(t)2\sum\limits_{k>0}{g_{k}^{2}\cos(2kd)}I_{k}(\Omega,t)\hfill\\\
\end{gathered}$ (101)

$\begin{gathered}\frac{{d\bar{\beta}_{2}}}{{dt}}=-\bar{\beta}_{2}(t)\sum\limits_{k}{\left({g_{k}^{(0)}}\right)^{2}}I_{k}(\Omega_{0},t)+i\frac{{\Omega-\Omega_{0}}}{2}\bar{\beta}_{2}(t)\hfill\\\
-\bar{\beta}_{1}(t)2\sum\limits_{k>0}{g_{k}^{(0)}g_{k}\cos(kd)}I_{k}(\Omega,t)\hfill\\\
-\bar{\beta}_{3}(t)2\sum\limits_{k>0}{g_{k}^{(0)}g_{k}\cos(kd)}I_{k}(\Omega,t)\hfill\\\
\end{gathered}$ (102)

$\begin{gathered}\frac{{d\bar{\beta}_{3}}}{{dt}}=-\bar{\beta}_{3}(t)\sum\limits_{k}{g_{k}^{2}}I_{k}(\Omega,t)-i\frac{{\Omega-\Omega_{0}}}{2}\bar{\beta}_{3}(t)\hfill\\\
-\bar{\beta}_{1}(t)2\sum\limits_{k>0}{g_{k}^{2}\cos(2kd)}I_{k}(\Omega,t)\hfill\\\
-\bar{\beta}_{2}(t)2\sum\limits_{k>0}{g_{k}g_{k}^{(0)}\cos(kd)}I_{k}(\Omega_{0},t)\hfill\\\
\end{gathered}$ (103)

The next step is to relate the coupling constants $g_{k}$ to the qubit decay
rate of spontaneous emission into waveguide mode. In accordance with Fermi
golden rule we define the qubit decay rates by the following expressions:

$\Gamma=2\pi\sum\limits_{k}{g_{k}^{2}\delta(\omega_{k}-\Omega)}$ (104)

$\Gamma_{0}=2\pi\sum\limits_{k}{\left({g_{k}^{(0)}}\right)^{2}\delta(\omega_{k}-\Omega_{0})}$
(105)

where

$g_{k}=\sqrt{\frac{{\omega_{k}D^{2}}}{{2\hbar\varepsilon_{0}V}}},$ (106)

$D$ is the matrix element of the qubit’s dipole moment operator, $V$ is the
effective volume where the interaction between qubit and electromagnetic field
takes place. For 1D case, a summation over $k$ is replaced by the integration
over $\omega$ in accordance with the prescription:

$\sum\limits_{k}{}\Rightarrow\frac{L}{{2\pi}}\int\limits_{-\infty}^{\infty}{dk}=\frac{L}{{2\pi}}2\int\limits_{0}^{\infty}{d\left|k\right|}=\frac{L}{{\pi\upsilon_{g}}}\int\limits_{0}^{\infty}{d\omega_{k}}$
(107)

where $L$ is a length of the waveguide, and we assumed a linear dispersion
law, $\omega_{k}=v_{g}k$. The application of (107) to, for example, (104),
allows the relation between a coupling constant $g_{k}$ and the decay rate
$\Gamma$.

$g_{\Omega}=\sqrt{\frac{{\Omega
D^{2}}}{{2\hbar\varepsilon_{0}V}}}=\left({\frac{{v_{g}\Gamma}}{{2L}}}\right)^{1/2}$
(108)

Therefore, we may relate the coupling constants
$g_{\Omega}^{(0)},g_{\Omega_{2}}^{(0)},g_{\Omega_{0}}$ with their respective
decay rates.

$g_{\Omega}^{(0)}=\sqrt{\frac{{\Omega
D_{0}^{2}}}{{2\hbar\varepsilon_{0}V}}}=\left({\frac{\Omega}{{\Omega_{0}}}}\right)^{1/2}\left({\frac{{v_{g}}}{{2L}}\Gamma_{0}}\right)^{1/2}$
(109)

$g_{\Omega_{0}}=\sqrt{\frac{{\Omega_{0}D^{2}}}{{2\hbar\varepsilon_{0}V}}}=\left({\frac{{\Omega_{0}}}{\Omega}}\right)^{1/2}\left({\frac{{v_{g}}}{{2L}}\Gamma}\right)^{1/2}$
(110)

$g_{\Omega_{0}}^{(0)}=\sqrt{\frac{{\Omega_{0}D_{0}^{2}}}{{2\hbar\varepsilon_{0}V}}}=\left({\frac{{v_{g}}}{{2L}}\Gamma_{0}}\right)^{1/2}$
(111)

Now we can calculate the different terms in (101)-(103). We begin with the sum
in the first line in (101).

$\begin{gathered}\sum\limits_{k}{g_{k}^{2}}I_{k}(\Omega,t)=\sum\limits_{k}{g_{k}^{2}}\left({\pi\delta(\omega_{k}-\Omega)-iP.v.\left({\frac{1}{{\omega_{k}-\Omega}}}\right)}\right)\hfill\\\
=\frac{\Gamma}{2}-iP.v.\sum\limits_{k}\left({\frac{{g_{k}^{2}}}{{\omega_{k}-\Omega}}}\right)\approx\frac{\Gamma}{2}\hfill\\\
\end{gathered}$ (112)

where $\Gamma$ is given in (104).

The second term in (112) gives rise to the shift of the qubit frequency.
Therefore, we incorporate it in the renormalized qubit frequency and will not
write it explicitly any more. The sum in the second line in (101) is
calculated as follows:

$\begin{gathered}2\sum\limits_{k>0}{g_{k}g_{k}^{(0)}\cos(kd)}I_{k}(\Omega_{0},t)=\hfill\\\
=\frac{L}{{\upsilon_{g}}}\int\limits_{0}^{\infty}{g_{k}^{(0)}g_{k}\cos(kd)\delta(\omega_{k}-\Omega_{2})d\omega_{k}}\hfill\\\
-2iP.v.\sum\limits_{k>0}\left({\frac{{g_{k}g_{k}^{(0)}\cos(kd)}}{{\omega_{k}-\Omega_{{}_{0}}}}}\right)\hfill\\\
=\frac{L}{{\upsilon_{g}}}g_{{}_{\Omega
0}}^{(0)}g_{\Omega_{0}}\cos(k_{0}d)-i\frac{L}{{v_{g}\pi}}g_{\Omega_{0}}g_{\Omega_{0}}^{(0)}P.v.\int\limits_{0}^{\infty}{d\omega}\frac{{\cos\left({\frac{\omega}{{v_{g}}}d}\right)}}{{\omega-\Omega_{0}}}\hfill\\\
\end{gathered}$ (113)

For principle value integral in (113) we obtain with a good accuracy (see
Appendix B):

$P.v.\int\limits_{0}^{\infty}{d\omega}\frac{{\cos\left({\frac{\omega}{{v_{g}}}d}\right)}}{{\omega-\Omega}}\approx-\pi\sin\left({\frac{\Omega}{{v_{g}}}d}\right)=-\pi\sin\left({k_{\Omega}d}\right)$
(114)

Therefore, we finally obtain:

$\begin{gathered}2\sum\limits_{k>0}{g_{k}g_{k}^{(0)}\cos(kd)}I_{k}(\Omega_{0},t)=\frac{L}{{\upsilon_{g}}}g_{{}_{\Omega_{0}}}^{(0)}g_{\Omega_{0}}e^{ik_{0}d}\hfill\\\
=\frac{1}{2}\left({\frac{{\Omega_{0}}}{\Omega}}\right)^{1/2}\sqrt{\Gamma\Gamma_{0}}e^{ik_{0}d}\hfill\\\
\end{gathered}$ (115)

Similar calculations give for the last sum in (101):

$2\sum\limits_{k>0}{g_{k}^{2}\cos(2kd)}I_{k}(\Omega,t)=\frac{\Gamma}{2}e^{2ikd}$
(116)

Collecting together (112), (115), and (116) we write the final form of
equation (101):

$\begin{gathered}\frac{{d\bar{\beta}_{2}}}{{dt}}=-\frac{{\Gamma_{0}}}{2}\bar{\beta}_{2}(t)+i\frac{{\Omega-\Omega_{0}}}{2}\bar{\beta}_{2}(t)\hfill\\\
-\frac{1}{2}\left({\frac{\Omega}{{\Omega_{0}}}}\right)^{1/2}\sqrt{\Gamma\Gamma_{0}}e^{ikd}\left({\bar{\beta}_{1}(t)+\bar{\beta}_{3}(t)}\right)\hfill\\\
\end{gathered}$ (117)

Similar calculations for the equations (102) and (103) yield the following
result:

$\begin{gathered}\frac{{d\bar{\beta}_{2}}}{{dt}}=-\frac{{\Gamma_{0}}}{2}\bar{\beta}_{2}(t)+i\frac{{\Omega-\Omega_{0}}}{2}\bar{\beta}_{2}(t)\hfill\\\
-\frac{1}{2}\left({\frac{\Omega}{{\Omega_{0}}}}\right)^{1/2}\sqrt{\Gamma\Gamma_{0}}e^{ikd}\left({\bar{\beta}_{1}(t)+\bar{\beta}_{3}(t)}\right)\hfill\\\
\end{gathered}$ (118)

$\begin{gathered}\frac{{d\bar{\beta}_{3}}}{{dt}}=-\frac{\Gamma}{2}\bar{\beta}_{3}(t)-i\frac{{\Omega-\Omega_{0}}}{2}\bar{\beta}_{3}(t)-\bar{\beta}_{1}(t)\frac{\Gamma}{2}e^{2ikd}\hfill\\\
-\bar{\beta}_{2}(t)\frac{1}{2}\left({\frac{{\Omega_{0}}}{\Omega}}\right)^{1/2}\sqrt{\Gamma\Gamma_{0}}e^{ik_{0}d}\hfill\\\
\end{gathered}$ (119)

##  Appendix B Proof of equation (114)

The integral

$\int\limits_{0}^{\infty}{d\omega}\frac{{\cos\left({\frac{\omega}{{v_{g}}}d}\right)}}{{\omega-\Omega}}=\int\limits_{0}^{\infty}{dx}\frac{{\cos\left({ax}\right)}}{{x-1}}$

, where $a=k_{\Omega}d$, can be expressed in terms of sine and cosine
integrals, $ci$ and $si$ Prud :

$\begin{gathered}\int\limits_{0}^{\infty}{dx}\frac{{\cos\left({ax}\right)}}{{x-1}}=-\cos(a)ci(a)-\sin(a)[si(a)+\pi]\hfill\\\
=-\cos(a)Ci(a)-\sin(a)[Si(a)+\frac{\pi}{2}]\hfill\\\ \end{gathered}$ (120)

where

$Ci(a)=-\int\limits_{a}^{\infty}{}\frac{{\cos
t}}{t}dt;\;Si(a)=\int\limits_{0}^{a}{}\frac{{\sin t}}{t}dt$ (121)

The integrand in left hand side in (120) has a singular point at $x=1$ which
manifests itself as a singularity of $Ci(a)$ at $a\rightarrow 0$ in right hand
side in (120). Therefore, we calculate integral (120) as Cauchy principal
value integral.

$\begin{gathered}P.v.\int\limits_{0}^{\infty}{\frac{{\cos
ax}}{{x-1}}}dx\approx P.v.\int\limits_{-\infty}^{\infty}{\frac{{\cos
ax}}{{x-1}}}dx=P.v.\int\limits_{-\infty}^{\infty}{\frac{{\cos
a(t+1)}}{t}}dt\hfill\\\ =\cos
a\;P.v.\int\limits_{-\infty}^{\infty}{\frac{{\cos at}}{t}}dt-\sin
a\;P.v.\int\limits_{-\infty}^{\infty}{\frac{{\sin at}}{t}}dt\hfill\\\
\hfill\\\ \end{gathered}$ (122)

$P.v.\int\limits_{-\infty}^{\infty}{\frac{{\cos
at}}{t}}dt=\mathop{\lim}\limits_{R\to\infty}\mathop{\lim}\limits_{\varepsilon\to
0}\left({\int\limits_{-R}^{-\varepsilon}{}\frac{{\cos
at}}{t}dt+\int\limits_{\varepsilon}^{R}{}\frac{{\cos at}}{t}dt}\right)$ (123)

$P.v.\int\limits_{-\infty}^{\infty}{\frac{{\sin
at}}{t}}dt=\mathop{\lim}\limits_{R\to\infty}\mathop{\lim}\limits_{\varepsilon\to
0}\left({\int\limits_{-R}^{-\varepsilon}{}\frac{{\sin
at}}{t}dt+\int\limits_{\varepsilon}^{R}{}\frac{{\sin at}}{t}dt}\right)$ (124)

Since

$\int\limits_{-R}^{-\varepsilon}{}\frac{{\cos
at}}{t}dt=-\int\limits_{\varepsilon}^{R}{}\frac{{\cos at}}{t}dt$

and

$\int\limits_{-R}^{-\varepsilon}{}\frac{{\sin
at}}{t}dt=\int\limits_{\varepsilon}^{R}{}\frac{{\sin at}}{t}dt$

we obtain

$P.v.\int\limits_{-\infty}^{\infty}{\frac{{\cos at}}{t}}dt=0$ (125)

$\begin{gathered}P.v.\int\limits_{-\infty}^{\infty}{\frac{{\sin
at}}{t}}dt=\mathop{\lim}\limits_{R\to\infty}\mathop{\lim}\limits_{\varepsilon\to
0}\left({2\int\limits_{\varepsilon}^{R}{}\frac{{\sin
at}}{t}dt}\right)=\hfill\\\
\mathop{\lim}\limits_{R\to\infty}\left({2\int\limits_{0}^{R}{}\frac{{\sin
at}}{t}dt}\right)=\mathop{\lim}\limits_{R\to\infty}\left({2Si(aR)}\right)=\pi\hfill\\\
\end{gathered}$ (126)

The final result in (126) is due to the known relation Prud :
$\mathop{\lim}\limits_{x\to\infty}Si(x)=\frac{\pi}{2}$.

Therefore, we finally obtain

$P.v.\int\limits_{0}^{\infty}{\frac{{\cos ax}}{{x-1}}}dx\approx-\pi\sin a$
(127)

Below, in Fig.13 we compare the $kd$\- dependence of (120) with that of (127).
We see a noticeable discrepancy for $kd<\pi/4$. For $kd>\pi/2$ two curves are
almost identical.

Figure 13: Comparison of principle value $F_{1}$ (right hand side of (127)
with the ”exact” expression $F_{2}$ (right hand side of (120)).

## References

  * (1) O. Astafiev, A. M. Zagoskin, A. A. Abdumalikov, Yu. A. Pashkin, T. Yamamoto, K. Inomata, Y. Nakamura, and J. S. Tsai, Resonance Fluorescence of a Single Artificial Atom. Science 327, 840 (2010).
  * (2) I.-C. Hoi, C. M. Wilson, G. Johansson, T. Palomaki, B. Peropadre, and P. Delsing, Phys. Demonstration of a Single-Photon Router in the Microwave Regime. Rev. Lett. 107, 073601 (2011).
  * (3) A. F. van Loo, A. Fedorov, K. Lalumi‘ere, B. C. Sanders, A. Blais, and A. Wallra, Photon-Mediated Interactions Between Distant Artificial Atoms. Science 342, 1494 (2013).
  * (4) D. Roy, C. M. Wilson, and O. Firstenberg, Strongly interacting photons in one-dimensional continuum Rev. Mod. Phys. 89, 021001 (2017).
  * (5) X. Gu, A. F. Kockum, A. Miranowicz, Yu- xi Liu, and F. Nori, Microwave photonics with superconducting quantum circuits. Phys. Rep. 718-719, 1 (2017).
  * (6) S. N. Shevchenko, Mesoscopic Physics Meets Quantum Engineering. World Scientific, Singapore (2019).
  * (7) A. Albrecht, L. Henriet, A. Asenjo-Garcia, P. B. Dieterle, O. Painter, and D. E. Chang, Subradiant states of quantum bits coupled to a one-dimensional waveguide, New J. Phys. 21, 025003 (2019).
  * (8) Y.-X. Zhang and K. Molmer, Theory of Subradiant States of a One-Dimensional Two-Level Atom Chain, Phys. Rev. Lett. 122, 203605 (2019).
  * (9) J. Ruostekoski and J. Javanainen, Arrays of strongly coupled atoms in a one-dimensional waveguide. Phys. Rev. A 96, 033857 (2017).
  * (10) K. Lalumi‘ere, B. C. Sanders, A. F. van Loo, A. Fedorov, A. Wallraff, and A. Blais, Input-output theory for waveguide QED with an ensemble of inhomogeneous atoms. Phys. Rev. A 88, 043806 (2013).
  * (11) D. E. Chang, L. Jiang, A. V. Gorshkov, and H. J. Kimble, Cavity QED with atomic mirrors. New J. Phys. 14, 063003 (2012).
  * (12) M. Mirhosseini, E. Kim, X. Zhang, A. Sipahigil, P. B. Dieterle, A. J. Keller,A. Asenjo-Garcia, D. E. Chang, and O. Painter, Cavity quantum electrodynamics with atom-like mirrors. Nature 569, 692 (2019).
  * (13) J. D. Brehm, A. N. Poddubny, A. Stehli, T. Wolz, H. Rotzinger, and A. V. Ustinov, Waveguide bandgap engineering with an array of superconducting qubits, npj Quantum Materials 6, 10 (2021).
  * (14) A. F. van Loo, A. Fedorov, K. Lalumi‘ere, B. C. Sanders, A. Blais, and A. Wallraff Photon-mediated interactions between distant artificial atoms. Science 342, 1494 (2014).
  * (15) E. Barnes, C. Arenz, A. Pitchford, and S. E. Economou, Fast microwave-driven three-qubit gates for cavity-coupled superconducting qubits. Phys. Rev B 96, 024504 (2017).
  * (16) L. T. Kenfack, M. Tchoffo, G. C. Fouokeng, and L. C. Fai, Dynamical evolution of entanglement of a three-qubit system driven by a classical environmental colored noise Quant. Inf. Process 17, 76 (2018).
  * (17) T. S. Tsoi and C. K. Law, Quantum interference effects of a single photon interacting with an atomic chain inside a one-dimensional waveguide. Phys. Rev. A 78, 063832 (2008).
  * (18) Yu Zhou, Zh. Zhang, Z. Yin, S. Huai, X. Gu, X. Xu, J. Allcock, F. Liu, G.Xi, Q. Yu, H. Zhang, M. Zhang, H. Li, X. Song, Zh. Wang, D. Zheng, Sh. An, Y. Zheng, and Sh. Zhang, Rapid and Unconditional Parametric Reset Protocol for Tunable Superconducting Qubits. ArXive: 2103.11315 [quant-ph] (2021).
  * (19) P. Magnard,, Ph. Kurpiers, B. Royer,, T, Walter, J.-C. Besse,, S. Gasparinetti,, M. Pechal,, J. Heinsoo,, S. Storz,, A. Blais, and, A, Wallraff, Fast and Unconditional All-Microwave Reset of a Superconducting Qubit. Phys. Rev. Lett. 121, 060502 (2018).
  * (20) R. H. Lehmberg, Radiation from an N-Atom System. I. General Formalism. Phys. Rev. A 2, 883 (1970).
  * (21) F. Keck, H. J. Korsch, and S. Mossmann, Unfolding a diabolic point: a generalized crossing scenario. J. Phys. A: Math. Gen. 36, 2125 (2003).
  * (22) D. C. Brody, Biorthogonal quantum mechanics. J. Phys. A: Math. Theor. 47, 035305 (2014).
  * (23) A. P. Prudnikov, Yu. A. Brychkov, and O. I. Marichev, Integrals and Series. Volume 1: Elementary Functions. Gordon and Breach Science Publishers / CRC Press. New York-London. 1986, 798 pages.