# On the Dynamics of the Tavis-Cummings Model

Zhiyuan Dong, Guofeng Zhang, Ai-Guo Wu, and Re-Bing Wu Zhiyuan Dong is with
School of Science, Harbin Institute of Technology, Shenzhen, Shenzhen, China
(e-mail: dongzhiyuan@hit.edu.cn).Guofeng Zhang is with the Department of
Applied Mathematics, The Hong Kong Polytechnic University, Hong Kong. He is
also with The Hong Kong Polytechnic University Shenzhen Research Institute,
Shenzhen, 518057, China (e-mail: guofeng.zhang@polyu.edu.hk).Ai-Guo Wu is with
School of Mechanical Engineering and Automation, Harbin Institute of
Technology, Shenzhen, Shenzhen, China (e-mail<EMAIL_ADDRESS>ag.wu@163.com).Re-Bing Wu is with the Center for Intelligent and Networked
Systems, Department of Automation, Tsinghua University, Beijing, 100084, China
(e-mail: rbwu@tsinghua.edu.cn).

###### Abstract

The purpose of this paper is to present a comprehensive study of the Tavis-
Cummings model from a system-theoretic perspective. A typical form of the
Tavis-Cummings model is composed of an ensemble of non-interacting two-level
systems (TLSs) that are collectively coupled to a common cavity resonator. The
associated quantum linear passive system is proposed, whose canonical form
reveals typical features of the Tavis-Cummings model, including $\sqrt{N}$\-
scaling, dark states, bright states, single-excitation superradiant and
subradiant states. The passivity of this linear system is related to the
vacuum Rabi mode splitting phenomenon in Tavis-Cummings systems. On the basis
of the linear model, an analytic form is presented for the steady-state output
state of the Tavis-Cummings model driven by a single-photon state. Master
equations are used to study the excitation properties of the Tavis-Cummings
model in the multi-excitation scenario. Finally, in terms of the transition
matrix for a linear time-varying system, a computational framework is proposed
for calculating the state of the Tavis-Cummings model, which is applicable to
the multi-excitation case.

###### Index Terms:

Quantum control, Tavis-Cummings model, two-level systems, open quantum
systems.

##  I Introduction

In 1954, Robert Dicke calculated [1] that an ensemble of gaseous molecules
interacting with a common radiation field could exhibit a coherent spontaneous
emission process, during which the molecules act as a giant molecule that
shows superradiation — cooperative radiation rate much faster than independent
individual radiation rates. This problem was further studied by Tavis and
Cummings [2] by means of a model of $N$ identical non-interacting two-level
systems (TLSs) coupled to a single-mode quantized radiation field, see Fig. 1
in section III for an example. An exact solution of the eigenstates of the
Hamiltonian for this model is derived. The model proposed in [2] is called the
Tavis-Cummings model in the subsequent literature. The Tavis-Cummings model
has been physically realized by quite a few experimental platforms, including
superconducting circuits [3, 4, 5, 6], NV spin ensembles [7], and double
quantum dots [8, 9, 10]. In the Tavis-Cummings model consisting of $N$ TLSs
equally coupled to a cavity resonator, under certain conditions the collective
coupling strength of the ensemble exhibits a $\sqrt{N}$-scaling, which can be
experimentally observed from the vacuum Rabi mode splitting of the resonator
transmission spectrum of the bright states. In addition to bright states, a
Tavis-Cummings model can also have dark states which contain single-excitation
subradiant states as a subclass. Applications of the Tavis-Cummings model can
be found in [11, 4, 9, 6] and references therein.

In this paper, we aim to introduce the Tavis-Cummings model to the quantum
control community and show that many of its typical properties can be
uncovered by means of systems theory. The main contributions are summarized
below.

In section III-B we propose a quantum linear system that is associated to the
Tavis-Cummings model. This linear model reveals the $\sqrt{N}$-scaling of the
coupling strength of the atomic ensemble. Moreover, a transfer function is
defined for a performance variable for which the system is passive, and the
transfer function reflects the vacuum Rabi mode spitting of the Tavis-Cummings
model. Finally, the structural decomposition of the linear model shows that
the bright states of the Tavis-Cummings model live in the controllable and
observable subspace, whereas the dark states reside in the uncontrollable and
unobservable subspace of the linear model.

In section IV-A we apply the quantum linear systems theory to derive an
analytic form of the output single-photon state of the Tavis-Cummings system
driven by a single-photon input. The simulations in Fig. 2 show that the input
photon tends not to interact with the atoms when the number of atoms is large,
and photon-atom interaction is easier when atoms are non-resonant. On the
other hand, when only one of the two-level atoms is initially excited and the
input field is vacuum, an analytic form of the joint system-field state is
given in section IV-C, which explains several experimental observations
including superradiance and subradiance [6].

In section V, we study the excitations of TLSs by means of master equations.
When all the $N$ atoms are initially in the excited state, we prove that
eventually they all settle to the ground state and the output field is in an
$N$-photon state, provided that the coupling strengths are identical.

In our preliminary study [13], a computational framework is proposed to
calculate the joint system-field state of a general open quantum system. Here,
we develop it further in section VI and apply it to the study of the Tavis-
Cummings model. In particular, we derive the exact form of 2- and 3-photon
states.

Notation. The reduced Planck constant $\hbar$ is set to $1$.
$\mathrm{i}=\sqrt{-1}$ is the imaginary unit, $\delta_{ij}$ denotes the
Kronecker delta function, and $\delta(t-r)$ is the Dirac delta function. $I$
is the identity matrix, and $\mathbf{0}$ is the zero vector or matrix whose
dimension can be easily determined from the context. Given a column vector of
complex numbers or operators $X=[x_{1},\ldots,x_{n}]^{\top}$, the complex
conjugate or adjoint operator of $X$ is denoted by
$X^{\\#}=[x_{1}^{\ast},\ldots,x_{n}^{\ast}]^{\top}$. Let
$X^{\dagger}=(X^{\\#})^{\top}$. Clearly, when $n=1$, $X^{\dagger}=X^{\ast}$.
Let $t_{0}$ be the initial time, i.e., the time when the system and its input
start interaction. $|g\rangle$ and $|e\rangle$ stand for the ground and
excited states of a two-level atom, respectively. The commutator between two
operators $A$ and $B$ is $[A,B]=AB-BA$. Define super-operators
$\mathcal{L}X\triangleq-\mathrm{i}[X,H]+\mathcal{D}_{L}X$ and
$\mathcal{L}^{\star}\rho\triangleq-\mathrm{i}[H,\rho]+\mathcal{D}_{L}^{\star}\rho$,
where
$\mathcal{D}_{L}X=L^{\dagger}XL-\frac{1}{2}L^{\dagger}LX-\frac{1}{2}XL^{\dagger}L$
and $\mathcal{D}_{L}^{\star}\rho=L\rho
L^{\dagger}-\frac{1}{2}L^{\dagger}L\rho-\frac{1}{2}\rho L^{\dagger}L$.

##  II Preliminaries

###  II-A Quantum systems and fields

Consider a quantum system driven by $m$ input fields. The inputs are optical
or microwave fields, which are represented by the annihilation operators
$b_{\mathrm{in},k}(t)$ and their adjoints $b_{\mathrm{in},k}^{\ast}(t)$
(creation operators), $k=1,\ldots,m$. If there are no photons in an input
channel, this input is in the vacuum state $\ket{\Phi_{0}}$. Annihilation and
creation operators satisfy

$\displaystyle b_{\mathrm{in},j}(t)\ket{\Phi_{0}}=0,\
[b_{\mathrm{in},j}(t),b_{\mathrm{in},k}(r)]=[b_{\mathrm{in},j}^{\ast}(t),b_{\mathrm{in},k}^{\ast}(r)]=0,$
(1)
$\displaystyle[b_{\mathrm{in},j}(t),b_{\mathrm{in},k}^{\ast}(r)]=\delta_{jk}\delta(t-r),\forall
j,k=1,\ldots,m,\ t,r\in\mathbb{R}.$

The integrated input annihilation and creation processes are respectively
$B_{\mathrm{in},k}(t)=\int_{-\infty}^{t}b_{\mathrm{in},k}(r)dr$ and
$B_{\mathrm{in},k}^{\ast}(t)=\int_{-\infty}^{t}b_{\mathrm{in},k}^{\ast}(r)dr$,
which are quantum Wiener processes. Define Itô increments
$dB_{\mathrm{in},k}(t)=B_{\mathrm{in},k}(t+dt)-B_{\mathrm{in},k}(t)$. Denote
$B_{\mathrm{in}}(t)=[B_{\mathrm{in},1}(t),\ldots,B_{\mathrm{in},m}(t)]^{\top}$.
In this paper, the input fields are assumed to be canonical fields which
include the vacuum, coherent, single- and multi-photon fields. Then
$dB_{\mathrm{in}}(t)dB_{\mathrm{in}}^{\top}(\tau)=dB_{\mathrm{in}}^{\\#}(t)dB_{\mathrm{in}}^{\dagger}(\tau)=dB_{\mathrm{in}}^{\\#}(t)dB_{\mathrm{in}}^{\top}(\tau)=\mathbf{0}$
and $dB_{\mathrm{in}}(t)dB_{\mathrm{in}}^{\dagger}(\tau)=I\delta_{t\tau}dt$.

The quantum system can be parametrized by a triple $(S,L,H)$ [15, 16]. Here,
$H$ is the inherent system Hamiltonian, $L=[L_{1},\ldots,L_{m}]^{\top}$
describes how the system is coupled to its environment, and $S$ is a
scattering operator (e.g., a beamsplitter or a phase shifter). In this paper,
it is assumed that $S=I$ (the identity operator). The temporal evolution of
the quantum system is governed by a unitary operator $U(t,t_{0})$, which is
the solution to the following Itô quantum stochastic differential equation
(QSDE):

$\displaystyle
dU(t,t_{0})=\left[-(\mathrm{i}H+\frac{1}{2}L^{\dagger}L)dt\right.$
$\displaystyle\hskip
60.27759pt\left.+dB_{\mathrm{in}}^{\dagger}(t)L-L^{\dagger}dB_{\mathrm{in}}(t)\right]U(t,t_{0})$

under the initial condition $U(t_{0},t_{0})=I$. Denote the joint system-field
state by $\ket{\Psi(t)}$. In the Schrödinger picture,
$\ket{\Psi(t)}=U(t,t_{0})|\Psi(t_{0})\rangle$, which is the solution to the
stochastic Schrödinger equation

$d\ket{\Psi(t)}=\left[-\mathrm{i}H_{\rm
eff}dt+dB_{\mathrm{in}}^{\dagger}(t)L-L^{\dagger}dB_{\mathrm{in}}(t)\right]\ket{\Psi(t)},$
(2)

where $H_{\rm eff}=H-\frac{\mathrm{i}}{2}L^{\dagger}L$ is the effective
Hamiltonian that is not self-adjoint [17, Chapter 11]. In particular, if the
input fields are initially in the vacuum state $\ket{\Phi_{0}}$, then (2)
reduces to, [17, Chapter 11],

$d\ket{\Psi(t)}=\left(-\mathrm{i}H_{\rm
eff}dt+dB_{\mathrm{in}}^{\dagger}(t)L\right)\ket{\Psi(t)}.$ (3)

On the other hand, in the Heisenberg picture, the time evolution of the system
operator $X$, denoted by $\mathfrak{j}_{t}(X)\equiv
X(t)=U^{\ast}(t,t_{0})(X\otimes I_{\mathrm{field}})U(t,t_{0})$, follows the
Itô QSDE

$\displaystyle
d\mathfrak{j}_{t}(X)=\mathfrak{j}_{t}(\mathcal{L}X)dt+\sum_{k=1}^{m}dB_{\mathrm{in},k}^{\ast}(t)\mathfrak{j}_{t}([X,L_{k}])$
$\displaystyle\hskip
34.44434pt+\sum_{k=1}^{m}\mathfrak{j}_{t}([L_{k}^{\ast},X])dB_{\mathrm{in},k}(t).$

Finally, the output field annihilation operators are
$B_{\mathrm{out},k}(t)=U^{\ast}(t,t_{0})B_{\mathrm{in},k}(t)U(t,t_{0})$,
$(k=1,\ldots,m)$, whose dynamical evolution is
$dB_{\mathrm{out}}(t)=L(t)dt+dB_{\mathrm{in}}(t)$. More discussions on open
quantum systems can be found in, e.g., [17, 18, 15, 16].

###  II-B Continuous-mode single-photon states

For each input channel $k=1,\ldots,m$, the creation operator
$b_{\mathrm{in},k}^{\ast}$ generates a photon from the vacuum. Mathematically,
$|1_{k,t}\rangle\triangleq b_{\mathrm{in},k}^{\ast}(t)\ket{\Phi_{0}}$ means a
photon is generated at time $t$ in the $k$th input channel. By (1), $\langle
1_{j,t}|1_{k,\tau}\rangle=\delta_{jk}\delta(t-\tau)$. Hence
$\\{|1_{k,t}\rangle:t\in\mathbb{R}\\}$ is an orthogonal basis of single-photon
states for each channel $k$. Indeed, a single-photon state with temporal pulse
shape $\xi(t)$ in the $k$th channel can be viewed as a superposition of a
continuum of $|1_{k,t}\rangle$, i.e.,

$|1_{\xi}\rangle=\int_{-\infty}^{\infty}\xi(t)|1_{k,t}\rangle
dt=\int_{-\infty}^{\infty}\xi(t)b_{\mathrm{in},k}^{\ast}(t)dt\ket{\Phi_{0}}.$
(4)

Physically, as a quantum state, $|1_{\xi}\rangle$ can be interpreted in the
following way: the probability of finding the photon in the time bin
$[t,t+dt)$ is $|\xi(t)|^{2}dt$. The normalization condition $\langle
1_{\xi}|1_{\xi}\rangle=1$ requires $\int_{-\infty}^{\infty}|\xi(t)|^{2}dt=1$.
As the single-photon state $|1_{\xi}\rangle$ is parameterized by an $L^{2}$
integrable function $\xi(t)$ over $\mathbb{C}$, it is called a continuous-mode
single-photon state [21, 22, 23, 24, 25, 26, 27, 28]. If a single photon is
superposed over $m$ channels, then the single-photon state is a superposition
state of the form

$|1_{\zeta}\rangle=\sum_{k=1}^{m}\int_{-\infty}^{\infty}\zeta_{k}(t)|1_{k,t}\rangle
dt,$

whose normalization condition is
$\sum_{k=1}^{m}\int_{-\infty}^{\infty}|\zeta_{k}(t)|^{2}dt=1$. In particular,
if $\zeta_{j}(t)\equiv 0$ for some $j=1,\ldots,m$, then it means that the
input channel $j$ is in the vacuum state and the photon is superposed over the
other input channels.

##  III The Tavis-Cummings model

We first present the Tavis-Cummings model in subsection III-A. In subsection
III-B, assuming that the atoms are initially in the ground state, the cavity
is empty and the input field is in the vacuum state, we present an associated
linear model. Physical interpretation of the linear model is discussed in
subsection III-C.

###  III-A The Tavis-Cummings model

Figure 1: Schematic of the Tavis-Cummings model. $N$ two-level atoms are
coupled to a single-mode cavity, which is driven by an input field. (Here, the
atoms are inside the cavity. However, in some physical implementations atoms
may be outside the cavity, but of course they must couple to the cavity; see
e.g., [9, Fig. 1]). The cavity can also be a transmission line in a
superconducting platform.

In the Tavis-Cummings model as shown in Fig. 1, the $N$ TLSs are not directly
coupled to each other; instead, they all couple to the common single-mode
cavity. The inherent system Hamiltonian of the Tavis-Cummings model is ([2,
(2.1)], [3, (1)], [4, (1)],[5, (1)],[9, (C7)],[6, (1)])

$H_{\mathrm{TC}}=\omega_{r}a^{\ast}a+\sum_{j=1}^{N}\left[\frac{\omega_{j}}{2}\sigma_{z,j}+\Gamma_{j}(a^{\ast}\sigma_{-,j}+\sigma_{+,j}a)\right].$
(5)

Here, $a$, $a^{\ast}$ denote the annihilation and creation operators of the
cavity mode satisfying $[a,a^{\ast}]=I$, $\omega_{r}$ is the frequency
detuning between the cavity mode and the input field. The two-level atom $j$
is coupled to the cavity with coupling strength $\Gamma_{j}$, $j=1,\ldots,N$,
which is assumed to a real number but can be negative [5, Table 1] or [7,
(2)]. The corresponding detuning between the transition frequency of the two-
level atom $j$ and the carrier frequency of the input field is denoted by
$\omega_{j}$. The lowering and raising operators of the two-level atom $j$ are
$\sigma_{-,j}=\ket{g_{j}}\bra{e_{j}}$ and
$\sigma_{+,j}=\ket{e_{j}}\bra{g_{j}}$, respectively. The Pauli $Z$ operator is
$\sigma_{z,j}=\sigma_{+,j}\sigma_{-,j}-\sigma_{-,j}\sigma_{+,j}$. The system
exchanges information with its environment by means of absorbing and emitting
photons, which is realized by the coupling operator $L=\sqrt{\kappa}a$. By the
development in subsection II-A, the Itô QSDEs for the Tavis-Cummings model in
the Heisenberg picture are

$\left\\{\begin{aligned}
d\sigma_{-,1}(t)&=-\mathrm{i}\omega_{1}\sigma_{-,1}(t)dt+\mathrm{i}\Gamma_{1}\sigma_{z,1}(t)a(t)dt,\\\
&\ \vdots\\\
d\sigma_{-,N}(t)&=-\mathrm{i}\omega_{N}\sigma_{-,N}(t)dt+\mathrm{i}\Gamma_{N}\sigma_{z,N}a(t)dt,\\\
da(t)&=-(\mathrm{i}\omega_{r}+\frac{\kappa}{2})a(t)dt\\\ &\hskip
8.61108pt-\mathrm{i}\sum_{j=1}^{N}\Gamma_{j}\sigma_{-,j}(t)dt-\sqrt{\kappa}dB_{\mathrm{in}}(t),\\\
dB_{\mathrm{out}}(t)&=\sqrt{\kappa}a(t)dt+dB_{\mathrm{in}}(t),\ t\geq
t_{0},\end{aligned}\right.$ (6)

which is bilinear.

###### Remark III.1

When all atoms are resonant with the resonator, i.e.,
$\omega_{1}=\cdots=\omega_{N}=\omega_{r}\equiv\omega_{s}$, rotations
$\sigma_{-,j}(t)\to e^{{\rm i}\omega_{s}t}\sigma_{-,j}(t)$, $a(t)\to e^{{\rm
i}\omega_{s}t}a(t)$, $B_{\mathrm{in}}(t)\to e^{{\rm
i}\omega_{s}t}B_{\mathrm{in}}(t)$, and $B_{\mathrm{out}}(t)\to e^{{\rm
i}\omega_{s}t}B_{\mathrm{out}}(t)$ convert system (6) to

$\left\\{\begin{aligned}
d\sigma_{-,1}(t)&=\mathrm{i}\Gamma_{1}\sigma_{z,1}(t)a(t)dt,\\\ \vdots&\\\
d\sigma_{-,N}(t)&=\mathrm{i}\Gamma_{N}\sigma_{z,N}a(t)dt,\\\
da(t)&=-\frac{\kappa}{2}a(t)dt-\mathrm{i}\sum_{j=1}^{N}\Gamma_{j}\sigma_{-,j}(t)dt-\sqrt{\kappa}dB_{\mathrm{in}}(t),\\\
dB_{\mathrm{out}}(t)&=\sqrt{\kappa}a(t)dt+dB_{\mathrm{in}}(t),\ t\geq
t_{0}.\end{aligned}\right.$ (7)

In other words, when all the atoms and the cavity are resonant, the
frequencies do not affect the system dynamics, except the central frequency of
the radiation field.

###  III-B The corresponding linear model

Assume that all the two-level atoms are initially in the ground state and the
cavity is in the vacuum state $\ket{0}$. That is, the initial state of the
Tavis-Cummings model is

$|\zeta\rangle=|g_{1}g_{2}\cdots g_{N}\rangle\otimes|0\rangle.$ (8)

Let

$X(t)=\left[\begin{array}[]{ccccc}\sigma_{-,1}(t)&\sigma_{-,2}(t)&\cdots&\sigma_{-,N}(t)&a(t)\\\
\end{array}\right]^{\top}.$ (9)

Notice that

$\sigma_{z,j}\ket{\zeta\Phi_{0}}=-\ket{\zeta\Phi_{0}}$ (10)

for all $j=1,\ldots,N$. From (6) we get

$\displaystyle dX(t)\ket{\zeta\Phi_{0}}$
$\displaystyle=AX(t)\ket{\zeta\Phi_{0}}dt+BdB_{\mathrm{in}}(t)\ket{\zeta\Phi_{0}},$
(11) $\displaystyle dB_{\mathrm{out}}(t)\ket{\zeta\Phi_{0}}$
$\displaystyle=CX(t)\ket{\zeta\Phi_{0}}dt+dB_{\mathrm{in}}(t)\ket{\zeta\Phi_{0}},$

where

$\displaystyle
A=-\mathrm{i}\left[\begin{array}[]{ccccc}\omega_{1}&0&\cdots&0&\Gamma_{1}\\\
0&\omega_{2}&\cdots&0&\Gamma_{2}\\\ \vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&\cdots&\omega_{N}&\Gamma_{N}\\\
\Gamma_{1}&\Gamma_{2}&\cdots&\Gamma_{N}&\omega_{r}-\frac{\kappa\mathrm{i}}{2}\end{array}\right],$
(12) $\displaystyle
B=\left[\begin{array}[]{ccccc}0&0&\cdots&0&-\sqrt{\kappa}\end{array}\right]^{\top},~{}~{}C=-B^{\top}.$

(11) is a linear system. Actually, a linear quantum system of $N+1$ quantum
harmonic oscillators $\bar{a}=[a_{1},\ldots,a_{N},a]^{\top}$ with system
Hamiltonian
$H=\omega_{r}a^{\ast}a+\sum_{j=1}^{N}\left[\omega_{j}a_{j}^{\ast}a_{j}+\Gamma_{j}(a^{\ast}a_{j}+a_{j}^{\ast}a)\right]$
and coupling operator $L=\sqrt{\kappa}a$ has the following linear Itô QSDEs

$\displaystyle d\bar{a}(t)$
$\displaystyle=A\bar{a}(t)dt+BdB_{\mathrm{in}}(t),$ (13) $\displaystyle
dB_{\mathrm{out}}(t)$ $\displaystyle=C\bar{a}(t)dt+dB_{\mathrm{in}}(t),$

where $A,B,C$ are exactly those in (12). More discussions on linear quantum
systems theory can be found in e.g., [29, 30, 31]. The transfer function of
the linear quantum system (13) is

$G[s]=1+C(sI-A)^{-1}B.$ (14)

Plugging (12) into (14) yields

$\displaystyle
G[s]=\frac{\displaystyle{\sum_{k=1}^{N}}\left(\Gamma_{k}^{2}+\frac{1}{N}(s+\mathrm{i}\omega_{k})(s+\mathrm{i}\omega_{r}-\frac{\kappa}{2})\right)\prod_{j\neq
k}^{N}(s+\mathrm{i}\omega_{j})}{\displaystyle{\sum_{k=1}^{N}}\left(\Gamma_{k}^{2}+\frac{1}{N}(s+\mathrm{i}\omega_{k})(s+\mathrm{i}\omega_{r}+\frac{\kappa}{2})\right)\prod_{j\neq
k}^{N}(s+\mathrm{i}\omega_{j})}.$ (15)

If $\omega_{1}=\cdots=\omega_{N}\equiv\omega_{s}$, (15) reduces to

$G[s]=\frac{(\sqrt{N}\bar{\Gamma})^{2}+(s+\mathrm{i}\omega_{s})(s+\mathrm{i}\omega_{r}-\frac{\kappa}{2})}{(\sqrt{N}\bar{\Gamma})^{2}+(s+\mathrm{i}\omega_{s})(s+\mathrm{i}\omega_{r}+\frac{\kappa}{2})},$
(16)

where $\bar{\Gamma}\triangleq\sqrt{\frac{1}{N}\sum_{k=1}^{N}\Gamma_{k}^{2}}$.
Let $\omega_{r}=\omega_{s}=0$, i.e., all atoms are resonant with the cavity
resonator. Define $T[s]\triangleq G[s]-1$. Then

$|T[{\rm
i}\omega]|^{2}=\frac{\kappa^{2}\omega^{2}}{[(\sqrt{N}\bar{\Gamma})^{2}-\omega^{2}]^{2}+\frac{\kappa^{2}}{4}\omega^{2}}.$
(17)

Clearly, $|T[0]|^{2}=0$ and $|T[{\rm i}\omega]|^{2}$ has two peaks attained at
$\omega=\pm\sqrt{N}\bar{\Gamma}$, respectively. (If $N=0$, then $|T[{\rm
i}\omega]|^{2}$ has only one peak attained at $\omega=0$, which is the empty
cavity case.)

###### Remark III.2

Suppose energy enters the system via the input field $B_{\rm{in}}$ and flows
out through the output field $B_{\rm{out}}$. $T[s]$ is related to the energy
stored in the system. In fact, define the performance variable $z\triangleq
C\bar{a}$. Then $T[s]$ is the transfer function from $B_{\rm{in}}$ to $z$. It
is easy to see that

$\left[\begin{array}[]{cc}A+A^{\dagger}+C^{\dagger}C&B-C^{\dagger}\\\
B^{\dagger}-C&0\end{array}\right]=0.$ (18)

Thus, by the positive real lemma in [33, Theorem 3], the system (13) with the
performance variable $z$ is passive. In subsection III-C, we will show that
$T[s]$ also reflects the vacuum Rabi mode splitting of the Tavis-Cummings
system (6). Hence, $T[s]$ bridges the passivity of the linear quantum system
(13) and the vacuum Rabi mode splitting phenomenon exhibited by the Tavis-
Cummings system (6).

In the following, we perform structural decomposition on the linear quantum
system (13). We partition $\omega_{1},\ldots,\omega_{N}$ into $M$ groups
according to their degeneracies, where the $j$th degenerated frequency is
denoted by $\tilde{\omega}_{j}$, ($j=1,\ldots,M$). Let the number of elements
be $n_{j}$ for the group $j$. In particular, if $M=N$, then
$\omega_{j}\neq\omega_{k}$ for all $1\leq j<k\leq N$. For convenience, we can
arrange the elements of $\bar{a}(t)$ in (13) so that the matrix $A$ in (12) is
of the form

$\displaystyle\tilde{A}$ $\displaystyle=$
$\displaystyle-\mathrm{i}\left[\begin{array}[]{cccccccc}\tilde{\omega}_{1}&\cdots&0&\cdots&0&\cdots&0&\Gamma_{1_{1}}\\\
\vdots&\ddots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&\cdots&\tilde{\omega}_{1}&\cdots&0&\cdots&0&\Gamma_{1_{n_{1}}}\\\
\vdots&\ddots&\vdots&\ddots&\vdots&\ddots&\vdots&\vdots\\\
0&\cdots&0&\cdots&\tilde{\omega}_{M}&\cdots&0&\Gamma_{M_{1}}\\\
\vdots&\ddots&\vdots&\ddots&\vdots&\ddots&\vdots&\vdots\\\
0&\cdots&0&\cdots&0&\cdots&\tilde{\omega}_{M}&\Gamma_{M_{n_{M}}}\\\
\Gamma_{1_{1}}&\cdots&\Gamma_{1_{n_{1}}}&\cdots&\Gamma_{M_{1}}&\cdots&\Gamma_{M_{n_{M}}}&\omega_{r}-\frac{\kappa\mathrm{i}}{2}\\\
\end{array}\right].$

In other words, we group the elements of $\bar{a}(t)$ according to the
partition of the detuned frequencies. It is easy to see that matrices $B$ and
$C$ remain the same under this re-arrangement.

###### Lemma III.1

Partition $\omega_{1},\ldots,\omega_{N}$ into groups as described above. There
is an orthogonal matrix $\tilde{T}$ that transforms the linear quantum system
(13) with system matrices $(\tilde{A},B,C)$ to another one, denoted $\Sigma$,
with system matrices

$\displaystyle\hat{A}=\tilde{T}^{\top}\tilde{A}\tilde{T}$ (19)
$\displaystyle=$
$\displaystyle-\mathrm{i}\left[\begin{array}[]{cccccc}\tilde{\omega}_{1}I_{n_{1}-1}&\mathbf{0}&\cdots&\cdots&\cdots&\mathbf{0}\\\
\mathbf{0}&\tilde{\omega}_{1}&\cdots&\cdots&\cdots&\sqrt{\tilde{\Gamma}_{1}}\\\
\vdots&\vdots&\ddots&\cdots&\cdots&\vdots\\\
\vdots&\vdots&\vdots&\tilde{\omega}_{M}I_{n_{M}-1}&\cdots&\mathbf{0}\\\
\vdots&\vdots&\vdots&\vdots&\tilde{\omega}_{M}&\sqrt{\tilde{\Gamma}_{M}}\\\
\mathbf{0}&\sqrt{\tilde{\Gamma}_{1}}&\cdots&\mathbf{0}&\sqrt{\tilde{\Gamma}_{M}}&\omega_{r}-\frac{\kappa\mathrm{i}}{2}\\\
\end{array}\right],$
$\displaystyle\hat{B}=\tilde{T}^{\top}B=B,~{}~{}\hat{C}=C\tilde{T}=C,$

where $\tilde{\Gamma}_{j}\triangleq\sum_{k=1}^{n_{j}}\Gamma_{j_{k}}^{2}$,
$(j=1,\ldots,M)$.

###### Proof:

The proof is constructive. Define a matrix $\tilde{T}$ as

$\tilde{T}=\left[\begin{array}[]{ccccc}\tilde{T}_{1}&\mathbf{0}&\cdots&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\tilde{T}_{2}&\cdots&\mathbf{0}&\mathbf{0}\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\
\mathbf{0}&\mathbf{0}&\cdots&\tilde{T}_{M}&\mathbf{0}\\\ 0&0&\cdots&0&1\\\
\end{array}\right],$ (20)

where for each $j=1,\ldots,M$, $\tilde{T}_{j}=[T_{j_{1}}\ T_{j_{2}}\ \cdots\
T_{j_{n_{j}}}]\in\mathbb{R}^{n_{j}\times n_{j}}$, in which

$\displaystyle T_{j_{1}}=$ $\displaystyle\
\sqrt{\frac{\Gamma_{j_{2}}^{2}}{\Gamma_{j_{1}}^{2}+\Gamma_{j_{2}}^{2}}}\left[1~{}\frac{-\Gamma_{j_{1}}}{\Gamma_{j_{2}}}~{}0~{}\cdots~{}0\right]^{\top},$
$\displaystyle T_{j_{2}}=$ $\displaystyle\
\sqrt{\frac{\Gamma_{j_{2}}^{2}\Gamma_{j_{3}}^{2}}{(\Gamma_{j_{1}}^{2}+\Gamma_{j_{2}}^{2})(\Gamma_{j_{1}}^{2}+\Gamma_{j_{2}}^{2}+\Gamma_{j_{3}}^{2})}}$
$\displaystyle\times\left[\frac{\Gamma_{j_{1}}}{\Gamma_{j_{2}}}~{}1~{}\frac{-(\Gamma_{j_{1}}^{2}+\Gamma_{j_{2}}^{2})}{\Gamma_{j_{2}}\Gamma_{j_{3}}}~{}0~{}\cdots~{}0\right]^{\top},$
$\displaystyle\vdots$ $\displaystyle T_{j_{(n_{j}-1)}}=$ $\displaystyle\
\sqrt{\frac{\Gamma_{j_{(n_{j}-1)}}^{2}\Gamma_{j_{n_{j}}}^{2}}{(\sum_{k=1}^{n_{j}-1}\Gamma_{j_{k}}^{2})(\sum_{k=1}^{n_{j}}\Gamma_{j_{k}}^{2})}}$
$\displaystyle\times\left[\frac{\Gamma_{j_{1}}}{\Gamma_{j_{(n_{j}-1)}}}~{}\frac{\Gamma_{j_{2}}}{\Gamma_{j_{(n_{j}-1)}}}~{}\cdots~{}1~{}\frac{-\sum_{k=1}^{n_{j}-1}\Gamma_{j_{k}}^{2}}{\Gamma_{j_{(n_{j}-1)}}g_{j_{n_{j}}}}\right]^{\top},$
$\displaystyle T_{j_{n_{j}}}=$ $\displaystyle\
\sqrt{\frac{1}{\sum_{k=1}^{n_{j}}\Gamma_{j_{k}}^{2}}}\left[\Gamma_{j_{1}}~{}\Gamma_{j_{2}}~{}\cdots~{}\Gamma_{j_{(n_{j}-1)}}~{}\Gamma_{j_{n_{j}}}\right]^{\top}.$

It can be easily verified that $\tilde{T}$ is orthogonal. Moreover, simple
algebraic manipulations yield that
$\tilde{T}^{\top}\tilde{A}\tilde{T}=\hat{A}$. ∎

The transformed linear quantum system with system matrices $(\hat{A},B,C)$ has
a nice structure. Denote by $b_{j}$ the system coordinate corresponding to the
row whose last entry is $-\mathrm{i}\sqrt{\bar{\Gamma}_{j}}$ in the matrix
$\hat{A}$, $j=1,\ldots,M$. Then from the structure of $(\hat{A},B,C)$ it can
be easily seen that this system has a subsystem of the form

$\left\\{\begin{aligned} db_{1}(t)=&\
-\mathrm{i}\tilde{\omega}_{1}b_{1}(t)dt-\mathrm{i}\sqrt{\bar{\Gamma}_{1}}a(t)dt,\\\
\vdots&\\\ db_{M}(t)=&\
-\mathrm{i}\tilde{\omega}_{M}b_{M}(t)dt-\mathrm{i}\sqrt{\bar{\Gamma}_{M}}a(t)dt,\\\
da(t)=&-(\mathrm{i}\omega_{r}+\frac{\kappa}{2})a(t)dt\\\
&-\mathrm{i}\sum_{j=1}^{M}\sqrt{\bar{\Gamma}_{j}}b_{j}(t)dt-\sqrt{\kappa}dB_{\mathrm{in}}(t),\\\
dB_{\mathrm{out}}(t)=&\sqrt{\kappa}a(t)dt+dB_{\mathrm{in}}(t).\end{aligned}\right.$
(21)

The other $M-1$ subsystems are all isolated systems, which are called
decoherence-free subsystems (DFSs) in the linear quantum control literature
[34, 29, 30, 31]. Of course, if $n_{j}=1$ for some $j=1,\ldots,M$, there is no
such subsystem as can be seen clearly from the matrix $\hat{A}$ in (19).
According to quantum linear systems theory, the subsystem (21) is a
controllable and observable subsystem.

The following result is an immediate consequence of Lemma III.1.

###### Corollary III.1

If $\omega_{1}=\cdots=\omega_{N}\equiv\omega_{s}$, then $M=1$ and $n_{1}=N$ in
(19). Accordingly, the transformation matrix $\tilde{T}$ in (20) reduces to

$\left[\begin{array}[]{c|c}\tilde{T}_{1}&\mathbf{0}\\\ 0&1\\\
\end{array}\right]\equiv T,$ (22)

which transforms the quantum linear system (13) to a new one with system
matrices

$\displaystyle\hat{A}=T^{\top}AT=\left[\begin{array}[]{c|c}\hat{A}_{\bar{c}\bar{o}}&\mathbf{0}\\\
\hline\cr\mathbf{0}&\hat{A}_{co}\\\ \end{array}\right],$ (23)
$\displaystyle\hat{B}=T^{\top}B=\left[\begin{array}[]{c}\hat{B}_{\bar{c}\bar{o}}\\\
\hline\cr\hat{B}_{co}\\\
\end{array}\right],\hat{C}=CT=[\hat{C}_{\bar{c}\bar{o}}\ |\ \hat{C}_{co}],$

where

$\displaystyle\hat{A}_{\bar{c}\bar{o}}=-\mathrm{i}\omega_{s}I_{N-1},\ \
\hat{A}_{co}=-\left[\begin{array}[]{cc}\mathrm{i}\omega_{s}&\mathrm{i}\sqrt{N}\bar{\Gamma}\\\
\mathrm{i}\sqrt{N}\bar{\Gamma}&\mathrm{i}\omega_{r}+\frac{\kappa}{2}\end{array}\right]$
(24) $\displaystyle\hat{B}_{\bar{c}\bar{o}}=\mathbf{0},\
\hat{B}_{co}=-\left[\begin{array}[]{c}0\\\
\sqrt{\kappa}\end{array}\right],\hat{C}_{\bar{c}\bar{o}}=\mathbf{0},\
\hat{C}_{co}=\left[\begin{array}[]{cc}0&\sqrt{\kappa}\end{array}\right].$

###### Remark III.3

In the linear quantum systems theory [34, 29, 30, 31] DFSs are defined in the
Heisenberg picture. On the other hand, decoherence-free subspaces widely used
in the quantum information community are characterized by density matrices in
the Schrödinger picture; see references 3-8 in [35]. It can be seen that the
eigenvectors $\\{T_{1},\ldots,T_{N-1}\\}$ of the DFS with system matrices
$(\hat{A}_{\bar{c}\bar{o}},\hat{B}_{\bar{c}\bar{o}},\hat{C}_{\bar{c}\bar{o}})$
constitute a basis of the decoherence-free subspace. Some examples will be
given in subsection III-C.

###### Remark III.4

The coupling operator $L=\sqrt{\kappa}a$ is called the noise operator in [35],
through which the system information leaks irreversibly into its surrounding
environment. It can be clearly seen from (23)-(24) that $L$ has no effect on
the DFS with system matrices
$(\hat{A}_{\bar{c}\bar{o}},\hat{B}_{\bar{c}\bar{o}},\hat{C}_{\bar{c}\bar{o}})$,
thus the DFS is robust with respect to the cavity decay rate $\kappa$. This is
the so-called $\gamma$-robustness in [35, Definition 3]. Moreover, from (24)
it can be seen that the DFS is robust with respect to the variations of the
atom-cavity coupling strengths $\Gamma_{j}$ as well. Finally, the DFS is not
attractive from [35, Proposition 3], which is for noiseless subsystems and
reduces to DFSs when ${\cal H}_{F}\simeq\mathbb{C}$; see [35, Definition 8 and
subsection II-A] for details.

###  III-C Physical interpretation

The discussions of the Tavis-Cummings model in the previous subsection in
terms of linear quantum systems theory appear purely mathematical; however,
these results can indeed reveal several typical features of the Tavis-Cummings
model.

Firstly, from (16), it can be seen that the collection of the $N$ atoms act as
a giant atom of detuned frequency $\omega_{s}$ and coupling strength
$\sqrt{N}\bar{\Gamma}$, which reflects the $\sqrt{N}$\- scaling of the
collective coupling strength; in other words, this giant atom decays $N$ times
as fast as a single atom. This is the physical basis of superradiance.

Secondly, the two peaks of the transfer function $|T[\mathrm{i}\omega]|^{2}$
in (17) at $\omega=\pm\sqrt{N}\bar{\Gamma}$ echo the vacuum Rabi mode
splitting in the Tavis-Cummings model, see, e.g., [4, Fig. 2(b)], [5, Fig. 3],
and [6, Fig. 2(b)].

Thirdly, the controllable and observable subsystem echos the bright states and
the uncontrollable and unobservable subsystems echo the dark states of the
Tavis-Cummings model. For simplicity, we look at the simplest case of
$\omega_{1}=\cdots=\omega_{N}\equiv\omega_{s}$. In this case, according to
Corollary III.1, there is an orthogonal matrix $T$ which yields a system with
system matrices given in (23). Let $N=3$ and assume the coupling constants
$\Gamma_{1}=-\Gamma_{2}=\Gamma_{3}$. (In [5],
$\Gamma_{1},\Gamma_{2},\Gamma_{3}$ are respectively
$g_{\mathrm{A}},g_{\mathrm{B}},g_{\mathrm{C}}$. As shown in [5, Table 1], the
actual values of $\Gamma_{1},-\Gamma_{2},\Gamma_{3}$ are not exactly
identical, but the discrepancy has negligible effect as can be seen from
consistency between the red region (for real data) and the white dashed curves
(for theoretical calculation) in [5, Fig. 3].) The orthogonal transformation
matrix $T$ can be calculated as
$T_{1}=\frac{1}{\sqrt{2}}\left[\begin{array}[]{cccc}1&1&0&0\\\
\end{array}\right]^{\top}$,
$T_{2}=\frac{1}{\sqrt{6}}\left[\begin{array}[]{cccc}-1&1&2&0\\\
\end{array}\right]^{\top}$,
$T_{3}=\frac{1}{\sqrt{3}}\left[\begin{array}[]{cccc}1&-1&1&0\\\
\end{array}\right]^{\top}$, and $T_{4}=\left[\begin{array}[]{cccc}0&0&0&1\\\
\end{array}\right]^{\top}$. Identify $g$ with $0$ and $e$ with $1$
respectively. It can be verified that the dark states
$|3,1d_{1}\rangle=\frac{1}{\sqrt{2}}(\ket{e,g,g,0}-\ket{g,g,e,0})$ and
$|3,1d_{2}\rangle=\frac{1}{\sqrt{2}}(\ket{g,e,g,0}+\ket{g,g,e,0})$ in [5] can
be expressed as $|3,1d_{1}\rangle=\frac{1}{2}T_{1}-\frac{\sqrt{3}}{2}T_{2}$,
and $|3,1d_{2}\rangle=\frac{1}{2}T_{1}+\frac{\sqrt{3}}{2}T_{2}$ respectively,
while the bright states
$|3,1\pm\rangle=\frac{1}{\sqrt{2}}\ket{g,g,g,1}\pm\frac{1}{\sqrt{6}}\ket{e,g,g,0}$
in [5] can be written as $|3,1\pm\rangle=\frac{1}{\sqrt{2}}(T_{4}\pm T_{3})$.
In other words, the dark states $|3,1d_{1}\rangle$ and $|3,1d_{2}\rangle$ live
in the decoherence-free subspace spanned by $T_{1}$ and $T_{2}$, while the
bright states $|3,1\pm\rangle$ live in the controllable and observable
subspace spanned by $T_{3}$ and $T_{4}$.

Take the case $N=2$ for another example. Similar correspondence can be found
between the eigenstates of the Tavis-Cummings model (6) and the vectors in
(22). Regard the states $|g,e,0\rangle$, $|e,g,0\rangle$, $|g,g,1\rangle$ in
[9] as vectors $[0~{}1~{}0]^{\top}$, $[1~{}0~{}0]^{\top}$, and
$[0~{}0~{}1]^{\top}$, respectively. In the resonant ($\omega_{r}=\omega_{s}$)
case, the dark state
$|0\rangle_{r3}=\frac{1}{\sqrt{N}\bar{\Gamma}}(\Gamma_{1}\ket{g,e,0}-\Gamma_{2}\ket{e,g,0})$
in [9, Appendix C-3] is $-T_{1}$, and the two bright states
$|\pm\rangle_{r3}=\frac{1}{N\bar{\Gamma}}(\Gamma_{2}\ket{g,e,0}+\Gamma_{1}\ket{e,g,0}\pm\sqrt{N}\bar{\Gamma}\ket{g,g,1})$
in [9, Appendix C-3] are $\frac{1}{\sqrt{2}}(T_{2}\pm T_{3})$. Also,
$|0\rangle_{r3}$ and $|\pm\rangle_{r3}$ are respectively $-\ket{2,1d}$ and
$-\ket{2,1\pm}$ in [5] when $\Gamma_{1}=-\Gamma_{2}=\Gamma_{3}$. Moreover, in
the dispersive regime ($\Delta_{r}=|\omega_{r}-\omega_{s}|\gg\Gamma_{j}$,
$j=1,\ldots,N$) considered in [9, Appendix C], the three states
$|+^{\prime}\rangle_{r3},|-^{\prime}\rangle_{r3},|1^{\prime}\rangle_{r3}$
given in [9, (C18)] can be expressed by $|+^{\prime}\rangle_{r3}=-T_{1}$,
$|-^{\prime}\rangle_{r3}=\frac{1}{\sqrt{N\bar{\Gamma}^{2}+\Delta_{r}^{2}}}(\Delta_{r}T_{2}-\sqrt{N}\bar{\Gamma}T_{3})\approx
T_{2}$, and
$|1^{\prime}\rangle_{r3}=\frac{1}{\sqrt{N\bar{\Gamma}^{2}+\Delta_{r}^{2}}}(\sqrt{N}\bar{\Gamma}T_{2}+\Delta_{r}T_{3})\approx
T_{3}$. Thus, in both cases, the dark states live in the space spanned by
$T_{1}$ while the bright states live in the space spanned by $T_{2}$ and
$T_{3}$.

Fourthly, the dynamics of two remote spin ensembles coupled by a cavity bus
are experimentally studied in [7]. The Hamiltonian of the system is given in
[7, (2)], which is of the form of $H_{\mathrm{TC}}$ for the Tavis-Cummings
model. When $\varphi\approx 48.1^{\mathrm{o}}$ in [7, Fig. 2(b)], both spin
ensembles resonate with the cavity mode. Thus, this particularly interesting
case can be analyzed by Corollary III.1. Indeed, the eigenstates $\ket{\pm}$
in [7, (3)] and the dark mode $\ket{D}$ in [7, (4)] can be obtained by means
of the orthonormal matrix $T$ in Corollary III.1.

Fifthly, in subsections IV-A and IV-B we study how the Tavis-Cummings model
(6) responses to a continuous-mode single-photon input state, where the linear
model developed in subsection III-B plays an essential role.

Finally, the single-excitation superradiant and subradiant states of the
Tavis-Cummings model can be analyzed by means of the quantum linear systems
theory presented in subsection III-B; see Remark IV.2 in subsection IV-C.

##  IV The single-excitation case

In this section, we investigate the dynamics of the Tavis-Cummings model when
there is only one excitation.

###  IV-A Response to single-photon inputs

In this subsection, we derive an analytic expression of the steady-state
output field state when the Tavis-Cummings model is initialized in the state
$\ket{\zeta}$ given in (8), and driven by a single-photon state.

We start with the following lemma which discusses the controllability ([34,
Sec. III-B], [29, Definition 1]) of the passive linear quantum system (13).

###### Lemma IV.1

The passive linear quantum system (13) is controllable if and only if
$\omega_{j}\neq\omega_{k}$ for all $1\leq j<k\leq N$.

###### Proof:

Let a scalar $\lambda\in\mathbb{C}$ and a vector
$x=[x_{1},\ldots,x_{N+1}]^{\top}\in\mathbb{C}^{N+1}$ satisfy

$Ax=\lambda x,\ \ x^{\dagger}B=0.$ (25)

We have

$x_{N+1}=0,\ \ \ -\mathrm{i}\omega_{k}x_{k}=\lambda x_{k},\ \ \ k=1,\ldots,N,$
(26)

and

$\sum_{k=1}^{N}\Gamma_{k}x_{k}=0.$ (27)

_Necessity_. We prove it by contradiction. Without loss of generality, assume
$\omega_{1}=\omega_{2}$. We choose $\lambda=-\mathrm{i}\omega_{1}$,
$x_{1}=\Gamma_{2}$, $x_{2}=-\Gamma_{1}$, and $x_{j}=0$ for all
$j=3,\ldots,N+1$. Clearly, the scalar $\lambda$ and the nonzero vector $x$
satisfy (25). This means that the system is not controllable. A contradiction
is reached.

_Sufficiency_. If $x_{j}\neq 0$ for some $j=1,\dots,N$, then by (26)
$\lambda=-\mathrm{i}\omega_{j}$. As $\omega_{j}\neq\omega_{k}$ for all $1\leq
j<k\leq N$, again by (26) $x_{k}=0$ for all $j\neq k$. Thus, (27) reduces to
$\Gamma_{j}x_{j}=0$, which yields $x_{j}=0$. This shows that any vector $x$
satisfying (25) must be a zero vector. Hence, the system is controllable. ∎

For system (21), $\tilde{\omega}_{j}\neq\tilde{\omega}_{k}$ for all $1\leq
j<k\leq M$. By Lemma IV.1, system (21) is Hurwitz stable.

With the aid of Lemma IV.1, the main result of this section can be derived.

###### Theorem IV.1

Assume that the Tavis-Cummings model (6) is initialized in the state
$\ket{\zeta}=|g_{1}g_{2}\cdots g_{N}\rangle\otimes|0\rangle$ and driven by a
single-photon input state with pulse shape $\xi$. The steady-state
($t\to\infty$ and $t_{0}\to-\infty$) output field state is a single-photon
state with the frequency-domain pulse shape

$\eta[\mathrm{i}\omega]=G[\mathrm{i}\omega]\xi[\mathrm{i}\omega],$ (28)

where the transfer function $G[s]$ is given by (15).

###### Proof:

By system (6) and (10), we have

$\displaystyle\bra{\zeta\Phi_{0}}b_{\rm out}(t)$ $\displaystyle=$
$\displaystyle Ce^{A(t-t_{0})}\bra{\zeta\Phi_{0}}X(t_{0})$
$\displaystyle+\int_{t_{0}}^{t}Ce^{A(t-\tau)}\bra{\zeta\Phi_{0}}b_{\rm
in}(\tau)d\tau+\bra{\zeta\Phi_{0}}b(t).$

If $A$ is Hurwitz stable, then $Ce^{A(t-t_{0})}\bra{\zeta\Phi_{0}}X(t_{0})\to
0$ as $t_{0}\to-\infty$. If $A$ is not Hurwitz stable, then by Corollary
III.1, the $\bar{c}\bar{o}$ subsystem does not affect the input-output
behavior, while the $co$ subsystem is Hurwitz stable; see Proposition VII.1 in
the Appendix. Hence we also have
$Ce^{A(t-t_{0})}\bra{\zeta\Phi_{0}}X(t_{0})\to 0$ as $t_{0}\to-\infty$. As a
result, sending $t_{0}\to-\infty$, we get

$\bra{\zeta\Phi_{0}}b_{\rm
out}(t)=\int_{-\infty}^{\infty}g_{G}(t-r)\bra{\zeta\Phi_{0}}b_{\rm in}(r)dr,$
(29)

where $g_{G}(t)$ is the impulse response function associated to the transfer
function $G[s]$ in (14). As $b_{\rm out}(t)=U^{\ast}b(t,t_{0})_{\rm
in}(t)U(t,t_{0})$, we have

$\bra{\zeta\Phi_{0}}b_{\rm
in}(t)=\int_{-\infty}^{\infty}g_{G}(t-r)\bra{\zeta\Phi_{0}}b^{-}(r,-\infty)dr,$
(30)

where $b^{-}(t,t_{0})\triangleq U(t,t_{0})b_{\rm in}(t)U^{\ast}(t,t_{0})$.
Then following the stable inverse technique in [24, Lemma 1], we can get

$\displaystyle b^{-\ast}(r,-\infty)\ket{\zeta\Phi_{0}}$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}g_{G^{-1}}^{\ast}(r-t)b_{\rm
in}^{\ast}(t)dt\ket{\zeta\Phi_{0}}$ (31) $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}g_{G}(t-r)b_{\rm
in}^{\ast}(t)dt\ket{\zeta\Phi_{0}},$

and thus

$\int_{-\infty}^{\infty}\xi(r)b^{-\ast}(r,-\infty)\ket{\zeta\Phi_{0}}dr=\int_{-\infty}^{\infty}\eta(t)b_{\rm
in}^{\ast}(t)\ket{\zeta\Phi_{0}}dt,$ (32)

where $\eta(t)$ is the time-domain counterpart of $\eta[\mathrm{i}\omega]$ in
(28). Consequently, the steady-state joint system-field state in [24, Eq.
(85)] is

$\rho_{\infty}=\ket{\zeta}\bra{\zeta}\otimes\ket{1_{\eta}}\bra{1_{\eta}}.$
(33)

The steady-state output field state $\rho_{\rm out}$ is then obtained by
tracing over the initial system state, i.e., $\rho_{\rm out}=\mathrm{Tr}_{\rm
sys}[\rho_{\infty}]=\ket{1_{\eta}}\bra{1_{\eta}}$, which is a pure state
$\ket{1_{\eta}}$. ∎

###  IV-B Simulation results for $N\leq 4$

In this subsection, we illustrate Theorem IV.1 by a special case of $N\leq 4$,
i.e., at most four two-level atoms are coupled to the cavity. The single-
photon input state is supposed to have a rising exponential pulse shape

$\xi(t)=\left\\{\begin{array}[]{cc}\sqrt{\gamma}e^{\frac{\gamma}{2}t},&t\leq
0,\\\ 0,&t>0,\end{array}\right.$ (34)

where $\gamma$ denotes the full width at half maximum (FWHM) of the Lorentzian
spectrum. The input and output photon probability distributions are shown in
Fig. 2. It is worthwhile to notice that the carrier frequency of the single-
photon input field is not shown in (34), the reason is that all the
frequencies in the system Hamiltonian $H_{\mathrm{TC}}$ in (5) are detuned
from this carrier frequency.

Figure 2: The input and output photon probability distributions. The input
photon probability distribution is plotted as the green curve; the red dot
dashed curve plots the output photon probability distribution in the case that
all the three two-level atoms are resonant with each other
($\omega_{1}=\omega_{2}=\omega_{3}=0$); the non-resonant case is plotted as
the blue dotted curve ($\omega_{1}=1$, $\omega_{2}=-1$, $\omega_{3}=0$); the
output photon probability distributions for the resonant cases of $N=1$
(Jaynes-Cummings system), $N=2$ and $N=4$ are shown as the black solid,
magenta dashed, and black dashed curves respectively. $(\gamma=\kappa=1$,
$\omega_{r}=0$, and $\Gamma_{j}=1$, $j=1,\ldots,4.$)

In Fig. 2, it can be observed that the emitted photon is more likely to be
found when $t<0$ as the number of atoms increases, which means it interacts
with the Jaynes-Cummings system (the $N=1$ case) more easily than with the
Tavis-Cummings model. On the other hand, the Rabi oscillation at $t>0$
indicates that the photon can be repeatedly absorbed and emitted by the two-
level atoms. The oscillation also becomes stronger when there are more atoms,
which indicates that the added atoms increase the time for the photon to
escape from the cavity. The Rabi oscillation monotonically decays when all the
atoms are resonant with each other, while revivals can be observed in the non-
resonant case (the blue dotted curve). Moreover, oscillation sustains much
longer in the non-resonant case than in the resonant case. Finally, simulation
shows that in the resonant case, the shapes of the input and output pulses are
quite close to each other when $N$ is large, showing that the input single
photon hardly interacts with the atomic ensemble. This phenomenon is confirmed
by $\lim_{N\to\infty}G[{\rm i}\omega]=1$ for any fixed $\omega$, where $G[s]$
is given in (16).

###  IV-C An analytic form of the superposition state

In this subsection, an analytic form of the joint system-field state is
derived.

The following theorem is the main result of this subsection.

###### Theorem IV.2

Assume that the $N$ two-level atoms of the Tavis-Cummings model are resonant
with each other, i.e., $\omega_{1}=\cdots=\omega_{N}\equiv\omega_{s}$, the
$k$th two-level atom is in the excited state, the others are in the ground
state, the cavity is empty and the Tavis-Cummings model is driven by the
vacuum input state. That is, the initial joint system-field state is
$|\Psi_{k}(0)\rangle=\ket{g_{1}g_{2}\ldots e_{k}\ldots
g_{N}0}\otimes\ket{\Phi_{0}}$. Then the joint system-field state is

$\displaystyle\ket{\Psi_{k}(t)}=$ $\displaystyle\ c_{k}(t)|g_{1}g_{2}\cdots
e_{k}\cdots g_{N}0\Phi_{0}\rangle$ $\displaystyle+\sum_{j\neq
k}^{N}c_{j}(t)|g_{1}g_{2}\cdots e_{j}\cdots g_{N}0\Phi_{0}\rangle$
$\displaystyle+c_{N+1,k}(t)\int_{0}^{t}\varphi(\tau)dB_{\mathrm{in}}^{\ast}(\tau)|g_{1}g_{2}\cdots
g_{N}0\Phi_{0}\rangle$ $\displaystyle+c_{N+2,k}(t)|g_{1}g_{2}\cdots
g_{N}1\Phi_{0}\rangle,$ (35)

where

$\displaystyle c_{k}(t)=\frac{e^{\frac{N-2}{2}\mathrm{i}\omega_{s}t}}{2\chi
N\bar{\Gamma}^{2}}\Bigg{[}2\chi\sum_{j\neq
k}^{N}\Gamma_{j}^{2}+\Gamma_{k}^{2}\Big{(}\lambda_{1}e^{-\frac{\lambda_{2}t}{4}}-\lambda_{2}e^{-\frac{\lambda_{1}t}{4}}\Big{)}\Bigg{]},$
$\displaystyle
c_{j}(t)=\frac{-\Gamma_{j}\Gamma_{k}e^{\frac{N-2}{2}\mathrm{i}\omega_{s}t}}{2\chi
N\bar{\Gamma}^{2}}\Bigg{[}2\chi-\left(\lambda_{1}e^{-\frac{\lambda_{2}t}{4}}-\lambda_{2}e^{-\frac{\lambda_{1}t}{4}}\right)\Bigg{]},j\neq
k,$ $\displaystyle
c_{N+1,k}(t)=\frac{\Gamma_{k}}{\sqrt{N}\bar{\Gamma}}e^{\frac{N-2}{2}\mathrm{i}\omega_{s}t},$
$\displaystyle
c_{N+2,k}(t)=\frac{-4\mathrm{i}\Gamma_{k}\sinh\left(\frac{1}{4}\chi
t\right)}{\chi}e^{-\frac{1}{4}[\kappa+2\mathrm{i}(\omega_{r}-(N-1)\omega_{s})]t},$
$\displaystyle\varphi(\tau)=\frac{-4\mathrm{i}\bar{\Gamma}\sqrt{\kappa
N}\sinh\left(\frac{1}{4}\chi\tau\right)}{\chi}e^{-\frac{\kappa\tau+2{\rm
i}[\omega_{s}(\tau-t)+\omega_{r}\tau-\omega_{s}t]}{4}},$ $\displaystyle\hskip
25.83325pt(0\leq\tau\leq t),$ (36)

with $\lambda_{1}=\kappa+\chi+2\mathrm{i}(\omega_{r}-\omega_{s})$,
$\lambda_{2}=\kappa-\chi+2\mathrm{i}(\omega_{r}-\omega_{s})$, and
$\chi=\sqrt{(\kappa+2\mathrm{i}\omega_{r}-2\mathrm{i}\omega_{s})^{2}-16N\bar{\Gamma}^{2}}$.
Moreover, in the steady state ($t=\infty$), by ignoring the rotating term
$e^{\frac{N-2}{2}\mathrm{i}\omega_{s}t}$ in the coefficients in (36)
(equivalently, setting $\omega_{s}=0$) the superposition state (35) becomes

$\displaystyle|\Psi_{k}(\infty)\rangle=$ $\displaystyle\
c_{k}(\infty)|g_{1}g_{2}\cdots e_{k}\cdots g_{N}0\Phi_{0}\rangle$ (37)
$\displaystyle+\sum_{j\neq k}^{N}c_{j}(\infty)|g_{1}g_{2}\cdots e_{j}\cdots
g_{N}0\Phi_{0}\rangle$ $\displaystyle+c_{N+1,k}(\infty)|g_{1}g_{2}\cdots
g_{N}01_{\varphi}\rangle,$

where

$\displaystyle c_{k}(\infty)=\frac{1}{N\bar{\Gamma}^{2}}\sum_{j\neq
k}^{N}\Gamma_{j}^{2},\ \
c_{j}(\infty)=\frac{-\Gamma_{j}\Gamma_{k}}{N\bar{\Gamma}^{2}},~{}j\neq k,$
(38) $\displaystyle
c_{N+1,k}(\infty)=\frac{\Gamma_{k}}{\sqrt{N}\bar{\Gamma}}.$

The proof of Theorem IV.2 follows the recursive formula for joint system-field
states of general quantum systems; see Remark VI.3 in Section VI.

###### Remark IV.1

By Theorem IV.2, it can be seen that the cavity is eventually empty
($c_{N+2,k}(\infty)=0$), which results in a superposition state of the two-
level atoms and the output field. Moreover, by (38), the steady-state
excitation probabilities of two-level atoms are independent of the cavity
resonant frequency $\omega_{r}$ as well as the atomic transition frequency
$\omega_{s}$.

It is well-known in quantum optics that the collective radiation of an
ensemble of two-level atoms can be accelerated by superradiance or inhibited
by subradiance [1, 2, 6], the following corollary gives the steady state of
the Tavis-Cummings model initialized in a superposition state of the
superradiant and subradiant states.

###### Corollary IV.1

Assume that the $N$ two-level atoms in the Tavis-Cummings model are resonant
with each other, i.e., $\omega_{1}=\cdots=\omega_{N}\equiv\omega_{s}$. Let the
system be driven by the vacuum input state and initialized in the following
superposition state $(\alpha\ket{B_{N}}+\beta\ket{D_{N}})\otimes\ket{0}$,
where $|\alpha|^{2}+|\beta|^{2}=1$, and the single-excitation superradiant
state $\ket{B_{N}}$ and subradiant state $\ket{D_{N}}$ are respectively

$\displaystyle\ket{B_{N}}=\frac{1}{\sqrt{N}}\sum_{k=1}^{N}\ket{g_{1}\cdots
e_{k}\cdots g_{N}},$ (39)
$\displaystyle\ket{D_{N}}=\frac{1}{\sqrt{N}}\sum_{k=1}^{N}e^{-\mathrm{i}\phi_{k}}\ket{g_{1}\cdots
e_{k}\cdots g_{N}}$

with $\phi_{k}=\frac{2\pi}{N}k$. Then the joint system-field state is

$|\Psi^{\prime}(t)\rangle=\frac{1}{\sqrt{N}}\sum_{k=1}^{N}\left(\alpha+\beta
e^{-\mathrm{i}\phi_{k}}\right)|\Psi_{k}(t)\rangle,$ (40)

where $|\Psi_{k}(t)\rangle$ is given in (35). Moreover, in the steady state
($t=\infty$), by ignoring the rotating term
$e^{\frac{N-2}{2}\mathrm{i}\omega_{s}t}$ in the coefficients in (36)
(equivalently, setting $\omega_{s}=0$), the superposition state (40) becomes

$|\Psi^{\prime}(\infty)\rangle=\frac{1}{\sqrt{N}}\sum_{k=1}^{N}\left(\alpha+\beta
e^{-\mathrm{i}\phi_{k}}\right)|\Psi_{k}(\infty)\rangle,$ (41)

where $|\Psi_{k}(\infty)\rangle$ is given in (37).

The proof of Corollary IV.1 is omitted.

Let $\alpha=0$, $\beta=1$, i.e., the Tavis-Cummings model is initialized in
the pure state $|D_{N}0\rangle$. By Corollary IV.1, the joint system-field
steady state is

$|\Psi^{\prime}(\infty)\rangle=\frac{1}{\sqrt{N}}\sum_{k=1}^{N}c_{k}^{\prime}(\infty)|g_{1}\cdots
e_{k}\cdots g_{N}0\Phi_{0}\rangle,$ (42)

where

$c_{k}^{\prime}(\infty)=\frac{e^{-\mathrm{i}\phi_{k}}\sum_{j\neq
k}^{N}\Gamma_{j}^{2}-\Gamma_{k}\sum_{j\neq
k}^{N}e^{-\mathrm{i}\phi_{j}}\Gamma_{j}}{N\bar{\Gamma}^{2}},$ (43)

provided that $\omega_{s}=0$. Assume further that the atoms are all equally
coupled with the cavity, i.e., $\Gamma_{1}=\Gamma_{2}=\cdots=\Gamma_{N}$, then
the steady state (42) reduces to

$|\Psi^{\prime}(\infty)\rangle=|D_{N}0\rangle\otimes\ket{\Phi_{0}}.$ (44)

In other words, in the steady state, the single excitation only exists in the
two-level atoms which have the same excitation probability $\frac{1}{N}$. The
Tavis-Cummings system cannot emit a photon into the cavity or the output
field. This theoretical result is consistent with the experimental results
given in [6, Fig. 2(c), Fig. 4(a)], where small fluctuations of the collective
swapping dynamics are due to the inhomogeneity of the coupling strengths. In
fact, when $\Gamma_{1}=\Gamma_{2}=\cdots=\Gamma_{N}$ and $\omega_{s}=0$, by
Corollary III.1, $\hat{A}_{\bar{c}\bar{o}}=0$, and thus the subsystem
$(\hat{A}_{\bar{c}\bar{o}},\hat{B}_{\bar{c}\bar{o}},\hat{C}_{\bar{c}\bar{o}})$
is static. On the other hand, it is easy to show that
$\ket{\Psi(t)}\equiv\ket{\Psi(0)}=|D_{N}0\rangle\otimes\ket{\Phi_{0}}$ for all
$t\geq 0$.

Let $\alpha=1$, $\beta=0$, i.e., the Tavis-Cummings model is initialized in
the pure state $|B_{N}0\rangle$. In this case, the steady state is

$\displaystyle|\Psi^{\prime}(\infty)\rangle=$
$\displaystyle\frac{1}{\sqrt{N}}\sum_{k=1}^{N}\bigg{(}c^{\prime}_{k}(\infty)|g_{1}\cdots
e_{k}\cdots g_{N}0\Phi_{0}\rangle$ (45)
$\displaystyle+c_{N+1,k}(\infty)|g_{1}g_{2}\cdots
g_{N}01_{\varphi}\rangle\bigg{)},$

where the pulse shape $\varphi$ is given in (36), and
$c^{\prime}_{k}(\infty)=\frac{\sum_{j\neq
k}^{N}\Gamma_{j}^{2}-\Gamma_{k}\sum_{j\neq
k}^{N}\Gamma_{j}}{N\bar{\Gamma}^{2}}$, provided that $\omega_{s}=0$, and
$c_{N+1,k}(\infty)$ is given in (38). Assume further that the atoms are all
equally coupled with the cavity, i.e.,
$\Gamma_{1}=\Gamma_{2}=\cdots=\Gamma_{N}$, the steady state (45) reduces to

$|\Psi^{\prime}(\infty)\rangle=|g_{1}g_{2}\cdots g_{N}01_{\varphi}\rangle.$
(46)

In other words, due to the existence of cavity decay rate $\kappa$, the Tavis-
Cummings system eventually emits a photon into the field.

###### Remark IV.2

The superradiant state $\ket{B_{N}}$ and subradiant state $\ket{D_{N}}$ can be
represented by the eigenvectors given in (22), respectively. For simplicity,
we assume the coupling strengths $\Gamma_{1}=\cdots=\Gamma_{N}\equiv\Gamma$.
If the state $|g_{1}\cdots e_{k}\cdots g_{N}0\rangle$ is viewed as a column
vector $[0\cdots 1\cdots 0]^{\top}$ (the $k$th element is $1$ and the others
are $0$), then $|B_{N}0\rangle$ is exactly the vector $T_{N}$ in (22).
Moreover, $A|B_{N}0\rangle=-\mathrm{i}\omega_{s}T_{N}-\mathrm{i}\sqrt{N}\Gamma
T_{N+1}$, which lives in the linear span of the two vectors $T_{N}$ and
$T_{N+1}$ of the controllable/observable subsystem $\hat{A}_{co}$ given in
(24). Hence, the single excitation swaps among all the two-level atoms and the
cavity. This echos the collective swapping dynamics shown in [6, Fig. 2(a)].
Moreover, as shown in (46), eventually all the oscillations will die out as
the photon is emitted into the external field due to the lossy nature of the
cavity. On the other hand, the subradiant state $|D_{N}0\rangle$ can be
represented by

$|D_{N}0\rangle=\sum_{j=1}^{N-1}\alpha_{j}T_{j},$ (47)

where
$\alpha_{j}=-\sqrt{\frac{j+1}{jN}}\left[e^{-\mathrm{i}\phi_{j+1}}+\frac{1}{j+1}\sum_{k=j+2}^{N}e^{-\mathrm{i}\phi_{k}}\right]$,
for $j=1,2,\ldots,N-2$, and
$\alpha_{N-1}=-\frac{1}{\sqrt{N-1}}e^{-\mathrm{i}\phi_{N}}$. Therefore, the
subradiant state given in (47) is a linear combination of the $N-1$
eigenvectors of the DFS given in (24). In this case, the Tavis-Cummings model
is neither reachable by its input nor detectable by its output. Thus, the two-
level atoms are decoupled from the cavity mode and the external field, and the
Tavis-Cummings model initialized in $|D_{N}0\rangle$ cannot emit a photon into
the cavity or even the external field, and remains in the subradiant state
(44).

##  V The multi-excitation case

In Section IV, we studied the single-excitation dynamics of the Tavis-Cummings
model. In this section, we present numerical studies of the multi-excitation
scenario. For ease of representation, we calculate the excitation
probabilities of the first two-level atom.

###  V-A The reduced density matrix of the first two-level atom

In this subsection, we present the master equation for the Tavis-Cummings
model; we also give the expression of the reduced density of the first two-
level atom.

According to [23], the master equation of the Tavis-Cummings model driven by a
single-photon input $\ket{1_{\xi}}$ is

$\displaystyle\dot{\bar{\rho}}^{11}(t)$
$\displaystyle=\mathcal{L}^{\star}\bar{\rho}^{11}(t)+\xi(t)[\bar{\rho}^{01}(t),L^{\ast}]+\xi^{\ast}(t)[L,\bar{\rho}^{10}(t)],$
(48) $\displaystyle\dot{\bar{\rho}}^{10}(t)$
$\displaystyle=\mathcal{L}^{\star}\bar{\rho}^{10}(t)+\xi(t)[\bar{\rho}^{00}(t),L^{\ast}],$
$\displaystyle\dot{\bar{\rho}}^{01}(t)$
$\displaystyle=\mathcal{L}^{\star}\bar{\rho}^{01}(t)+\xi^{\ast}(t)[L,\bar{\rho}^{00}(t)],$
$\displaystyle\dot{\bar{\rho}}^{00}(t)$
$\displaystyle=\mathcal{L}^{\star}\bar{\rho}^{00}(t),\ t\geq t_{0},$

where the initial states are

$\bar{\rho}^{11}(t_{0})=\bar{\rho}^{00}(t_{0})=|\zeta\rangle\langle\zeta|,~{}\bar{\rho}^{10}(t_{0})=\bar{\rho}^{01}(t_{0})=0$
(49)

with $\ket{\zeta}$ being the initial system state. As we will focus on the
excitation probability of the first two-level atom, we use the partial trace
to get its reduced density operator

$\displaystyle\rho_{A_{1}}(t)\triangleq$ $\displaystyle\
\mathrm{Tr}_{A_{N}}[\cdots\mathrm{Tr}_{A_{2}}[\mathrm{Tr}_{\mathrm{cav}}[\bar{\rho}^{11}(t)]]]$
(50) $\displaystyle=$ $\displaystyle\sum\langle z_{2}z_{3}\ldots
z_{N}n|\bar{\rho}^{11}(t)|z_{2}z_{3}\ldots z_{N}n\rangle,$

where $|z_{j}\rangle=|g_{j}\rangle$ or $|e_{j}\rangle$, $j=2,3,\ldots,N$, and
$\ket{n}$ is the state of the cavity.

###### Remark V.1

When the input is in the vacuum state $\ket{\Phi_{0}}$, (48)-(49) reduce to
the commonly used master equation

$\dot{\bar{\rho}}^{00}(t)=\mathcal{L}^{\star}\bar{\rho}^{00}(t),\ \ \
\bar{\rho}^{00}(t_{0})=\ket{\zeta}\bra{\zeta},\ t\geq t_{0}.$ (51)

Accordingly, the excitation probability of the first two-level atom is given
by $\braket{e_{1}}{\rho_{A_{1}}(t)}{e_{1}}$, where $\rho_{A_{1}}(t)$ should be
computed via (50) by replacing $\bar{\rho}^{11}(t)$ with $\bar{\rho}^{00}(t)$.

In the following simulations, we assume that there are three resonant two-
level atoms that are equally coupled to the cavity, i.e.,
$\omega_{1}=\omega_{2}=\omega_{3}=\omega_{r}=1$ and
$\Gamma_{1}=\Gamma_{2}=\Gamma_{3}=1$. Also let $\kappa=1.5$. When the input
field is in the single photon state $\ket{\Phi_{1}}$, we use a Gaussian pulse
shape

$\xi(t)=\left(\frac{\Omega^{2}}{2\pi}\right)^{1/4}\exp\left[-\frac{\Omega^{2}}{4}(t-t_{p})^{2}\right],$
(52)

where $t_{p}$ is the peak arrival time of the photon, and $\Omega$ is the
frequency bandwidth. Set $t_{p}=3$ and $\Omega=2\kappa$. See the green solid
curve in Fig. 3 for the plot of $|\xi(t)|^{2}$. As $|\xi(t)|^{2}\approx 0$
when $t=0$, it is safe to let the initial time $t_{0}$ be 0.

###  V-B One initially excited atom plus a single-photon input

Assume that one of the two-level atoms is in the excited state, the cavity is
empty, and the system is driven by a single-photon input. Then there are two
excitations in the whole system. According to the reduced density matrix (50),
the excitation probability of the first two-level atom in this two-excitation
case can be calculated as

$\displaystyle P_{\mathrm{TLS}1}(t)=\braket{e_{1}}{\rho_{A_{1}}(t)}{e_{1}}$
$\displaystyle=$ $\displaystyle\langle
e_{1}g_{2}g_{3}0|\bar{\rho}^{11}(t)|e_{1}g_{2}g_{3}0\rangle+\langle
e_{1}e_{2}g_{3}0|\bar{\rho}^{11}(t)|e_{1}e_{2}g_{3}0\rangle$
$\displaystyle+\langle
e_{1}g_{2}e_{3}0|\bar{\rho}^{11}(t)|e_{1}g_{2}e_{3}0\rangle+\langle
e_{1}g_{2}g_{3}1|\bar{\rho}^{11}(t)|e_{1}g_{2}g_{3}1\rangle,$

where the term $\langle
e_{1}g_{2}g_{3}0|\bar{\rho}^{11}(t)|e_{1}g_{2}g_{3}0\rangle$ indicates that a
photon is in the external field.

Figure 3: The excitation probabilities of the first two-level atom. Green
solid curve: $|\xi(t)|^{2}$; red solid curve (TLS1:
$\ket{e_{1}g_{2}g_{3}0}\otimes\ket{\Phi_{0}}$); blue solid curve (TLS2/3:
$\ket{g_{1}e_{2}g_{3}0}\otimes\ket{\Phi_{0}}$ or
$\ket{g_{1}g_{2}e_{3}0}\otimes\ket{\Phi_{0}}$); red dashed curve
(TLS1+$\ket{\Phi_{1}}$: $\ket{e_{1}g_{2}g_{3}0}\otimes\ket{\Phi_{1}}$); blue
dashed curve (TLS2/3+$\ket{\Phi_{1}}$:
$\ket{g_{1}e_{2}g_{3}0}\otimes\ket{\Phi_{1}}$ or
$\ket{g_{1}g_{2}e_{3}0}\otimes\ket{\Phi_{1}}$).

We consider the effect of single-photon input state on the excitation
probability, which is simulated in Fig. 3. The red solid curve is for the case
when the initial joint system-field state
$\ket{\Psi(0)}=\ket{e_{1}g_{2}g_{3}0}\otimes\ket{\Phi_{0}}$, the blue solid
curve is for $\ket{\Psi(0)}=\ket{g_{1}e_{2}g_{3}0}\otimes\ket{\Phi_{0}}$ or
$\ket{\Psi(0)}=\ket{g_{1}g_{2}e_{3}0}\otimes\ket{\Phi_{0}}$, the red dashed
curve is for $\ket{\Psi(0)}=\ket{e_{1}g_{2}g_{3}0}\otimes\ket{\Phi_{1}}$, and
the blue dashed curve is for
$\ket{\Psi(0)}=\ket{g_{1}e_{2}g_{3}0}\otimes\ket{\Phi_{1}}$ or
$\ket{\Psi(0)}=\ket{g_{1}g_{2}e_{3}0}\otimes\ket{\Phi_{1}}$. Firstly, the
incident single photon hardly interacts with the first two-level atom when the
time $t$ is far earlier than the peak arrival time $t_{p}=3$, which is similar
to the case when a single photon of Gaussian pulse shape is used to excite a
two-level atom; see [36, Fig. 4(a)]. Secondly, the solid red curve is above
the solid blue cure for all time. However there is crossover between the
dashed curves. Thirdly, the input photon does not affect the steady-state
excitation probability of the first TLS (TLS1); however, it does affect the
transient dynamics. Finally, from the dashed curves it can be seen that the
incident photon tends to increase (or decrease) the excitation probability of
the first two-level atom when it is initialized in the ground (or excited)
state.

###  V-C Two initially excited atoms plus vacuum input

Figure 4: The excitation probabilities of the first two-level atom. Red solid
curve (TLS1: $\ket{e_{1}g_{2}g_{3}0}\otimes\ket{\Phi_{0}}$); blue solid curve
(TLS2/3: $\ket{g_{1}e_{2}g_{3}0}\otimes\ket{\Phi_{0}}$ or
$\ket{g_{1}g_{2}e_{3}0}\otimes\ket{\Phi_{0}}$); red dashed curve (TLS2 + TLS3:
$\ket{g_{1}e_{2}e_{3}0}\otimes\ket{\Phi_{0}}$); blue dashed curve (TLS1+
TLS2/3: $\ket{e_{1}e_{2}g_{3}0}\otimes\ket{\Phi_{0}}$ or
$\ket{e_{1}g_{2}e_{3}0}\otimes\ket{\Phi_{0}}$); green solid curve (TLS1+TLS2 +
TLS3:$\ket{e_{1}e_{2}e_{3}0}\otimes\ket{\Phi_{0}}$; green dashed curve
(TLS1+TLS2+TLS3: $\ket{e_{1}e_{2}e_{3}0}\otimes\ket{\Phi_{0}}$ and
$\Gamma_{1}\neq\Gamma_{2}\neq\Gamma_{3}$).

In the following simulations, we assume that the Tavis-Cummings system is
driven by the vacuum input $\ket{\Phi_{0}}$ and one, two, or even three atoms
are initially excited.

We present the following theoretical result first.

###### Theorem V.1

Suppose that the Tavis-Cummings system (6) is initialized in the state
$\ket{\zeta}=\ket{e_{1}\ldots e_{N}0}$ and driven by the vacuum input field.
If $\omega_{1}=\ldots=\omega_{N}$ and $\Gamma_{1}=\ldots=\Gamma_{N}$, then the
steady-state ($t=\infty$) output field must be in an $N$-photon state.

###### Proof:

We use the master equation (51) to prove this result. The steady-state system
state can be found by solving

$\mathcal{L}^{\star}\bar{\rho}^{00}(\infty)=0.$ (53)

Clearly, in the steady state ($t=\infty$), the cavity cannot contain any
photon; otherwise the photon leaks out from the lossy cavity ($\kappa\neq 0$).
In other words, such a state cannot be a steady state. Consequently, the
steady-state system state is of the form
$\bar{\rho}^{00}(\infty)=\rho_{A}\otimes\ket{0}\bra{0}$, where $\rho_{A}$ is
the steady state of the two-level atoms. By the form of the system Hamiltonian
$H_{\mathrm{TC}}$ and the coupling operator $L$ given in subsection III-A,
(53) is actually

$\displaystyle\frac{1}{2}\sum_{j}\omega_{j}[\sigma_{z,j},\rho_{A}]\otimes\ket{0}\bra{0}$
$\displaystyle+$
$\displaystyle\sum_{j}\Gamma_{j}\bigg{(}\sigma_{-,j}\rho_{A}\otimes\ket{1}\bra{0}-\rho_{A}\sigma_{+,j}\ket{0}\bra{1}\bigg{)}=0.$

Given $\omega_{1}=\ldots=\omega_{N}\equiv\omega_{s}$ and
$\Gamma_{1}=\ldots=\Gamma_{N}\equiv\Gamma$, we have

$\displaystyle\frac{\omega_{s}}{2}\sum_{j}[\sigma_{z,j},\rho_{A}]\otimes\ket{0}\bra{0}$
$\displaystyle+\Gamma\sum_{j}\bigg{(}\sigma_{-,j}\rho_{A}\otimes\ket{1}\bra{0}-\rho_{A}\sigma_{+,j}\ket{0}\bra{1}\bigg{)}=0,$

which yields

$\sum_{j}\sigma_{-,j}\rho_{A}=0.$ (54)

In general, $\rho_{A}$ is of the form

$\sum_{i_{1},\ldots,i_{N};j_{1},\ldots,j_{N}}\alpha_{i_{1},\ldots,i_{N};j_{1},\ldots,j_{N}}\ket{f_{i_{1}}\ldots
f_{i_{N}}}\bra{g_{j_{1}}\ldots g_{j_{N}}},$

where $f_{i_{k}}$ and $g_{j_{k}}$ are either $e_{k}$ or $g_{k}$. Because
$\omega_{1}=\ldots=\omega_{N}$ and $\Gamma_{1}=\ldots=\Gamma_{N}$, all the
atoms are indistinguishable. As a result, all the coefficients must be
identical; i.e., $\alpha_{i_{1},\ldots,i_{N};j_{1},\ldots,j_{N}}\equiv\alpha$
for some $\alpha$. Consequently,

$\rho_{A}=\alpha\sum_{i_{1},\ldots,i_{N};j_{1},\ldots,j_{N}}\ket{f_{i_{1}}\ldots
f_{i_{N}}}\bra{g_{j_{1}}\ldots g_{j_{N}}}.$

Substituting this form of $\rho_{A}$ into (54) we get
$\sigma_{-,j}\ket{f_{i_{1}}\ldots f_{i_{N}}}=0,$ for all $j=1,\ldots,m$.
Clearly, all $f_{i_{k}}$ must be $g_{k}$. Consequently, $\rho_{A}$ contains a
single term and is of the form $\rho_{A}=\ket{g_{1}\ldots
g_{N}}\bra{g_{1}\ldots g_{N}}$. In other words, all atoms are in their ground
state. Therefore, in the steady state, the output field is in an $N$-photon
state. ∎

The simulation results are shown in Fig. 4. We have the following
observations. (i) The first atom’s excitation probability when it is
initialized in the excited state (the red solid curve) is greater than that
when either the second atom or the third atom is initialized in the excited
state (the blue solid curve). (ii) The first atom’s excitation probability
when both the second and third atoms are initialized in the excited state (the
red dashed curve) is greater than that when both the first and second atoms
(or both the first and third atoms) are initialized in the excited state (the
blue dashed curve). (iii) The red solid and dashed curves have the same final
value $\approx 0.44$, while the blue solid and dashed curves have the same
final value $\approx 0.11$. (iv) The excitation probability of the first two-
level atom settles to $0$ when all the three atoms are initially excited, see
the green solid curve in Fig. 4. Actually, by Theorem V.1, all the three atoms
eventually settle in their ground state and a 3-photon output state is
generated. However, if the coupling strengths are not identical, the atoms may
not settle to their ground state, as demonstrated by the green dashed curve in
Fig. 4, where $\Gamma_{1}=1$, $\Gamma_{2}=1.5$, and $\Gamma_{3}=2$.

###### Remark V.2

In Figs. 3 and 4, when the first atom is initialized in the excited state and
the other two are in the ground state, the final value of the excitation
probability of the first atom is approximately $0.44$. By Theorem IV.2,
$|c_{1}(\infty)|^{2}=\frac{4}{9}$, which explains this simulation result. The
other final value $0.11$ in Figs. 3 and 4 can be explained by
$|c_{2}(\infty)|^{2}=|c_{3}(\infty)=\frac{1}{9}$. Moreover, $|c_{k}(t)|^{2}\to
1$ and $|c_{j}(t)|^{2}\to 0$ ($j\neq k$) as $N\to\infty$. This means that the
system tends to remain intact when the number of atoms is sufficiently large.
Finally, if $\Gamma_{1}=\cdots=\Gamma_{N}$, then by (38) we have
$|c_{N+1,k}(\infty)|^{2}=\frac{1}{N}$, i.e., only $1/N$ of the excitation
energy is radiated. This confirms the analysis in the second paragraph above
(29) in [1].

###### Remark V.3

Although the red solid and dashed curves have identical stationary excitation
probability, their crests and troughs are almost symmetrical during the
transient process. Similar phenomena can be observed when $N=2$. This is
consistent with the experimental result given in [4, Fig. 4c]. Specifically,
write states $|e_{1}g_{2}\rangle$ and $|g_{1}e_{2}\rangle$ as
$\ket{\uparrow\downarrow}$ and $\ket{\downarrow\uparrow}$, respectively. In
[4, Fig. 4c], the solid red (green) circles plot the evolution of
$\ket{\downarrow\uparrow}$ ( $\ket{\uparrow\downarrow}$). The coherent state
transfer of the two states emerges and oscillates symmetrically.

##  VI A general form of the joint system-field states

In Section IV, analytical results are presented for the single-excitation
Tavis-Cummings model. These results are not applicable to the multi-excitation
case discussed in Section V. Motivated by this, in this section we derive a
recursive relation for the computation of the joint system-field state when
the system initially contains multiple excitations $R>1$ and is driven by $m$
input fields initialized in the vacuum state.

###  VI-A Modeling

Let $H$ be the internal system Hamiltonian of a quantum system. As the system
Hamiltonian may be tuned by a time-varying classical signal, for example,
transition frequencies of superconducting qubits can be tuned via Josephson
energy in Josephson-junction-based superconducting circuits [37], we use
$H(t)$ to emphasize the explicit dependence of the system Hamiltonian on time.
The coupling between the quantum system and the $k$th ($k=1,\ldots,m$) input
channel can be described by a coupling operator $L_{k}(t)$; again, here we
allow explicit time dependence of the coupling as in some quantum systems
couplings can be tuned in real time [38].

The system of interest may be an ensemble of TLSs, or resonators, or even an
ensemble of TLSs residing in a resonator like the Tavis-Cummings model studied
in this paper. If the number of TLSs is finite, the number of photons in the
resonator is finite, and the system is driven by a finite number of photons,
then the number of total excitations $R$ in the whole system (the system plus
the external field) is finite too. In this case, the state of the quantum
system has an orthonormal basis of the form

$\\{\ket{0_{s}},\ldots,\ket{{(K-1)}_{s}}\\},$ (55)

where the subscript “$s$” indicates that the basis states $\ket{j_{s}}$,
$j=0,1,\ldots,K-1$, are for the system. For our Tavis-Cummings mode in Fig. 1,
$K=2^{N}(R+1)$ as each of the $N$ atoms has two basis states and the cavity
can contain at most $R$ photons.

###### Remark VI.1

If the system of interest is coherently driven, then the number of excitations
can be arbitrary large as the drive may generate an arbitrary number of
excitations. However, if the coherent drive is not very strong and lasts not
long, the number of excitations in the system will be upper bounded. Hence, it
can still be assumed that the system admits a basis of finitely many elements.

For notational convenience, the following notation will be adopted. For each
$i=1,\ldots,m$, we use $t^{i}_{1\to k_{i}}$ to denote a set of ordered real
numbers $\\{t^{i}_{1},\ldots,t^{i}_{k_{i}}:t^{i}_{1}<\cdots<t^{i}_{k_{i}}\\}$.
Similarly, $1_{t^{j}_{1\to k_{j}}}$ is a shorthand of
$1_{t^{j}_{1}},\cdots,1_{t^{j}_{k_{j}}}$. Let $\int_{r}^{t^{i}_{2\to k_{i}}}$
be the abbreviation for $\int_{r}^{t^{i}_{k_{i}}}\cdots\int_{r}^{t^{i}_{2}}$.
Finally, for each given non-negative integer $n$,
$\sum_{k_{1},\ldots,k_{m}}^{n}$ means the summation over all possible
combinations of the non-negative integers $k_{1},\ldots,k_{m}$ that satisfy
$\sum_{j=1}^{m}k_{j}=n$.

The temporal photon-number basis of $n$ photons superposed over $m$ channels
is

$\displaystyle\Bigg{\\{}$
$\displaystyle|1_{t^{1}_{1}}\rangle\otimes\cdots\otimes|1_{t^{1}_{k_{1}}}\rangle\otimes\cdots\otimes|1_{t^{m}_{1}}\rangle\otimes\cdots\otimes|1_{t^{m}_{k_{m}}}\rangle:$
(56) $\displaystyle
t^{1}_{1}<\cdots<t^{1}_{k_{1}},\cdots,t^{m}_{1}<\cdots<t^{m}_{k_{m}},\sum_{i=1}^{m}k_{i}=n\Bigg{\\}},$

where $k_{1},\ldots,k_{m}$ are non-negative integers. Particularly, if
$k_{i}=0$, then there are no photons in channel $i$ and correspondingly the
term $|1_{t^{i}_{1}}\rangle\otimes\cdots\otimes|1_{t^{i}_{k_{i}}}\rangle$
reduces to $\ket{\Phi_{0}}$. As a common practice in quantum physics, the
tensor product state in (56) is often written as
$|1_{t^{1}_{1}},\ldots,1_{t^{1}_{k_{1}}},\ldots,1_{t^{m}_{1}},\ldots,1_{t^{m}_{k_{m}}}\rangle$,
which is $\ket{1_{t^{1}_{1\to k_{1}}},\ldots,1_{t^{m}_{1\to k_{m}}}}$ by means
of the notation given in the last paragraph.

Using the temporal photon-number basis in (56), an $m$-channel $n$-photon
state can be written as

$\displaystyle|n_{\xi}\rangle=\sum_{k_{1},\ldots,k_{m}}^{n}\int_{-\infty}^{\infty}\cdots\int_{-\infty}^{\infty}\xi(t^{1}_{1\to
k_{1}},\ldots,t^{m}_{1\to k_{m}})$ (57) $\displaystyle\ \ \ \ \ \ \ \ \ \ \ \
\ \ \ |1_{t^{1}_{1\to k_{1}}},\ldots,1_{t^{m}_{1\to k_{m}}}\rangle
dt^{1}_{1\to k_{1}}\cdots dt^{m}_{1\to k_{m}},$

where the $j$th channel has $k_{j}$ photons, and $\sum_{j=1}^{m}k_{j}=n$ is
the total photon number. Clearly, $\xi(t^{1}_{1\to k_{1}},\ldots,t^{m}_{1\to
k_{m}})$ is the probability amplitude of the component that the $j$th channel
has photons at time instants $t^{j}_{1},\ldots,t^{j}_{k_{j}}$ for all
$j=1,\ldots,m$.

Due to the indistinguishability of photons in each channel, for each fixed
$j=1,\ldots,m$, $\xi(t^{1}_{1\to k_{1}},\ldots,t^{m}_{1\to k_{m}})$ is
permutation-invariant with respect to indices $\\{t^{j}_{1\to k_{j}}\\}$.
Thus, under scaling $\frac{1}{k_{1}!\cdots k_{m}!}$, the state
$|n_{\xi}\rangle$ can be rewritten as

$\displaystyle|n_{\xi}\rangle$ $\displaystyle=$
$\displaystyle\sum_{k_{1},\ldots,k_{m}}^{n}\int_{-\infty}^{\infty}\int_{-\infty}^{t^{1}_{2\to
k_{1}}}\cdots\int_{-\infty}^{\infty}\int_{-\infty}^{t^{m}_{2\to
k_{m}}}\xi(t^{1}_{1\to k_{1}},\ldots,t^{m}_{1\to k_{m}})$ $\displaystyle\ \ \
\ \ \ \ \ \ \ |1_{t^{1}_{1\to k_{1}}},\ldots,1_{t^{m}_{1\to k_{m}}}\rangle
dt^{1}_{1\to k_{1}}\cdots dt^{m}_{1\to k_{m}}.$ (58)

In fact, one can always define an $n$-photon state using (58) with an
arbitrary multivariate function $\xi(t^{1}_{1\to k_{1}},\ldots,t^{m}_{1\to
k_{m}})$ (under normalization), then obtain the form in (57) by permutating
indices $\\{t^{j}_{1\to k_{j}}\\}$ for all $j=1,\ldots,m$. The form of
$n$-photon states (58) is more often to see in practice than that in (57). For
example, by applying a coherent drive to a superconducting qubit embedded in a
chiral waveguide, the qubit may generate photon states of the form (58).
Therefore, in what follows we use (58) to describe $n$-photon states.

###  VI-B The general form

In this subsection, we present a computational procedure that can be used to
compute the joint system-field state. Some other computational framework can
be found in e.g., [40, 41].

The basis state of the joint system-field state when the system is at level
$\ket{j_{s}}$ and channel $i$ has $k_{i}$ photons at time instants
$t^{i}_{1\to k_{i}}$ is $\ket{j_{s}}\otimes|1_{t^{1}_{1\to
k_{1}}},\ldots,1_{t^{m}_{1\to k_{m}}}\rangle$. Hence, the joint system-field
state is of the general form

$\displaystyle\ket{\Psi}=$
$\displaystyle\sum_{j=0}^{K-1}\sum_{n=0}^{\infty}\sum_{k_{1},\ldots,k_{m}}^{n}\int_{-\infty}^{\infty}\int_{-\infty}^{t^{1}_{2\to
k_{1}}}\cdots\int_{-\infty}^{\infty}\int_{-\infty}^{t^{m}_{2\to k_{m}}}$ (59)
$\displaystyle\hskip 17.22217pt\xi^{j_{s}}(t^{1}_{1\to
k_{1}},\ldots,t^{m}_{1\to k_{m}})\ket{j_{s}}$ $\displaystyle\hskip
17.22217ptdB_{\mathrm{in},1}^{\ast}(t^{1}_{1\to k_{1}})\cdots
dB_{\mathrm{in},m}^{\ast}(t^{m}_{1\to k_{m}})\ket{\Phi_{0}},$

where we write informally
$dB_{\mathrm{in},k}^{\ast}(t)=b_{\mathrm{in},k}^{\ast}(t)dt$ and
$dB_{\mathrm{in},i}^{\ast}(t^{i}_{1\to k_{i}})$ as the product
$dB^{\ast}_{\mathrm{in},i}(t^{i}_{1})\cdots
dB^{\ast}_{\mathrm{in},i}(t^{i}_{k_{i}})$. The normalization condition
$\braket{\Psi}{\Psi}=1$ gives

$\displaystyle\sum_{j=0}^{K-1}\sum_{n=0}^{\infty}\sum_{k_{1},\ldots,k_{m}}^{n}\int_{-\infty}^{\infty}\int_{-\infty}^{t^{1}_{2\to
k_{1}}}\cdots\int_{-\infty}^{\infty}\int_{-\infty}^{t^{m}_{2\to k_{m}}}$
$\displaystyle\hskip 17.22217pt|\xi^{j_{s}}(t^{1}_{1\to
k_{1}},\ldots,t^{m}_{1\to k_{m}})|^{2}dt^{1}_{1\to k_{1}}\cdots dt^{m}_{1\to
k_{m}}=1.$

Denote

$\ket{\eta(t^{1}_{1\to k_{1}},\ldots,t^{m}_{1\to
k_{m}})}=\sum_{j=0}^{K-1}\xi^{j_{s}}(t^{1}_{1\to k_{1}},\ldots,t^{m}_{1\to
k_{m}})\ket{j_{s}}.$

It is worthwhile to note that $\ket{\eta(t^{1}_{1\to k_{1}},\ldots,t^{m}_{1\to
k_{m}})}$ in general is not normalized. The state $\ket{\Psi}$ in (59) can be
re-written as

$\displaystyle\ket{\Psi}=$
$\displaystyle\sum_{n=0}^{\infty}\sum_{k_{1},\ldots,k_{m}}^{n}\int_{-\infty}^{\infty}\int_{-\infty}^{t^{1}_{2\to
k_{1}}}\cdots\int_{-\infty}^{\infty}\int_{-\infty}^{t^{m}_{2\to k_{m}}}$
$\displaystyle\hskip 17.22217pt\ket{\eta(t^{1}_{1\to k_{1}},\ldots,t^{m}_{1\to
k_{m}})}$ $\displaystyle\hskip 17.22217ptdB_{\mathrm{in},1}^{\ast}(t^{1}_{1\to
k_{1}})\cdots dB_{\mathrm{in},m}^{\ast}(t^{m}_{1\to k_{m}})\ket{\Phi_{0}}.$

If the input field is initially in the vacuum state $\ket{\Phi_{0}}$, then the
initial joint system-field state is
$\ket{\Psi(t_{0})}=\sum_{j=0}^{K-1}\xi^{j_{s}}\ket{j_{s}}\otimes\ket{\Phi_{0}}\equiv\ket{\eta_{0}}\otimes\ket{\Phi_{0}}$,
where $\sum_{j=0}^{K-1}|\xi^{j_{s}}|^{2}=1$. Because the photon generation
time is from the initial time $t_{0}=0$ to the present time $t>0$, the joint
system-field state at time $t$ is

$\displaystyle\ket{\Psi(t)}=$
$\displaystyle\sum_{n=0}^{\infty}\sum_{k_{1},\ldots,k_{m}}^{n}\int_{0}^{t}\int_{0}^{t^{1}_{2\to
k_{1}}}\cdots\int_{0}^{t}\int_{0}^{t^{m}_{2\to k_{m}}}$ (60)
$\displaystyle\ket{\eta_{t}(t^{1}_{1\to k_{1}},\ldots,t^{m}_{1\to k_{m}})}$
$\displaystyle dB_{\mathrm{in},1}^{\ast}(t^{1}_{1\to k_{1}})\cdots
dB_{\mathrm{in},m}^{\ast}(t^{m}_{1\to k_{m}})\ket{\Phi_{0}},$

where

$\ket{\eta_{t}(t^{1}_{1\to k_{1}},\ldots,t^{m}_{1\to
k_{m}})}=\sum_{j=0}^{K-1}\xi^{j_{s}}_{t}(t^{1}_{1\to k_{1}},\ldots,t^{m}_{1\to
k_{m}})\ket{j_{s}}.$ (61)

(Here, the subscript “$t$” indicates that the coefficients are time
dependent.) In particular, the term in $\ket{\Psi(t)}$ that corresponds to
$n=0$ is
$\ket{\eta_{t}}\ket{\Phi_{0}}=\sum_{j=0}^{K-1}\xi^{j_{s}}_{t}\ket{j_{s}}\ket{\Phi_{0}}$.
In other words, the field is in the vacuum state and the system is in the
superposition state $\sum_{j=0}^{K-1}\xi^{j_{s}}_{t}\ket{j_{s}}$.

In what follows, we aim to derive formulas for computing the joint system-
field state $\ket{\Psi(t)}$. Differentiating both sides of (60) yields

$\displaystyle d\ket{\Psi(t)}$ (62) $\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}\sum_{k_{1},\ldots,k_{m}}^{n}dB_{\mathrm{in},1}^{\ast}(t)\int_{0}^{t}\int_{0}^{t^{1}_{2\to
k_{1}-1}}\cdots\int_{0}^{t}\int_{0}^{t^{m}_{2\to k_{m}}}$ $\displaystyle\hskip
51.6665pt\ket{\eta_{t}(t^{1}_{1\to k_{1}-1},t,t^{m}_{2\to
k_{2}},\ldots,t^{m}_{1\to k_{m}})}$ $\displaystyle\hskip
51.6665ptdB_{\mathrm{in},1}^{\ast}(t^{1}_{1\to k_{1}-1})\cdots
dB_{\mathrm{in},m}^{\ast}(t^{m}_{1\to k_{m}})\ket{\Phi_{0}}$
$\displaystyle\vdots$
$\displaystyle+\sum_{n=1}^{\infty}\sum_{k_{1},\ldots,k_{m}}^{n}dB_{\mathrm{in},m}^{\ast}(t)\int_{0}^{t}\int_{0}^{t^{1}_{2\to
k_{1}}}\cdots\int_{0}^{t}\int_{0}^{t^{m}_{2\to k_{m}-1}}$ $\displaystyle\hskip
60.27759pt\ket{\eta_{t}(t^{1}_{1\to k_{1}},\ldots,t^{m-1}_{1\to
k_{m-1}},t^{m}_{1\to k_{m}-1},t)}$ $\displaystyle\hskip
60.27759ptdB_{\mathrm{in},1}^{\ast}(t^{1}_{1\to k_{1}})\cdots
dB_{\mathrm{in},m}^{\ast}(t^{m}_{1\to k_{m}-1})\ket{\Phi_{0}}$
$\displaystyle+dt\sum_{n=0}^{\infty}\sum_{k_{1},\ldots,k_{m}}^{n}\int_{-\infty}^{t}\int_{-\infty}^{t^{1}_{2\to
k_{1}}}\cdots\int_{-\infty}^{t}\int_{-\infty}^{t^{m}_{2\to k_{m}}}$
$\displaystyle\hskip 60.27759pt\ket{\dot{\eta}_{t}(t^{1}_{1\to
k_{1}},\ldots,t^{m}_{1\to k_{m}})}$ $\displaystyle\hskip
60.27759ptdB_{\mathrm{in},1}^{\ast}(t^{1}_{1\to k_{1}})\cdots
dB_{\mathrm{in},m}^{\ast}(t^{m}_{1\to k_{m}})\ket{\Phi_{0}},$

where $\dot{\eta}_{t}$ means the derivative of $\eta$ with respect to $t$.
Comparing the right-hand sides of (62) and (3) (for $\ket{\Psi(t)}$ in (60)),
we get

$\displaystyle\ket{\dot{\eta}_{t}(t^{1}_{1\to k_{1}},\ldots,t^{m}_{1\to
k_{m}})}$ (63) $\displaystyle=$
$\displaystyle-\mathrm{i}H_{\mathrm{eff}}(t)\ket{\eta_{t}(t^{1}_{1\to
k_{1}},\ldots,t^{m}_{1\to k_{m}})},$

and

$\displaystyle\ket{\eta_{t}(t^{1}_{1\to k_{1}},\ldots,t^{i}_{1\to
k_{i}},t,\ldots,t^{m}_{1\to k_{m}})}$ (64) $\displaystyle=$ $\displaystyle
L_{i}(t)\ket{\eta_{t}(t^{1}_{1\to k_{1}},\ldots,t^{i}_{1\to
k_{i}},\ldots,t^{m}_{1\to k_{m}})}$

for all $i=1,\ldots,m$. Under the basis (55), the effective Hamiltonian
$H_{\mathrm{eff}}$ of the system has a matrix representation. Hence, similar
to what has been done in [13], define the propagator $V(t)$ (a matrix
function) that solves a system of deterministic homogeneous ODEs
$\dot{V}(t)=-\mathrm{i}H_{\mathrm{eff}}(t)V(t)$ under the initial condition
$V(0)=I$. Then one can define the transition matrix

$G(t,\tau)\triangleq V(t)V(\tau)^{-1},~{}t,\tau\geq 0.$ (65)

Clearly, for the $0$ photon state case, (63) and (65) yield

$\ket{\eta_{t}}=G(t,0)\ket{\eta_{0}},~{}t\geq 0.$ (66)

As is well-known in linear systems theory, see e.g. [42, Chapter 4],
iteratively using the transition matrix $G(t,\tau)$ in (65) and (64) we have

$\displaystyle\ket{\eta_{t}(t^{1}_{1\to k_{1}},t^{2}_{1\to
k_{2}},\ldots,t^{m}_{1\to k_{m}})}$ $\displaystyle=$ $\displaystyle\
G(t,t^{1}_{k_{1}})\ket{\eta_{t^{1}_{k_{1}}}(t^{1}_{1\to
k_{1}-1},t^{1}_{k_{1}},t^{2}_{1\to k_{2}},\ldots,t^{m}_{1\to k_{m}})}$
$\displaystyle=$ $\displaystyle\
G(t,t^{1}_{k_{1}})L_{1}(t^{1}_{k_{1}})\ket{\eta_{t^{1}_{k_{1}}}(t^{1}_{1\to
k_{1}-1},t^{2}_{1\to k_{2}},\ldots,t^{m}_{1\to k_{m}})}$ $\displaystyle=$
$\displaystyle\
G(t,t^{1}_{k_{1}})L_{1}(t^{1}_{k_{1}})G(t^{1}_{k_{1}},t^{1}_{k_{1}-1})$
$\displaystyle\times\ket{\eta_{t^{1}_{k_{1}-1}}(t^{1}_{1\to
k_{1}-1},t^{2}_{1\to k_{2}},\ldots,t^{m}_{1\to k_{m}})}$ $\displaystyle\vdots$
$\displaystyle=$ $\displaystyle\ G(t,t^{1}_{k_{1}})L_{1}(t^{1}_{k_{1}})\cdots
G(t^{1}_{2},t^{1}_{1})L_{1}(t^{1}_{1})$ $\displaystyle\times
G(t^{1}_{1},t^{2}_{k_{2}})\ket{\eta_{t^{2}_{k_{2}}}(t^{2}_{1\to
k_{2}},\ldots,t^{m}_{1\to k_{m}})}$ $\displaystyle\vdots$ $\displaystyle=$
$\displaystyle\ G(t,t^{1}_{k_{1}})L_{1}(t^{1}_{k_{1}})\cdots
G(t^{1}_{2},t^{1}_{1})L_{1}(t^{1}_{1})G(t^{1}_{1},t^{2}_{k_{2}})L_{2}(t^{2}_{k_{2}})\times\cdots$
$\displaystyle\times
G(t^{2}_{2},t^{2}_{1})L_{2}(t^{2}_{1})G(t^{2}_{1},t^{3}_{k_{3}})\ket{\eta_{t^{3}_{k_{3}}}(t^{3}_{1\to
k_{3}},\ldots,t^{m}_{1\to k_{m}})}$ $\displaystyle\vdots$ $\displaystyle=$
$\displaystyle\ G(t,t^{1}_{k_{1}})L_{1}(t^{1}_{k_{1}})\cdots
G(t^{1}_{2},t^{1}_{1})L_{1}(t^{1}_{1})$ $\displaystyle\times
G(t^{1}_{1},t^{2}_{k_{2}})L_{2}(t^{2}_{k_{2}})\cdots
G(t^{2}_{2},t^{2}_{1})L_{2}(t^{2}_{1})\times\cdots$ $\displaystyle\times
G(t^{m-1}_{1},t^{m}_{k_{m}})L_{m}(t^{m}_{k_{m}})\cdots
G(t^{m}_{2},t^{m}_{1})L_{m}(t^{m}_{1})G(t^{m}_{1},0)\ket{\eta_{0}}.$

Consequently, we have the following general formula

$\displaystyle\ket{\eta_{t}(t^{1}_{1\to k_{1}},t^{2}_{1\to
k_{2}},\ldots,t^{m}_{1\to k_{m}})}$ (67) $\displaystyle=$ $\displaystyle
G(t,t^{1}_{k_{1}})L_{1}(t^{1}_{k_{1}})\cdots
G(t^{1}_{2},t^{1}_{1})L_{1}(t^{1}_{1})$ $\displaystyle\times
G(t^{1}_{1},t^{2}_{k_{2}})L_{2}(t^{2}_{k_{2}})\cdots
G(t^{2}_{2},t^{2}_{1})L_{2}(t^{2}_{1})\cdots$ $\displaystyle\times
G(t^{m-1}_{1},t^{m}_{k_{m}})L_{m}(t^{m}_{k_{m}})\cdots
G(t^{m}_{2},t^{m}_{1})L_{m}(t^{m}_{1})$ $\displaystyle\times
G(t^{m}_{1},0)\ket{\eta_{0}}.$

###### Remark VI.2

To better understand the recursive algorithm given in (67), we take the two-
channel case as an example, i.e., $m=2$. In this case, the first channel
contains $k_{1}$ photons and the second channel contains $k_{2}$ photons. The
recursive functions (63) and (64) are reduced to be

$\displaystyle\ket{\dot{\eta}_{t}^{k_{1},k_{2}}(t^{1}_{1\to k_{1}},t^{2}_{1\to
k_{2}})}=-\mathrm{i}H_{\mathrm{eff}}(t)\ket{\eta_{t}^{k_{1},k_{2}}(t^{1}_{1\to
k_{1}},t^{2}_{1\to k_{2}})},$ (68)

and

$\displaystyle\ket{\eta_{t}^{k_{1}+1,k_{2}}(t^{1}_{1\to k_{1}},t,t^{2}_{1\to
k_{2}})}=L_{1}(t)\ket{\eta_{t}^{k_{1},k_{2}}(t^{1}_{1\to k_{1}},t^{2}_{1\to
k_{2}})},$ (69) $\displaystyle\ket{\eta_{t}^{k_{1},k_{2}+1}(t^{1}_{1\to
k_{1}},t^{2}_{1\to k_{2}},t)}=L_{2}(t)\ket{\eta_{t}^{k_{1},k_{2}}(t^{1}_{1\to
k_{1}},t^{2}_{1\to k_{2}})}.$

According to (67), we have

$\displaystyle\ket{\eta_{t}^{k_{1},k_{2}}(t^{1}_{1\to k_{1}},t^{2}_{1\to
k_{2}})}$ (70) $\displaystyle=$ $\displaystyle
G(t,t^{1}_{k_{1}})L_{1}(t^{1}_{k_{1}})\cdots
G(t^{1}_{2},t^{1}_{1})L_{1}(t^{1}_{1})$ $\displaystyle\times
G(t^{1}_{1},t^{2}_{k_{2}})L_{2}(t^{2}_{k_{2}})\cdots
G(t^{2}_{2},t^{2}_{1})L_{2}(t^{2}_{1})$ $\displaystyle\times
G(t^{2}_{1},0)\ket{\eta_{0}},$

which yields

$\displaystyle\ket{\eta_{t}}=$ $\displaystyle G(t,0)\ket{\eta_{0}},$ (71)
$\displaystyle\ket{\eta_{t}^{1,0}(t^{1}_{1})}=$ $\displaystyle
G(t,t_{1}^{1})L_{1}(t^{1}_{1})G(t_{1}^{1},0)\ket{\eta_{0}},$
$\displaystyle\ket{\eta_{t}^{0,1}(t^{1}_{2})}=$ $\displaystyle
G(t,t^{2}_{1})L_{2}(t^{2}_{1})G(t^{2}_{1},0)\ket{\eta_{0}},$
$\displaystyle\ket{\eta_{t}^{1,1}(t^{1}_{1},t^{2}_{1})}=$ $\displaystyle
G(t,t^{1}_{1})L_{1}(t^{1}_{1})G(t^{1}_{1},t^{2}_{1})L_{2}(t^{2}_{1})G(t^{2}_{1},0)\ket{\eta_{0}}$
$\displaystyle=$ $\displaystyle
G(t,t^{2}_{1})L_{2}(t^{2}_{1})G(t^{2}_{1},t^{1}_{1})L_{1}(t^{1}_{1})G(t^{1}_{1},0)\ket{\eta_{0}},$
$\displaystyle\vdots$

That is, there are multiple sequential orders to generate a
($k_{1}+k_{2}$)-photon state distributed in two channels ($k_{1},k_{2}\geq
1$), which allow the two channels being created crosswise.

The following recursive relation turns out useful for deriving all the states.

$\displaystyle\ket{\eta_{t}(t^{1}_{1\to k_{1}},\ldots,t^{i}_{1\to
k_{i}+1},\ldots,t^{m}_{1\to k_{m}})}$ (72) $\displaystyle=$ $\displaystyle\
G(t,t^{i}_{k_{i}+1})L_{i}(t^{i}_{k_{i}+1})G(t,t^{i}_{k_{i}+1})^{-1}$
$\displaystyle\ket{\eta_{t}(t^{1}_{1\to k_{1}},\ldots,t^{i}_{1\to
k_{i}},\ldots,t^{m}_{1\to k_{m}})}.$

###### Remark VI.3

Theorem IV.2 in Section IV-C can be proved by applying the recursive relation
(72) to the special case when $m=1$ and the initial system-field state is

$\ket{\Psi(0)}=\ket{\Psi_{k}(0)}=|\eta_{0}\rangle\otimes\ket{\Phi_{0}},$ (73)

where $|\eta_{0}\rangle=|g_{1}\cdots e_{k}\cdots g_{N}0\rangle$. In what
follows we sketch the proof. Because there is only one input channel and the
number of excitation is 1, i.e., $m=1$ and $k_{1}=1$, the joint system-field
state $\ket{\Psi(t)}$ in (60) is

$\ket{\Psi(t)}=|\eta_{t}\rangle\otimes\ket{\Phi_{0}}+\int_{0}^{t}|\eta_{t}(t_{1})\rangle
dB_{\mathrm{in}}^{\ast}(t_{1})\otimes\ket{\Phi_{0}}.$ (74)

The first term in (74) indicates that the single excitation exists in the
Tavis-Cummings model. $|\eta_{t}\rangle$ can be described by

$|\eta_{t}\rangle=\sum_{j=1}^{N}c_{j}(t)|g_{1}\cdots e_{j}\cdots
g_{N}0\rangle+c_{N+2,k}(t)|g_{1}g_{2}\cdots g_{N}1\rangle$

with the initial conditions $c_{k}(0)=1$ and $c_{j}(0)=c_{N+2,k}(0)=0$ for
$1\leq j\leq N,~{}j\neq k$. On the other hand, the second term in (74) means
that the Tavis-Cummings model emits a photon into the output field, which can
be rewritten as

$\displaystyle\int_{0}^{t}|\eta_{t}(t_{1})\rangle
dB_{\mathrm{in}}^{\ast}(t_{1})\otimes\ket{\Phi_{0}}$ $\displaystyle=$
$\displaystyle
c_{N+1}(t)\int_{0}^{t}\varphi(t_{1})dB_{\mathrm{in}}^{\ast}(t_{1})|g_{1}g_{2}\cdots
g_{N}0\rangle\otimes\ket{\Phi_{0}}.$

Here, $\varphi(t_{1})$ is the pulse shape of the single-photon output state.
When $t\rightarrow\infty$, $c_{N+1}(\infty)$ denotes the probability amplitude
corresponding to the reduced joint state $|g_{1}g_{2}\cdots
g_{N}0\rangle\otimes|\Phi_{1}\rangle$. By the recursive relation (72), we have
$\ket{\eta_{t}(t_{1})}=G(t,t_{1})L_{1}(t_{1})G(t,t_{1})^{-1}\ket{\eta_{t}}$,
where $\ket{\eta_{t}}$ can be calculated via (66) and (73). Theorem IV.2
follows.

In the following subsections, we apply the above theory to our Tavis-Cummings
model. In this case, $m=1$, $L_{i}(t)$ in (64) is $L=\sqrt{\kappa}a$, and
$H_{\mathrm{eff}}(t)$ in (63) is
$H_{\mathrm{TC}}-\frac{\mathrm{i}}{2}L^{\ast}L$, where $H_{\mathrm{TC}}$ is
given in (5). For simplicity, we set
$\omega_{1}=\ldots=\omega_{N}\equiv\omega_{r}=0$ and
$\Gamma_{1}=\ldots=\Gamma_{N}=\kappa=1$. Assume that the Tavis-Cummings model
contains $R$ initially excited two-level atoms, the cavity is initially empty,
and the input is in the vacuum state. In such setting, the number of basis
states of the Tavis-Cummings model is $K=2^{N}(R+1)$ as discussed in
subsection VI-A. Finally, for better illustration we plot the symmetric pulse
shape in (57), instead of that in (58).

###  VI-C $N=3$ and $R=1$.

In this case, the total number of basis states is $K=16$. The joint system-
field state can be expressed as
$\ket{\Psi(t)}=|\eta_{t}\rangle\otimes\ket{\Phi_{0}}+\int_{0}^{t}|\eta_{t}(t_{1})\rangle
dB_{\mathrm{in}}^{\ast}(t_{1})\ket{\Phi_{0}}$, where by (61),
$|\eta_{t}(t_{1},\ldots,t_{k})\rangle=\sum_{j=0}^{15}\xi^{j_{s}}_{t}(t_{1},\ldots,t_{k})\ket{j_{s}},~{}k=0,1$.
Let the initial joint system-field state be
$\ket{\Psi(0)}=|e_{1}g_{2}g_{3}0\rangle\otimes\ket{\Phi_{0}}$. By the
recursive relation (72), we have

$\displaystyle\xi^{1_{s}}_{t}=\frac{2}{\sqrt{47}}\left(e^{\frac{-1-\sqrt{47}\mathrm{i}}{4}t}-e^{\frac{-1+\sqrt{47}\mathrm{i}}{4}t}\right),$
$\displaystyle\xi^{2_{s}}_{t}=-\frac{1}{3}+\frac{2\left(\frac{\sqrt{47}+\mathrm{i}}{4}e^{\frac{-1-\sqrt{47}\mathrm{i}}{4}t}+\frac{\sqrt{47}-\mathrm{i}}{4}e^{\frac{-1+\sqrt{47}\mathrm{i}}{4}t}\right)}{3\sqrt{47}},$
$\displaystyle\xi^{4_{s}}_{t}=\xi^{2_{s}}_{t},$
$\displaystyle\xi^{8_{s}}_{t}=\frac{2}{3}+\frac{2\left(\frac{\sqrt{47}+\mathrm{i}}{4}e^{\frac{-1-\sqrt{47}\mathrm{i}}{4}t}+\frac{\sqrt{47}-\mathrm{i}}{4}e^{\frac{-1+\sqrt{47}\mathrm{i}}{4}t}\right)}{3\sqrt{47}},$
$\displaystyle\xi^{0_{s}}_{t}(t_{1})=-\frac{4\mathrm{i}}{\sqrt{47}}e^{-\frac{1}{4}t_{1}}\sin\left(\frac{\sqrt{47}}{4}t_{1}\right),~{}0\leq
t_{1}\leq t,$

and all the others are 0. In the steady state ($t=\infty$), we have
$\xi^{2_{s}}_{\infty}=\xi^{4_{s}}_{\infty}=-\frac{1}{3},\xi^{8_{s}}_{\infty}=\frac{2}{3}$,
which correspond to the basis vectors
$|g_{1}g_{2}e_{3}0\rangle,|g_{1}e_{2}g_{3}0\rangle,|e_{1}g_{2}g_{3}0\rangle$
respectively. Moreover, as
$\lim_{t\to\infty}\int_{0}^{\infty}|\xi^{0_{s}}_{t}(t_{1})|^{2}dt_{1}=\frac{1}{3}$,
the steady-state joint system-field state is

$\displaystyle|\Psi(\infty)\rangle=\frac{2}{3}|e_{1}g_{2}g_{3}0\rangle\otimes\ket{\Phi_{0}}-\frac{1}{3}|g_{1}e_{2}g_{3}0\rangle\otimes\ket{\Phi_{0}}$
(75) $\displaystyle\hskip
34.44434pt-\frac{1}{3}|g_{1}g_{2}e_{3}0\rangle\otimes\ket{\Phi_{0}}+\frac{1}{\sqrt{3}}|g_{1}g_{2}g_{3}0\rangle\otimes|1_{\eta}\rangle,$

where the pulse shape of the single photon is
$\eta(t)=\sqrt{3}\xi^{0_{s}}_{\infty}(t)$ up to a global phase. Finally, the
square of the probability amplitude $2/3$ in (75) is consistent with the
limiting value of the red solid curve in Fig. 4.

###  VI-D $N=2$ and $R=2$.

In this case, the total number of basis states is $K=12$. The joint system-
field state can be expressed as

$\displaystyle\ket{\Psi(t)}=$
$\displaystyle|\eta_{t}\rangle\otimes\ket{\Phi_{0}}+\int_{0}^{t}|\eta_{t}(t_{1})\rangle
dB_{\mathrm{in}}^{\ast}(t_{1})\ket{\Phi_{0}}$
$\displaystyle+\int_{0}^{t}\int_{0}^{t_{2}}|\eta_{t}(t_{1},t_{2})\rangle
dB_{\mathrm{in}}^{\ast}(t_{1})dB_{\mathrm{in}}^{\ast}(t_{2})\ket{\Phi_{0}},$

where
$|\eta_{t}(t_{1},\ldots,t_{k})\rangle=\sum_{j=0}^{11}\xi^{j_{s}}_{t}(t_{1},\ldots,t_{k})\ket{j_{s}},~{}k=0,1,2$.
Let the initial joint system-field state be
$\ket{\Psi(0)}=|e_{1}e_{2}0\rangle\otimes\ket{\Phi_{0}}$. By Theorem V.1, the
steady-state joint system-field state ($t\rightarrow\infty$) is

$|\Psi(\infty)\rangle=|g_{1}g_{2}0\rangle\otimes|2_{\eta}\rangle,$ (76)

where $|2_{\eta}\rangle$ is the 2-photon output state. By the recursive
relation (72), the pulse shape of the 2-photon state $|2_{\eta}\rangle$ is

$\displaystyle\eta(t_{1},t_{2})=$
$\displaystyle\xi^{0_{s}}_{\infty}(t_{1},t_{2})$ $\displaystyle=$
$\displaystyle(0.0228787-0.00785441\mathrm{i})e^{\mu_{1}t_{1}-\mu_{4}t_{2}}$
$\displaystyle+$
$\displaystyle(0.319599+0.0291387\mathrm{i})e^{\mu_{2}t_{1}-\mu_{4}t_{2}}$
$\displaystyle-$
$\displaystyle(0.342477+0.0212843\mathrm{i})e^{\mu_{3}t_{1}-\mu_{4}t_{2}}+\mathrm{c.c.},$

where “$\mathrm{c.c.}$” means complex conjugate, and
$\mu_{1}=-0.336506+3.79453\mathrm{i}$, $\mu_{2}=-0.336506-1.01065\mathrm{i}$,
$\mu_{3}=-0.076987+1.39194\mathrm{i}$, and $\mu_{4}=0.25+1.39194\mathrm{i}$.
It can be easily verified that
$\int_{0}^{\infty}\int_{0}^{t_{2}}|\eta(t_{1},t_{2})|^{2}dt_{1}dt_{2}=1$. The
probability distribution of the steady-state output two-photon state
$|2_{\eta}\rangle$ in (76) is simulated in Fig. 5.

Figure 5: Probability distribution of $|2_{\eta}\rangle$ in (76).

###  VI-E $N=3$ and $R=2$.

In this case, the total number of basis states is $K=24$. The joint system-
field state can be expressed as

$\displaystyle\ket{\Psi(t)}=$
$\displaystyle|\eta_{t}\rangle\otimes\ket{\Phi_{0}}+\int_{0}^{t}|\eta_{t}(t_{1})\rangle
dB_{\mathrm{in}}^{\ast}(t_{1})\ket{\Phi_{0}}$
$\displaystyle+\int_{0}^{t}\int_{0}^{t_{2}}|\eta_{t}(t_{1},t_{2})\rangle
dB_{\mathrm{in}}^{\ast}(t_{1})dB_{\mathrm{in}}^{\ast}(t_{2})\ket{\Phi_{0}},$

where
$|\eta_{t}(t_{1},\ldots,t_{k})\rangle=\sum_{j=0}^{23}\xi^{j_{s}}_{t}(t_{1},\ldots,t_{k})\ket{j_{s}},~{}k=0,1,2$.
Let the initial joint system-field state be
$\ket{\Psi(0)}=|g_{1}e_{2}e_{3}0\rangle\otimes\ket{\Phi_{0}}$. According to
the recursive relation (72), we have the steady-state joint system-field state

$\displaystyle|\Psi(\infty)\rangle=$
$\displaystyle\frac{2}{3}|e_{1}g_{2}g_{3}0\rangle\otimes|1_{\eta^{{12}_{s}}}\rangle+\frac{1}{3}|g_{1}e_{2}g_{3}0\rangle\otimes|1_{\eta^{{6}_{s}}}\rangle$
(77)
$\displaystyle+\frac{1}{3}|g_{1}g_{2}e_{3}0\rangle\otimes|1_{\eta^{{3}_{s}}}\rangle+\frac{\sqrt{3}}{3}|g_{1}g_{2}g_{3}0\rangle\otimes|2_{\eta^{{0}_{s}}}\rangle,$

where
$\eta^{{12}_{s}}(t)=\frac{3}{2}\xi^{12_{s}}_{\infty}(t)=0.5163975\left(e^{\lambda^{\prime}t}-e^{\lambda^{\prime\ast}t}\right)$
with $\lambda^{\prime}=-0.25+0.968246\mathrm{i}$,
$\eta^{{6}_{s}}(t)=\eta^{{3}_{s}}(t)=-\eta^{{12}_{s}}(t)$, and

$\displaystyle\eta^{{0}_{s}}(t_{1},t_{2})=$
$\displaystyle\sqrt{3}\xi^{{0}_{s}}_{\infty}(t_{1},t_{2})$ $\displaystyle=$
$\displaystyle(0.338571+0.0140041\mathrm{i})e^{\lambda_{1}^{\prime}t_{1}-\lambda_{4}^{\prime}t_{2}}$
$\displaystyle+(0.0155773-0.00341903\mathrm{i})e^{\lambda_{2}^{\prime}t_{1}-\lambda_{4}^{\prime}t_{2}}$
$\displaystyle-(0.354149+0.0105851\mathrm{i})e^{\lambda_{3}^{\prime}t_{1}-\lambda_{4}^{\prime}t_{2}}+\mathrm{c.c.}$

with $\lambda_{1}^{\prime}=-0.301227-1.40985\mathrm{i}$,
$\lambda_{2}^{\prime}=-0.301227+4.83767\mathrm{i}$,
$\lambda_{3}^{\prime}=-0.147546+1.71391\mathrm{i}$, and
$\lambda_{4}^{\prime}=0.25+1.71391\mathrm{i}$. The square of the probability
amplitude $\frac{2}{3}$ in (77) is consistent with the limiting value of the
red dashed curve in Fig. 4. It is clear in Fig. 6 that $\ket{1_{\eta}}$ has
more oscillations than $\ket{1_{\eta^{{12}_{s}}}}$.

Figure 6: The comparison between probability distributions of the steady-state
single-photon output states $\ket{1_{\eta}}$ in (75) and
$\ket{1_{\eta^{{12}_{s}}}}$ in (77). Figure 7: Probability distribution of the
two-photon output state $|2_{\eta^{{0}_{s}}}\rangle$ in (77).

Both Figs. 5 and 7 show the probability distributions of a two-photon output
state. The one generated by the Tavis-Cummings system with $3$ atoms (shown in
Fig. 7) oscillates more rapidly than that with $2$ atoms (see Fig. 5). This
result indicates that the Rabi oscillation is enhanced by adding more atoms,
which is consistent with the discussions in Fig. 2 and the observations in
[43, Fig. 5].

###  VI-F $N=3$ and $R=3$.

In this case, the total number of basis states is $K=32$ and the joint system-
field state can be expressed as

$\displaystyle\ket{\Psi(t)}=|\eta_{t}\rangle\otimes\ket{\Phi_{0}}+\int_{0}^{t}|\eta_{t}(t_{1})\rangle
dB_{\mathrm{in}}^{\ast}(t_{1})\ket{\Phi_{0}}$ $\displaystyle+$
$\displaystyle\int_{0}^{t}\int_{0}^{t_{2}}|\eta_{t}(t_{1},t_{2})\rangle
dB_{\mathrm{in}}^{\ast}(t_{1})dB_{\mathrm{in}}^{\ast}(t_{2})\ket{\Phi_{0}}$
$\displaystyle+$
$\displaystyle\int_{0}^{t}\int_{0}^{t_{3}}\int_{0}^{t_{2}}|\eta_{t}(t_{1},t_{2},t_{3})\rangle
dB_{\mathrm{in}}^{\ast}(t_{1})dB_{\mathrm{in}}^{\ast}(t_{2})dB_{\mathrm{in}}^{\ast}(t_{3})\ket{\Phi_{0}},$

where
$|\eta_{t}(t_{1},\ldots,t_{k})\rangle=\sum_{j=0}^{31}\xi^{j_{s}}_{t}(t_{1},\ldots,t_{k})\ket{j_{s}},~{}k=0,1,2,3$.
Since the three atoms are all initially excited, the initial joint system-
field state is $\ket{\Psi(0)}=|e_{1}e_{2}e_{3}0\rangle\otimes\ket{\Phi_{0}}$.
By Theorem V.1, the steady-state joint system-field state is

$|\Psi(\infty)\rangle=|g_{1}g_{2}g_{3}0\rangle\otimes|3_{\eta}\rangle,$ (78)

where $|3_{\eta}\rangle$ is the three-photon output state, whose probability
distribution is shown in Fig. 8.

Figure 8: Probability distribution of the three-photon output state
$|3_{\eta}\rangle$ in (78).

##  VII Conclusion

In this paper, we have studied the Tavis-Cummings model. Specifically, we have
applied quantum linear systems theory to reveal typical features of the Tavis-
Cummings model, the analytical expression has been derived for the output
single-photon state of a Tavis-Cummings system in response to a single-photon
input. We have also proposed a computational framework to derive the analytic
form of the superposition state of the system and field. Note that the linear
quantum systems’ approach is only applicable when there is only one
excitation. Future studies are to be done on subradiance and superradiance in
the multi-excitation case.

## Appendix

As concepts and properties of quantum linear passive systems are used in this
paper, in particular Section IV, in this appendix we collect some of them for
readability of this paper.

For the quantum linear system (13), from (18) we know the transfer function
$G[s]$ is an all-pass function [44, pp. 357]. Actually, as shown in Remark
III.2, system (13) is passive. If we take expectation on both sides of (13)
with respect to the initial joint system-field state (which is a unit vector
in a Hilbert space), we get a classical linear system

$\displaystyle\left\langle\dot{\bar{a}}(t)\right\rangle$
$\displaystyle=A\left\langle\bar{a}(t)\right\rangle+B\left\langle
b_{\mathrm{in}}(t)\right\rangle,$ (79) $\displaystyle\left\langle
b_{\mathrm{out}}(t)\right\rangle$
$\displaystyle=C\left\langle\bar{a}(t)\right\rangle dt+\left\langle
b_{\mathrm{in}}(t)\right\rangle.$

Thus we can define controllability, observability, and Hurwitz stability for
the quantum linear system (13) using those for the classical linear system
(79).

###### Definition VII.1

The quantum linear passive system (13) is said to be Hurwitz stable (resp.
controllable, observable) if the corresponding classical linear system (79) is
Hurwitz stable (resp. controllable, observable).

The following result reveals the structure of quantum linear passive systems;
see [29, 31] for more details.

###### Proposition VII.1

The quantum linear passive system (13) has the following properties.

(i)

Its Hurwitz stability, controllability and observability are equivalent to
each other.

(ii)

It has only controllable and observable ($co$) subsystem and uncontrollable
and unobservable ($\bar{c}\bar{o}$) subsystem. Moreover, each $\bar{c}\bar{o}$
subsystem is a closed (namely isolated) quantum system; see Eq. (21) and
Corollary III.1.

(iii)

Its poles corresponding to a $\bar{c}\bar{o}$ subsystem are all on the
imaginary axis, while its poles corresponding to a $co$ subsystem are on the
open left half of the complex plane.

## References

  * [1] R. H. Dicke, “Coherence in spontaneous radiation processes,” _Physical Review_ , vol. 93, no. 1, p. 99, 1954.
  * [2] M. Tavis and F. W. Cummings, “Exact solution for an $n$-molecule—radiation-field hamiltonian,” _Physical Review_ , vol. 170, pp. 379–384, 1968.
  * [3] M. A. Sillanpää, J. I. Park, and R. W. Simmonds, “Coherent quantum state storage and transfer between two phase qubits via a resonant cavity,” _Nature_ , vol. 449, no. 7161, pp. 438–442, 2007.
  * [4] J. Majer, J. Chow, J. Gambetta, J. Koch, B. Johnson, J. Schreier, L. Frunzio, D. Schuster, A. A. Houck, A. Wallraff _et al._ , “Coupling superconducting qubits via a cavity bus,” _Nature_ , vol. 449, no. 7161, pp. 443–447, 2007.
  * [5] J. M. Fink, R. Bianchetti, M. Baur, M. Göppl, L. Steffen, S. Filipp, P. J. Leek, A. Blais, and A. Wallraff, “Dressed collective qubit states and the tavis-cummings model in circuit QED,” _Physical Review Letters_ , vol. 103, p. 083601, 2009.
  * [6] Z. Wang, H. Li, W. Feng, X. Song, C. Song, W. Liu, Q. Guo, X. Zhang, H. Dong, D. Zheng, H. Wang, and D.-W. Wang, “Controllable switching between superradiant and subradiant states in a 10-qubit superconducting circuit,” _Physical Review Letters_ , vol. 124, p. 013601, 2020.
  * [7] T. Astner, S. Nevlacsil, N. Peterschofsky, A. Angerer, S. Rotter, S. Putz, J. Schmiedmayer, and J. Majer, “Coherent coupling of remote spin ensembles via a cavity bus,” _Physical Review Letters_ , vol. 118, no. 14, p. 140502, 2017.
  * [8] G.-W. Deng, D. Wei, S.-X. Li, J. Johansson, W.-C. Kong, H.-O. Li, G. Cao, M. Xiao, G.-C. Guo, F. Nori _et al._ , “Coupling two distant double quantum dots with a microwave resonator,” _Nano Letters_ , vol. 15, no. 10, pp. 6620–6625, 2015.
  * [9] D. J. van Woerkom, P. Scarlino, J. H. Ungerer, C. Müller, J. V. Koski, A. J. Landig, C. Reichl, W. Wegscheider, T. Ihn, K. Ensslin, and A. Wallraff, “Microwave photon-mediated interactions between semiconductor qubits,” _Physical Review X_ , vol. 8, p. 041018, 2018.
  * [10] B. Wang, T. Lin, H. Li, S. Gu, M. Chen, G. Guo, H. Jiang, X. Hu, G. Cao, and G. Guo, “Correlated spectrum of distant semiconductor qubits coupled by microwave photons,” _Science Bulletin_ , vol. 66, no. 4, pp. 332–338, 2021\.
  * [11] M. Gross and S. Haroche, “Superradiance: An essay on the theory of collective spontaneous emission,” _Physics Reports_ , vol. 93, no. 5, pp. 301–396, 1982\.
  * [12] Q.-M. Chen, Y.-x. Liu, L. Sun, and R.-B. Wu, “Tuning the coupling between superconducting resonators with collective qubits,” _Physical Review A_ , vol. 98, p. 042328, 2018.
  * [13] W.-L. Li, G. Zhang, and R.-B. Wu, “The dynamical model of flying-qubit control systems,” in _20th IFAC World Congress_ , vol. 50, no. 1, 2020, pp. 11 755–11 759.
  * [14] Y. Pan and G. Zhang, “Scattering of few photons by a ladder-type quantum system,” _Journal of Physics A: Mathematical and Theoretical_ , vol. 50, no. 34, p. 345301, 2017.
  * [15] J. E. Gough and M. R. James, “The series product and its application to quantum feedforward and feedback networks,” _IEEE Transactions on Automatic Control_ , vol. 54, no. 11, pp. 2530–2544, 2009.
  * [16] J. Combes, J. Kerckhoff, and M. Sarovar, “The SLH framework for modeling quantum input-output networks,” _Advances in Physics: X_ , vol. 2, no. 3, pp. 784–888, 2017.
  * [17] C. Gardiner and P. Zoller, _Quantum Noise: A Handbook of Markovian and Non-Markovian Quantum Stochastic Methods with Applications to Quantum Optics_. Springer Science & Business Media, 2004, vol. 56.
  * [18] H. M. Wiseman and G. J. Milburn, _Quantum Measurement and Control_. Cambridge university press, 2010.
  * [19] J. Zhang, Y.-X. Liu, R.-B. Wu, K. Jacobs, and F. Nori, “Quantum feedback: theory, experiments, and applications,” _Physics Reports_ , vol. 679, pp. 1–60, 2017.
  * [20] R. Loudon, _The Quantum Theory of Light_. OUP Oxford, 2000.
  * [21] G. J. Milburn, “Coherent control of single photon states,” _The European Physical Journal Special Topics_ , vol. 159, no. 1, pp. 113–117, Jun 2008.
  * [22] S. Fan, S. E. Kocabas, and J. T. Shen, “Input-output formalism for few-photon transport in one-dimensional nanophotonic waveguides coupled to a qubit,” _Physical Review A_ , vol. 82, p. 063821, 2010.
  * [23] J. E. Gough, M. R. James, H. I. Nurdin, and J. Combes, “Quantum filtering for systems driven by fields in single-photon states or superposition of coherent states,” _Physical Review A_ , vol. 86, no. 4, p. 043819, 2012.
  * [24] G. Zhang and M. R. James, “On the response of quantum linear systems to single photon input fields,” _IEEE Transactions on Automatic Control_ , vol. 58, no. 5, pp. 1221–1235, 2013.
  * [25] J. E. Gough and G. Zhang, “Generating nonclassical quantum input field states with modulating filters,” _EPJ Quantum Technology_ , vol. 2, pp. 2–15, 2015\.
  * [26] Y. Pan, G. Zhang, and M. R. James, “Analysis and control of quantum finite-level systems driven by single-photon input states,” _Automatica_ , vol. 69, pp. 18–23, 2016.
  * [27] Y. Pan, D .Dong, and G. Zhang, “Exact analysis of the response of quantum systems to two photons using a QSDE approach,” _New Journal of Physics_ , vol. 18, pp. 033004, 2016.
  * [28] A. O. Davis, V. Thiel, M. Karpiński, and B. J. Smith, “Measuring the single-photon temporal-spectral wave function,” _Physical Review Letters_ , vol. 121, no. 8, p. 083602, 2018.
  * [29] J. E. Gough and G. Zhang, “On realization theory of quantum linear systems,” _Automatica_ , vol. 59, pp. 139–151, 2015.
  * [30] H. I. Nurdin and N. Yamamoto, Linear Dynamical Quantum Systems: Analysis, Synthesis, and Control, Springer, 2017.
  * [31] G. Zhang, S. Grivopoulos, I. R. Petersen, and J. E. Gough, “The Kalman decomposition for linear quantum systems,” _IEEE Transactions on Automatic Control_ , vol. 63, no. 2, pp. 331–346, 2018.
  * [32] G. Zhang, I. R. Petersen, and J. Li, “Structural characterization of linear quantum systems with application to back-action evading measurement,” _IEEE Transactions on Automatic Control_ , vol. 65, pp. 3157–3163, 2020.
  * [33] G. Zhang and M. R. James, “Direct and indirect couplings in coherent feedback control of linear quantum systems,” _IEEE Transactions on Automatic Control_ , vol. 56, pp. 1535–1550, 2011.
  * [34] N. Yamamoto, “Decoherence-free linear quantum subsystems,” _IEEE Transactions on Automatic Control_ , vol. 59, no. 7, pp. 1845–1857, 2014.
  * [35] F. Ticozzi and L. Viola, “Quantum markovian subsystems: invariance, attractivity, and control,” _IEEE Transactions on Automatic Control_ , vol. 53, no. 9, pp. 2048–2063, 2008.
  * [36] Y. Wang, J. Minář, L. Sheridan, and V. Scarani, “Efficient excitation of a two-level atom by a single photon in a propagating mode,” _Physical Review A_ , vol. 83, no. 6, p. 063842, 2011.
  * [37] A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, “Circuit quantum electrodynamics,” _Review of Modern Physics_ , vol. 93, p. 025005, 2021.
  * [38] M. Pierre, I.-M. Svensson, S. Raman Sathyamoorthy, G. Johansson, and P. Delsing, “Storage and on-demand release of microwaves using superconducting resonators with tunable coupling,” _Applied Physics Letters_ , vol. 104, no. 23, p. 232604, 2014.
  * [39] K. A. Fischer, L. Hanschke, J. Wierzbowski, T. Simmet, C. Dory, J. J. Finley, J. Vučković, and K. Müller, “Signatures of two-photon pulses from a quantum two-level system,” _Nature Physics_ , vol. 13, no. 7, pp. 649–654, 2017.
  * [40] K. A. Fischer, R. Trivedi, V. Ramasesh, I. Siddiqi, and J. Vučković, “Scattering into one-dimensional waveguides from a coherently-driven quantum-optical system,” _Quantum_ , vol. 2, p. 69, 2018.
  * [41] R. Trivedi, K. Fischer, S. Xu, S. Fan, and J. Vuckovic, “Few-photon scattering and emission from low-dimensional quantum systems,” _Physical Review B_ , vol. 98, no. 14, p. 144112, 2018.
  * [42] W. J. Rugh, _Linear System Theory_. Prentice-Hall, Inc., 1996.
  * [43] O. O. Chumak and E. V. Stolyarov, “Photon distribution function for propagation of two-photon pulses in waveguide-qubit systems,” _Physical Review A_ , vol. 90, p. 063832, 2014.
  * [44] K. Zhou, J. Doyle, and K. Glover,  _Robust and Optimal Control,_ Upper Saddle River, NJ:Prentice-Hall, 1996.