# Generation of time-frequency entangled photon pairs propagating in separate
waveguides in circuit QED setup

E. V. Stolyarov  0000-0001-7531-5953<EMAIL_ADDRESS>Bogolyubov
Institute for Theoretical Physics, National Academy of Sciences of Ukraine,
Metrolohichna Street 14-b, 03143 Kyiv, Ukraine

###### Abstract

Time-frequency entangled photons constitute an important resource for a
plethora of applications across the diverse quantum technology landscape.
Thus, efficient and tunable setups for the generation of entangled photons are
requisite for modern quantum technologies. In this work, we propose a generic
cavity QED setup designed for on-demand generation of time-frequency entangled
photon pairs, with each photon propagating in a separate waveguide. We outline
a potential incarnation of this setup in the microwave superconducting circuit
QED architecture. We derive and numerically solve the set of equations of
motion governing the evolution of the quantum state of the system, allowing us
to examine the photon emission dynamics. Using the Schmidt decomposition of
the joint spectral amplitude of the emitted photon pair, we compute the
entanglement entropy analyzing its dependence on the system parameters. We
outline the potential extension of the proposed scheme for the generation of
multiphoton time-frequency entangled states.

##  I Introduction

The photonic time-frequency (TF) degree of freedom can serve as a basis for
encoding quantum information, including time [1], frequency [2], and temporal-
mode [3, *ansari2018] encodings. Moreover, the TF degree is intrinsically
robust for long-distance transmission of quantum information via microwave and
optical waveguides. Nonclassical spectral and temporal correlations of
photonic TF entangled states can be harnessed for boosting the sensitivity and
resolution of the two-photon absorption spectroscopy [5, 6, *jav1990,
*tabakaev2021, 9, *oka2020] and ultrafast TF-resolved Raman spectroscopy [11]
and enhancing the detection efficiency of low-reflectivity objects in the
bright thermal environment [12, 13, *barz2020]. The large dimensionality of TF
entangled states is instrumental in the realization of high-dimensional
quantum key distribution protocols [15, 16, *nunn2013]. Besides, TF entangled
states proved to be useful in a range of other applications, such as testing
Bell inequalities [18, 19, 20], quantum lithography [21], quantum-enhanced
clock synchronization [22], and biomedical imaging [23]. However, most of the
experimental implementations of these quantum protocols and techniques use
infrared and visible-light photons. Having an efficient and scalable scheme
for the generation of TF entangled photons in the microwave domain, one can
leverage that diverse quantum-optical toolbox in the elaborate microwave
photonic platform based on superconducting quantum circuits [24]. Rapid
progress and ongoing advances in the development of microwave superconducting
circuit QED systems demonstrated the high potential of these quantum systems
as a versatile hardware architecture for the realization of quantum
information processing (QIP) devices [25, *blais2020].

A number of theoretical proposals [27, 28, 29, 30] and experimental
demonstrations [31, 32, 33, 34, 35, 36, 37] of various designs of
deterministic and probabilistic circuit QED sources of TF entangled microwave
photons were put forward. These systems exploit the nonlinear properties of
the Josephson junctions. One of the approaches consists in using voltage-
biased Josephson junctions interacting with modes of open microwave
resonators. Due to strong nonlinear light-charge interaction, the inelastic
tunneling of Cooper pairs across the Josephson junction can generate
nonclassical states of light [31, 32, 33, 38, 39, *grimm2019, 41, *ma2022],
particularly TF entangled photons [31, 32, 33]. Alternatively, one can employ
the nondegenerate Josephson mixer – a circuit QED device, which parametrically
couples two modes with different frequencies by a pump at their sum-frequency
[43, 44]. In Ref. [34], the authors demonstrated the generation of entangled
microwave radiation over two separate transmission lines using the
nondegenerate Josephson mixer. A Josephson metamaterial operating in a Kerr-
free three-wave mixing mode was utilized for the generation of TF entangled
microwave photons in Ref. [37]. Another approach is leveraging superconducting
artificial atoms, which are anharmonic multilevel quantum systems [45, 46]. By
populating higher excited states of an artificial atom via driving the two-
photon transition, one can trigger the process of spontaneous cascaded
relaxation leading to the emission of TF entangled photons [27, 28, 30, 29,
35, 36].

Figure 1: (a) Schematic representation of the system comprised of a pair of
2LEs coupled to a mode of an open resonator. (b) The equivalent representation
of the system depicted in Fig. 1(a) – a 2LE interacting with a V3LE coupled to
a waveguide. (c) Scheme of the system composed of a resonator interacting with
a $\displaystyle\textstyle\Xi$-configuration 3LE and the diagram showing the
interrelation between the $\displaystyle\textstyle\Xi$3LE transition
frequencies and the resonator frequency. (d) The generic schematics of the
proposed source of TF entangled photon pairs. (e) Diagram of the interrelation
between the frequencies of the resonators and the transition frequencies of
the emitters in the setup demonstrated in Fig. 1(d).

In the paper, we start with proposing a general scheme of the setup designed
for the generation of TF entangled photons and then outline its potential
experimental implementation as a microwave photonic device based on a
superconducting circuit QED platform. In this scheme, the correlated
relaxation of a pair of excited quantum emitters is realized when the
relaxation of one of the emitters into the open resonator triggers the
relaxation of its counterpart leading to the emission of a pair of TF
entangled photons. One of the advantages of the proposed scheme for generation
of TF entangled photon pairs is that it does not require driving the two-
photon transitions of the artificial atoms. Moreover, it can be extended for
the generation of multiphoton entangled states. We focus on studying the
proposed device in the pulsed on-demand regime, when exactly only one pair of
entangled photons are emitted. Using a state-vector approach, we provide a
real-time picture of photon emission dynamics. We use the Schmidt
decomposition of the joint spectral amplitude of the emitted photons for
computing the von Neumann (entanglement) entropy, which is employed as a
measure of their bipartite entanglement.

The structure of the paper is as follows. In Sec. II, we discuss the general
principles determining the operation of the proposed source of TF entangled
photon pairs. Its potential implementation within the microwave circuit QED
architecture is outlined in Sec. III. In Sec. IV, we present the model
Hamiltonian describing the system. In Sec. V, we derive the equations of
motion governing the quantum-state evolution of the system and study the
emission dynamics. The dependence of the entanglement degree of the emitted
photons on the system parameters is analyzed in Sec. VI. We summarize our
results and outline possible applications and extensions of the proposed
source of entangled photon pairs in Sec. VII. The additional considerations
and details of derivations are delegated to the Appendixes.

##  II General scheme and operational principle

To explain the key principles defining the operation of the proposed source of
TF entangled photon pairs, let us consider a paradigmatic system schematically
presented in Fig. 1(a). This system consists of an open resonator interacting
with a pair of two-level emitters (2LEs). The first 2LE is coupled to the
resonator with strength $\displaystyle\textstyle g_{1}$, and the second 2LE
interacts with the resonator with strength $\displaystyle\textstyle g_{2}$.
There is no direct interaction between the emitters. For simplicity, we assume
that the resonator and the 2LEs have identical frequencies
$\displaystyle\textstyle\omega$. The resonator is open due to its coupling to
a waveguide, which leads to the photon leakage from the resonator with rate
$\displaystyle\textstyle\kappa$. Initially, the first 2LE resides in the
ground state, while the second 2LE is prepared in the excited state. Due to
interaction with the open resonator, the second 2LE eventually decays,
delivering the excitation (photon) to the waveguide.

The system, shown in Fig. 1(a), can be represented as the excited-state 2LE
interacting with the V-configuration three-level emitter (3LE) coupled to the
waveguide. Figure 1(b) illustrates the equivalent representation, whose
detailed justification is provided in Appendix A. The V-type 3LE is formed by
the ground state $\displaystyle\textstyle|\mathcal{G}\rangle$ and a pair of
excited states $\displaystyle\textstyle|\mathcal{E}^{\pm}\rangle$. The latter
decay to the waveguide with the rate $\displaystyle\textstyle\kappa/2$. The
frequency of
$\displaystyle\textstyle|\mathcal{G}\rangle\leftrightarrow|\mathcal{E}^{\pm}\rangle$
transition is $\displaystyle\textstyle\omega\pm g_{1}$. Both transitions are
coupled to the 2LE with strength $\displaystyle\textstyle g_{2}/\sqrt{2}$. The
increase of coupling $\displaystyle\textstyle g_{1}$ leads to stronger
detuning of
$\displaystyle\textstyle|\mathcal{G}\rangle\leftrightarrow|\mathcal{E}^{\pm}\rangle$
transition from the 2LE frequency $\displaystyle\textstyle\omega$. For
$\displaystyle\textstyle g_{1}>g_{2}/\sqrt{2}$, the excitation exchange
between the 2LE and the V3LE becomes inefficient, resulting in a suppression
of the photon leakage into the waveguide. These qualitative considerations are
confirmed by the calculations revealing the suppression of the photon leakage
into the waveguide with the increase of the ratio $\displaystyle\textstyle
g_{1}/g_{2}$. Details of calculations are given in Appendix A.

The above result implies that the process of photon emission can be
manipulated by controlling the coupling parameter $\displaystyle\textstyle
g_{1}$. In this regard, let us consider a system comprised by a 3LE
interacting with a single-mode resonator. This system is sketched in Fig.
1(c). The ground state $\displaystyle\textstyle|g\rangle$ and the excited
states $\displaystyle\textstyle|e\rangle$ and
$\displaystyle\textstyle|f\rangle$ of the 3LE are arranged in a ladder
($\displaystyle\textstyle\Xi$-type) configuration implying that only
$\displaystyle\textstyle|g\rangle\leftrightarrow|e\rangle$ and
$\displaystyle\textstyle|e\rangle\leftrightarrow|f\rangle$ transitions are
allowed. The frequencies of the resonator and
$\displaystyle\textstyle|e\rangle\leftrightarrow|f\rangle$ transition
coincide, while $\displaystyle\textstyle|g\rangle\leftrightarrow|e\rangle$
transition is strongly detuned from the resonator. In this case, the
interaction between the 3LE and the resonator is approximately described by
the Hamiltonian
$\displaystyle\textstyle{\mathcal{H}}_{\mathrm{int}}=g\,(c^{\dagger}|e\rangle\langle
f|+|f\rangle\langle e|c)$, where $\displaystyle\textstyle g$ denotes the
coupling parameter, and the operator $\displaystyle\textstyle c$
($\displaystyle\textstyle c^{\dagger}$) annihilates (creates) a resonator
photon. One can formally rewrite the resonator-emitter interaction Hamiltonian
as $\displaystyle\textstyle{\mathcal{H}}_{\mathrm{int}}=g(|e\rangle\langle
e|+|f\rangle\langle f|)(c^{\dagger}|e\rangle\langle f|+|f\rangle\langle e|c)$.
Using the property $\displaystyle\textstyle|g\rangle\langle
g|+|e\rangle\langle e|+|f\rangle\langle f|=\mathbb{1}$, with
$\displaystyle\textstyle\mathbb{1}$ being the unity operator, leads to
$\displaystyle\textstyle{\mathcal{H}}_{\mathrm{int}}=\tilde{g}(c^{\dagger}|e\rangle\langle
f|+|f\rangle\langle e|c)$, where $\displaystyle\textstyle\tilde{g}\equiv
g(\mathbb{1}-|g\rangle\langle g|)$ can be interpreted as a resonator-emitter
coupling dependent on the $\displaystyle\textstyle\Xi$3LE ground state
population. Thus, if the 3LE resides in the ground state
$\displaystyle\textstyle|g\rangle$, its
$\displaystyle\textstyle|f\rangle\leftrightarrow|e\rangle$ transition is
decoupled from the resonator mode. This effect was utilized, e.g., in Ref.
[47] for the implementation of a continuous-wave single-photon transistor and
Ref. [48] for the implementation of a qubit-state-controlled single-photon
switch.

Harnessing the effects outlined above, we propose a setup consisting of two
emitters interacting with two resonators (marked as $\displaystyle\textstyle
A$ and $\displaystyle\textstyle B$) coupled to the output waveguides (marked
as $\displaystyle\textstyle\alpha$ and $\displaystyle\textstyle\beta$). The
considered setup is schematically shown in Fig. 1(e). In this setup, the first
emitter is coupled to both resonators. This emitter is represented by a
$\displaystyle\textstyle\Xi$-configuration 3LE with the ground state
$\displaystyle\textstyle|g\rangle_{1}$ and a pair of excited states
$\displaystyle\textstyle|e\rangle_{1}$ and
$\displaystyle\textstyle|f\rangle_{1}$. The frequency of
$\displaystyle\textstyle|g\rangle_{1}\leftrightarrow|e\rangle_{1}$ transition
coincides with the frequency $\displaystyle\textstyle\omega_{a}$ of resonator
$\displaystyle\textstyle A$, while the frequency of
$\displaystyle\textstyle|e\rangle_{1}\leftrightarrow|f\rangle_{1}$ transition
coincides with the frequency $\displaystyle\textstyle\omega_{b}$ of resonator
$\displaystyle\textstyle B$. The anharmonicity of the eigenlevels of the
$\displaystyle\textstyle\Xi$3LE results in the inhibition of the excitation
exchange between resonator $\displaystyle A$ and
$\displaystyle\textstyle|e\rangle_{1}\leftrightarrow|f\rangle_{1}$ transition
and between resonator $\displaystyle\textstyle B$ and
$\displaystyle\textstyle|g\rangle_{1}\leftrightarrow|e\rangle_{1}$ transition.
The second emitter is represented by a 2LE with the ground state
$\displaystyle\textstyle|g\rangle_{2}$ and the excited state
$\displaystyle\textstyle|e\rangle_{2}$. This emitter is coupled only to
resonator $\displaystyle\textstyle B$, and the frequency of its
$\displaystyle\textstyle|g\rangle_{2}\leftrightarrow|e\rangle_{2}$ transition
coincides with the resonator frequency $\displaystyle\textstyle\omega_{b}$.
The interrelation between the frequencies of the resonators and the
transitions frequencies of the emitters is shown in Fig. 1(d).

In the system outlined above, the _correlated relaxation_ of the emitters
emerges. The mechanism of this relaxation is as follows. When the 3LE resides
in its excited state $\displaystyle\textstyle|e\rangle_{1}$, the relaxation of
the 2LE into waveguide $\displaystyle\textstyle\beta$ via resonator
$\displaystyle\textstyle B$ is inhibited due to the mechanism outlined in the
second paragraph of this section. In the course of the 3LE relaxation from the
excited state $\displaystyle\textstyle|e\rangle_{1}$ to the ground state
$\displaystyle\textstyle|g\rangle_{1}$,
$\displaystyle\textstyle|f\rangle_{1}\leftrightarrow|e\rangle_{1}$ transition
of the 3LE is getting decoupled from the resonator $\displaystyle\textstyle
B$, as discussed earlier in this section, which triggers the relaxation of the
2LE excited state $\displaystyle\textstyle|e\rangle_{2}$. Hence, the photons
tend to be emitted into the waveguides one right after another, implying that
they exhibit _time correlation_ which, in turn, indicates their _frequency
(energy) anti-correlation_.

Figure 2: Sketch of the potential experimental circuit QED implementation of
the generic scheme illustrated in Fig. 1(e). The indices of the elements are
in correspondence to the scheme in Fig. 1(e). Each output transmission line is
additionally equipped with the Purcell filter (PF). The excited states of the
artificial atoms (Xmons) are individually prepared using the classical pulses
sent over XY control lines.

##  III Circuit QED setup

Having discussed the general principle of operation of the proposed source of
TF entangled photon pairs, let us now outline its potential implementation in
the circuit QED architecture. The superconducting circuit QED platform
provides all necessary components for implementation of the generic scheme
discussed in Sec. II, such as microwave resonators and transmission lines, as
well as Josephson-junction artificial atoms featuring tunable multilevel
structure and long coherence times [46].

The schematic illustration of the proposed circuit QED setup is shown in Fig.
2. In the considered setup, the quantum emitters are represented by the Xmon
version [49] of the transmon – a charge-type artificial atom derived from a
Cooper-pair box and featuring resilience to charge fluctuations [50,
*houck2009]. Transmon artificial atoms, in particular, their Xmon
modifications offer a simple design, long coherence times, and tunable
parameters [52], making them a popular choice for the realization of various
QIP systems [53, 54, 55]. Moreover, using tantalum as a base superconductor
along with improving processing techniques led to a recent impressive increase
in the transmon excited state lifetime [56, 57]. In our setup, the Xmons are
capacitively coupled to the coplanar waveguide (CPW) resonators which, in
turn, are coupled to the transmission-line waveguides serving as output
channels for the emitted photons.

The states of the transmons are manipulated using the classical signals
delivered via individual (XY) control lines [58]. Due to the weak
anharmonicity of transmons, for their preparation in the lower excited state,
one can rely on one of the elaborate techniques, providing high fidelities of
preparation while avoiding excitation of the higher excited states [59, 60,
61]. Alternatively, for transmon control, one can use a single-flux quantum
digital logic [62]. This approach allows integration of the control circuitry
on a chip bringing most of the setup components into the cryogenic stage.

The output transmission lines are supplemented with the Purcell filters for
suppression of the superfluous emission from the off-resonant transitions of
the transmons (more details on the role of these elements are given in
Appendix B.2). The setup presented in Fig. 2 adopts the approach from Ref.
[63] with the quarter-wavelength CPW resonators side-coupled to the output
transmission lines. These additional resonators serve as notch filters.

In the setup version shown in Fig. 2, the fixed couplings between the circuit
elements are used. For achieving _in-situ_ tunability of the setup parameters
for better control over the characteristics of the emitted photons, instead of
fixed couplings set on a chip fabrication stage, one can implement tunable
couplings using, for example, superconducting quantum interference device
couplers [64, 65, *wu2018] controlled by an external flux.

##  IV The Model

The circuit QED system described above is modeled by the Hamiltonian

$\begin{split}{\mathcal{H}}_{\mathrm{sys}}=&{\mathcal{H}}_{\mathrm{t}}+{\mathcal{H}}_{\mathrm{r}}+{\mathcal{H}}_{\mathrm{r}-\mathrm{t}}+{\mathcal{H}}_{\mathrm{w}}+{\mathcal{H}}_{\mathrm{r}-\mathrm{w}},\end{split}$
(1)

where the first term describes the transmon-type artificial atoms, the second
term is the Hamiltonian of the resonators, the third term describes the
interaction between the resonators and transmons, and the last pair of terms
describes the waveguides and their coupling to the resonators. In what follows
we set $\displaystyle\textstyle\hbar=1$, thus, measuring all energies in
frequency units.

The Hamiltonian $\displaystyle\textstyle{\mathcal{H}}_{\mathrm{t}}$
represented in the basis of the transmons eigenstates reads as

${\mathcal{H}}_{\mathrm{t}}=\sum_{j=1}^{2}\left(\omega_{j}^{ge}\sigma^{ee}_{j}+\omega^{gf}_{j}\sigma^{ff}_{j}\right),$
(2)

where $\displaystyle\textstyle\omega^{ge}_{j}$,
$\displaystyle\textstyle\omega^{gf}_{j}=\omega_{j}^{ge}+\omega_{j}^{ef}$, and
$\displaystyle\omega^{ef}_{j}$ stand for the frequency of
$\displaystyle\textstyle|g\rangle_{j}\leftrightarrow|e\rangle_{j}$,
$\displaystyle\textstyle|g\rangle_{j}\leftrightarrow|f\rangle_{j}$, and
$\displaystyle\textstyle|e\rangle_{j}\leftrightarrow|f\rangle_{j}$ transition
of the $\displaystyle\textstyle j$-th transmon ($\displaystyle\textstyle
j\in\\{1,2\\}$), respectively. For convenience, in Eq. (2) we introduced the
transmon operator
$\displaystyle\textstyle\sigma_{j}^{ss^{\prime}}=|s\rangle_{j}\langle
s^{\prime}|_{j}$ where $\displaystyle\textstyle s,s^{\prime}\in\\{g,e,f\\}$ is
the index of the transmon eigenlevel. Since we consider the case of only two
excitations in the system, we restricted the subspace of the transmon states
to the ground state $\displaystyle\textstyle|g\rangle_{j}$ and a pair of the
lowest excited states $\displaystyle\textstyle|e\rangle_{j}$ and
$\displaystyle\textstyle|f\rangle_{j}$.

The Hamiltonian of the resonators
$\displaystyle\textstyle{\mathcal{H}}_{\mathrm{r}}$ has the form

${\mathcal{H}}_{\mathrm{r}}=\omega_{a}{\mathcal{N}}_{a}+\omega_{b}{\mathcal{N}}_{b},$
(3)

where $\displaystyle\textstyle\omega_{a}$ and
$\displaystyle\textstyle\omega_{b}$ denote the frequencies of resonator
$\displaystyle\textstyle A$ and $\displaystyle\textstyle B$, respectively.
Operators $\displaystyle\textstyle{\mathcal{N}}_{a}=a^{\dagger}a$ and
$\displaystyle\textstyle{\mathcal{N}}_{b}=b^{\dagger}b$ are the operators of
the photon number in resonator $\displaystyle\textstyle A$ and
$\displaystyle\textstyle B$, where $\displaystyle\textstyle
a$($\displaystyle\textstyle a^{\dagger}$) and $\displaystyle\textstyle b$
($\displaystyle\textstyle b^{\dagger}$) stand for the photon annihilation
(creation) operators in the corresponding resonator.

For clarity, we represent the resonator-transmons interaction Hamiltonian as
$\displaystyle\textstyle{\mathcal{H}}_{\mathrm{r}-\mathrm{t}}={\mathcal{H}}^{\mathrm{rsn}}_{\mathrm{r}-\mathrm{t}}+{\mathcal{H}}^{\mathrm{dsp}}_{\mathrm{r}-\mathrm{t}}$,
where the first term describes the _resonant_ interactions determining the
dynamics of the system, while the second term describes the _dispersive_
interactions giving rise to shifts in the frequencies of the resonators and
the transmons. In the rotating-wave approximation, the term
$\displaystyle\textstyle{\mathcal{H}}^{\mathrm{res}}_{\mathrm{r-a}}$ has the
form

${\mathcal{H}}^{\mathrm{rsn}}_{\mathrm{r}-\mathrm{t}}=g_{1a}a^{\dagger}\sigma_{1}^{ge}+g_{1b}b^{\dagger}\sigma_{1}^{ef}+g_{2b}b^{\dagger}\sigma_{2}^{ge}+\textrm{H.c.},$
(4)

and the dispersive term
$\displaystyle\textstyle{\mathcal{H}}^{\mathrm{disp}}_{\mathrm{r-a}}$ reads as

${\mathcal{H}}^{\mathrm{dsp}}_{\mathrm{r}-\mathrm{t}}=\eta_{1a}a^{\dagger}\sigma_{1}^{ef}+\eta_{1b}b^{\dagger}\sigma_{1}^{ge}+\eta_{2b}b^{\dagger}\sigma_{2}^{ef}+\textrm{H.c.}.$
(5)

Parameters $\displaystyle\textstyle g_{1a}$ and $\displaystyle\textstyle
g_{1b}$ are the coupling strengths of
$\displaystyle\textstyle|g\rangle_{1}\leftrightarrow|e\rangle_{1}$ and
$\displaystyle\textstyle|e\rangle_{1}\leftrightarrow|f\rangle_{1}$ transitions
of the first transmon to resonators $\displaystyle\textstyle A$ and
$\displaystyle\textstyle B$, respectively. Parameter $\displaystyle\textstyle
g_{2b}$ stands for the coupling strength of
$\displaystyle\textstyle|g\rangle_{2}\leftrightarrow|e\rangle_{2}$ transition
of the second transmon to resonator $\displaystyle\textstyle B$. Parameter
$\displaystyle\textstyle\eta_{1a}$ is the coupling strength of resonator
$\displaystyle\textstyle A$ to the first transmon
$\displaystyle\textstyle|e\rangle_{1}\leftrightarrow|f\rangle_{1}$ transition,
$\displaystyle\textstyle\eta_{1b}$ stands for the coupling strength of
resonator $\displaystyle\textstyle B$ to
$\displaystyle\textstyle|g\rangle_{1}\leftrightarrow|e\rangle_{1}$ transition
of the first transmon, and $\displaystyle\textstyle\eta_{2b}$ denotes the
coupling strength of resonator $\displaystyle\textstyle B$ to
$\displaystyle\textstyle|e\rangle_{2}\leftrightarrow|f\rangle_{2}$ transition
of the second transmon. For the transmon-type artificial atoms, one has
$\displaystyle\textstyle\eta_{1a}\approx\sqrt{2}g_{1a}$,
$\displaystyle\textstyle\eta_{1b}\approx g_{1b}/\sqrt{2}$, and
$\displaystyle\textstyle\eta_{2b}\approx\sqrt{2}g_{2b}$ [50].

The waveguides are modeled by the continua of non-interacting bosonic modes
with the Hamiltonian $\displaystyle\textstyle{\mathcal{H}}_{\mathrm{w}}$ given
by

${\mathcal{H}}_{\mathrm{w}}=\int^{\infty}_{0}\mathrm{d}\nu\nu\left(\alpha^{\dagger}_{\nu}\alpha_{\nu}+\beta^{\dagger}_{\nu}\beta_{\nu}\right),$
(6)

where $\displaystyle\textstyle\alpha_{\nu}(a^{\dagger}_{\nu})$ and
$\displaystyle\textstyle\beta_{\nu}(\beta^{\dagger}_{\nu})$ stand for the
annihilation (creation) operators of the photon with frequency
$\displaystyle\nu$ propagating in the waveguide $\displaystyle\alpha$ and
$\displaystyle\beta$, respectively.

The Hamiltonian $\displaystyle{\mathcal{H}}_{\mathrm{r}-\mathrm{w}}$,
describing the coupling between the resonators and waveguides, has the form:

${\mathcal{H}}_{\mathrm{r}-\mathrm{w}}=\int^{\infty}_{0}\mathrm{d}\nu\left[f_{a}(\nu)a^{\dagger}\alpha_{\nu}+f_{b}(\nu)b^{\dagger}\beta_{\nu}+\textrm{H.c.}\right],$
(7)

where $\displaystyle\textstyle f_{a}(\nu)$ and $\displaystyle\textstyle
f_{b}(\nu)$ are the coupling strengths of resonator $\displaystyle\textstyle
A$ to waveguide $\displaystyle\textstyle\alpha$ and resonator
$\displaystyle\textstyle B$ to waveguide $\displaystyle\textstyle\beta$,
correspondingly.

Note that the model Hamiltonian in Eq. (1) does not explicitly include the
Purcell filters. Their contribution is absorbed into the frequency-dependent
resonator-waveguide couplings $\displaystyle f_{a}(\nu)$ and $\displaystyle
f_{b}(\nu)$. The effect of the Purcell filters is outlined in Appendix B.2.

Figure 3: The snapshots of the joint spectra of the emitted photon pairs for
the different moments of time. Parameters of the system are as follows:
$\displaystyle\textstyle(\omega^{ge}_{1}-\omega^{ef}_{1})/2\pi=400\,\mathrm{MHz}$,
$\displaystyle\textstyle\kappa_{a}=\kappa_{b}=\kappa$,
$\displaystyle\textstyle\kappa/2\pi=25\,\mathrm{MHz}$,
$\displaystyle\textstyle g_{1a}/2\pi=5\,\mathrm{MHz}$,
$\displaystyle\textstyle g_{2b}/2\pi=10\,\mathrm{MHz}$, (upper row)
$\displaystyle\textstyle g_{1b}=0$; (middle row) $\displaystyle\textstyle
g_{1b}/2\pi=10\,\mathrm{MHz}$; (bottom row) $\displaystyle\textstyle
g_{1b}/2\pi=25\,\mathrm{MHz}$. For aiding the visualization, each plot were
normalized on the corresponding maximum value of the joint spectrum
$\displaystyle\textstyle\max\\{S_{\nu,\nu^{\prime}}(t)\\}$ and then multiplied
on the probability $\displaystyle\textstyle p_{\alpha\beta}(t)$ of finding
both waveguides containing photons at time $\displaystyle\textstyle t$. The
corresponding values of $\displaystyle\textstyle p_{\alpha\beta}(t)$ are shown
in each plot.

### The effective Hamiltonian

The dispersive coupling term
$\displaystyle\textstyle{\mathcal{H}}^{\mathrm{dsp}}_{\mathrm{r-a}}$ given by
Eq. (5) can be treated as a perturbation of the Hamiltonian
$\displaystyle\textstyle{\mathcal{H}}_{\mathrm{sys}}$ provided that
$\displaystyle\textstyle|\lambda_{1a}|,|\lambda_{1b}|,|\lambda_{2b}|\ll 1$,
where

$\lambda_{1a}=\frac{\eta_{1a}}{\omega_{1}^{ef}-\omega_{a}},\,\lambda_{1b}=\frac{\eta_{1b}}{\omega_{1}^{ge}-\omega_{b}},\,\lambda_{2b}=\frac{\eta_{2b}}{\omega_{2}^{ef}-\omega_{b}}.$

In this case, one can eliminate
$\displaystyle\textstyle{\mathcal{H}}^{\mathrm{dsp}}_{\mathrm{r-a}}$ via the
Schrieffer-Wolff-type transformation
$\displaystyle\textstyle{\mathcal{H}}_{\mathrm{sys}}\rightarrow{\mathcal{U}}^{\dagger}{\mathcal{H}}_{\mathrm{sys}}{\mathcal{U}}$
with $\displaystyle{\mathcal{U}}$ being the unitary operator given by [67, 68]

${\mathcal{U}}=\exp\left(\lambda_{1a}a^{\dagger}\sigma_{1}^{ef}+\lambda_{1b}b^{\dagger}\sigma_{1}^{ge}+\lambda_{2b}b^{\dagger}\sigma_{2}^{ef}-\textrm{H.c.}\right).$
(8)

Keeping the terms contributing up to the first order in the small parameters
$\displaystyle\textstyle\lambda_{1a}$, $\displaystyle\textstyle\lambda_{1b}$,
and $\displaystyle\textstyle\lambda_{2b}$, one arrives at the effective
Hamiltonian $\displaystyle\textstyle\mathcal{H}_{\mathrm{eff}}$ of the form
(see details of derivation in Appendix B):

$\mathcal{H}_{\mathrm{eff}}={\bar{\mathcal{H}}}_{\mathrm{t}}+{\bar{\mathcal{H}}}_{\mathrm{r}}+{\mathcal{H}}^{\mathrm{rsn}}_{\mathrm{r}-\mathrm{t}}+{\mathcal{H}}_{\mathrm{rr-
tm}}+{\mathcal{H}}_{\mathrm{w}}+{\mathcal{H}}_{\mathrm{r}-\mathrm{w}}.$ (9)

In the dressed basis, the Hamiltonian describing the transmons reads as

$\begin{split}{\bar{\mathcal{H}}}_{\mathrm{t}}=\bar{\omega}_{1}^{ge}\sigma_{1}^{ee}+\bar{\omega}_{1}^{gf}\sigma_{1}^{ff}+\omega_{2}^{ge}\sigma_{1}^{ee}+\bar{\omega}_{2}^{gf}\sigma_{2}^{ff},\end{split}$
(10)

with $\displaystyle\textstyle\bar{\omega}_{1}^{ge}=\omega_{1}^{ge}+\chi_{1b}$,
$\displaystyle\textstyle\bar{\omega}_{1}^{gf}=\bar{\omega}_{1}^{ge}+\bar{\omega}_{1}^{ef}=\omega_{1}^{ef}+\chi_{1a}$,
and
$\displaystyle\textstyle\bar{\omega}_{2}^{gf}=\omega_{2}^{ge}+\bar{\omega}_{2}^{ef}=\omega_{2}^{gf}+\chi_{2b}$
being the dressed (renormalized) frequencies of
$\displaystyle\textstyle|g\rangle_{1}\leftrightarrow|e\rangle_{1}$ and
$\displaystyle\textstyle|g\rangle_{1}\leftrightarrow|f\rangle_{1}$ transitions
of the first transmon, and
$\displaystyle\textstyle|g\rangle_{2}\leftrightarrow|f\rangle_{2}$ transition
of the second transmon, respectively. Here, we introduced the notations
$\displaystyle\textstyle\chi_{1a}=\lambda_{1a}\eta_{1a}$,
$\displaystyle\textstyle\chi_{1b}=\lambda_{1b}\eta_{1b}$, and
$\displaystyle\textstyle\chi_{2b}=\lambda_{2b}\eta_{2b}$.

The Hamiltonian of the resonators acquires the form

$\begin{split}{\bar{\mathcal{H}}}_{\mathrm{r}}=&\,\big{(}\omega_{a}+\chi_{1a}\mathcal{Z}_{1}^{fe}\big{)}a^{\dagger}a\\\
&\,+\big{(}\omega_{b}+\chi_{1b}\mathcal{Z}_{1}^{eg}+\chi_{2b}\mathcal{Z}_{2}^{fe}\big{)}b^{\dagger}b,\end{split}$
(11)

where we introduced the notation
$\displaystyle\textstyle\mathcal{Z}_{j}^{ss^{\prime}}\equiv\sigma_{j}^{ss}-\sigma_{j}^{s^{\prime}s^{\prime}}$.

The term
$\displaystyle\textstyle{\mathcal{H}}_{\mathrm{r}\mathrm{r}-\mathrm{t}}$ in
Eq. (9) is given by

${\mathcal{H}}_{\mathrm{r}\mathrm{r}-\mathrm{t}}=\varUpsilon\big{(}a^{\dagger}b^{\dagger}\sigma_{1}^{gf}+\sigma_{1}^{fg}a\,b\big{)},$
(12)

where $\displaystyle\textstyle
2\varUpsilon=\lambda_{1a}\eta_{1b}-\lambda_{1b}\eta_{1a}$. This term describes
the processes of simultaneous exchange of two excitations between the first
transmon and resonators $\displaystyle\textstyle A$ and
$\displaystyle\textstyle B$.

In the further analysis, we work in the dressed basis and describe the
dynamics of the system using the Hamiltonian
$\displaystyle\textstyle\mathcal{H}_{\mathrm{eff}}$ expressed by Eqs.
(9)–(12).

Note that here we focus on the unitary dynamics of the system, neglecting the
dissipation. We assume that dissipation processes occur on the timescales much
longer than the timescales of the coherent processes in the system. With such
an approach, we simplify the theoretical treatment of the considered system
while still grasping the essential features of its quantum dynamics.

Indeed, for the system parameters we work with, the photon leakage rates from
the resonators to the output waveguides are
$\displaystyle\textstyle\kappa/2\pi>1\,\mathrm{MHz}$, while the internal
quality factors of the CPW resonators can exceed $\displaystyle\textstyle
10^{6}$ [69, 70], giving the intrinsic resonator photon loss rates
$\displaystyle\textstyle\gamma_{\mathrm{res}}/2\pi\lesssim
0.01\,\mathrm{MHz}$. The coaxial resonators [71, 72] and the three-dimensional
microwave cavities [73, 74, 75] offer even higher internal quality factors
reaching up to $\displaystyle 10^{8}-10^{9}$. Thus, the process of photon
leakage to the output waveguides dominates over the photon dissipation inside
the resonators. The typical excited-state lifetimes of the state-of-the-art
transmons reach $\displaystyle\textstyle 0.05-0.1\,\mathrm{ms}$ [49, 55], with
the recent experiments [56, 57] reporting lifetimes up to
$\displaystyle\textstyle 0.5\,\mathrm{ms}$, which are more than an order of
magnitude longer than the times the transmons dwell in their excited states
for the system parameters we consider in this work. Thus, the decoherence of
transmons has a negligible effect on the emission dynamics.

##  V Emission dynamics

We assume that initially (at $\displaystyle\textstyle t=0$), the system is
prepared in the state $\displaystyle\textstyle|\Psi_{\mathrm{in}}\rangle$ with
the transmons residing in their excited states
$\displaystyle\textstyle|e\rangle_{1}$ and
$\displaystyle\textstyle|e\rangle_{2}$, while the resonators and the
waveguides are void of excitation. Thus, the initial state of the system reads
as
$\displaystyle\textstyle|\Psi_{\mathrm{in}}\rangle=\sigma_{1}^{eg}\sigma_{2}^{eg}|\varnothing\rangle$,
where
$\displaystyle\textstyle|\varnothing\rangle=|\varnothing\rangle_{\alpha}|\varnothing\rangle_{\beta}|0\rangle_{a}|0\rangle_{b}|g\rangle_{1}|g\rangle_{2}$
stands for the vacuum state of the system – a state with a void of photons in
the resonators and the waveguides and the transmons residing in their ground
states. In our analysis, we do not account for the effect of thermal
excitations on the dynamics of the system. This simplification is justified
since the superconducting circuit QED systems usually operate at temperatures
$\displaystyle\textstyle T_{\mathrm{sys}}\sim 10-30\,\mathrm{mK}$, while the
typical working frequencies of microwave circuit QED setups are
$\displaystyle\omega_{\mathrm{sys}}/2\pi\sim 5-20\,\mathrm{GHz}$ [24, 58]. For
these parameters, one obtains the estimate $\displaystyle\textstyle
n_{\mathrm{th}}<10^{-3}$ for the average number of thermal excitations in the
system.

At the arbitrary moment of time $\displaystyle\textstyle t$, the state of the
system
$\displaystyle\textstyle|\varPsi(t)\rangle=\mathrm{e}^{-\mathrm{i}\mathcal{H}_{\mathrm{eff}}t}|\Psi_{\mathrm{in}}\rangle$
is expressed as

$\begin{split}|\varPsi(t)\rangle=&\,\int^{\infty}_{0}\mathrm{d}\nu\int^{\infty}_{0}\mathrm{d}\nu^{\prime}\varPhi_{\nu,\nu^{\prime}}(t)\alpha^{\dagger}_{\nu}\beta^{\dagger}_{\nu^{\prime}}|\varnothing\rangle+\int^{\infty}_{0}\mathrm{d}\nu\left[\varXi^{\alpha}_{\nu}(t)\alpha^{\dagger}_{\nu}b^{\dagger}+\varXi^{\beta}_{\nu}(t)\beta^{\dagger}_{\nu}a^{\dagger}+\varTheta^{\alpha}_{\nu}(t)\alpha^{\dagger}_{\nu}\sigma_{2}^{eg}+\varTheta^{\beta}_{\nu}(t)\beta^{\dagger}_{\nu}\sigma_{1}^{eg}\right]|\varnothing\rangle\\\
&\,+R(t)a^{\dagger}b^{\dagger}|\varnothing\rangle+Q(t)\sigma_{1}^{fg}|\varnothing\rangle+Y_{a}(t)a^{\dagger}\sigma_{2}^{eg}|\varnothing\rangle+Y_{b}(t)b^{\dagger}\sigma_{1}^{eg}|\varnothing\rangle+X(t)\sigma_{1}^{eg}\sigma_{2}^{eg}|\varnothing\rangle.\end{split}$
(13)

The first term in Eq. (13) corresponds to the state of the system with both
waveguides hosting a photon, with
$\displaystyle\textstyle\varPhi_{\nu,\nu^{\prime}}$ being the joint spectral
amplitude. The quantity $\displaystyle\textstyle
S_{\nu,\nu}(t)=|\varPhi_{\nu,\nu^{\prime}}(t)|^{2}$ determines the joint
spectrum – the probability density distribution of finding a photon with
frequency $\displaystyle\textstyle\nu$ propagating in waveguide
$\displaystyle\textstyle\alpha$ and a photon with frequency
$\displaystyle\textstyle\nu^{\prime}$ propagating in waveguide
$\displaystyle\textstyle\beta$ at the moment of time $\displaystyle\textstyle
t$. The probability of finding the waveguides both accommodating photons is
given by

$p_{\alpha\beta}(t)=\int^{\infty}_{0}\mathrm{d}\nu\int^{\infty}_{0}\mathrm{d}\nu^{\prime}\,S_{\nu,\nu^{\prime}}(t).$
(14)

The remaining terms in the upper line of Eq. (13) correspond to the states
with photon propagating in only one of the waveguides:
$\displaystyle\textstyle\varXi_{\nu}^{\alpha}(\varXi_{\nu}^{\beta})$ is the
amplitude of the state with the photon in resonator $\displaystyle\textstyle
A(B)$ and the transmons in their ground states,
$\displaystyle\textstyle\varTheta_{\nu}^{\alpha}(\varTheta_{\nu}^{\beta})$ is
the amplitude of the state with the resonators void of photons and the first
(second) transmon residing in the excited state
$\displaystyle\textstyle|e\rangle_{1(2)}$. The terms in the bottom line of Eq.
(13) correspond to the states of the system with both waveguides void of
excitations.

The equations of motion governing the evolution of the probability amplitudes
in Eq. (13) can be compactly written as

$\frac{\partial\mathbf{\boldsymbol{\mu}}(t)}{\partial
t}=-\mathrm{i}\,\boldsymbol{\Omega}\,\boldsymbol{\mu}(t),$ (15)

where
$\displaystyle\textstyle\boldsymbol{\mu}(t)=[\varPhi_{\nu,\nu^{\prime}}(t),\varXi^{\alpha}_{\nu}(t),\varXi^{\beta}_{\nu}(t),\varTheta^{\alpha}_{\nu}(t),\varTheta^{\beta}_{\nu}(t),R(t),Q(t),Y_{a}(t),Y_{b}(t),X(t)]^{\mathrm{T}}$
and the matrix $\displaystyle\textstyle\boldsymbol{\Omega}$ reads as

$\boldsymbol{\Omega}=\begin{pmatrix}\nu+\nu^{\prime}&f_{b}(\omega_{b})&f_{a}(\omega_{a})&0&0&0&0&0&0&0\\\
0&\nu+\widetilde{\omega}_{b}-\chi_{1b}&0&g_{2b}&0&f_{a}(\omega_{a})&0&0&0&0\\\
0&0&\nu^{\prime}+\widetilde{\omega}_{a}&0&g_{1a}&f_{b}(\omega_{b})&0&0&0&0\\\
0&g_{2b}&0&\nu+\omega_{2}^{ge}&0&0&0&f_{a}(\omega_{a})&0&0\\\
0&0&g_{1a}&0&\nu^{\prime}+\bar{\omega}_{1}^{ge}&0&0&0&f_{b}(\omega_{b})&0\\\
0&0&0&0&0&\widetilde{\omega}_{a}+\widetilde{\omega}_{b}-\chi_{1b}&\varUpsilon&g_{2b}&g_{1a}&0\\\
0&0&0&0&0&\varUpsilon&\bar{\omega}^{gf}_{1}&0&g_{1b}&0\\\
0&0&0&0&0&g_{2b}&0&\widetilde{\omega}_{a}+\omega_{2}^{ge}&0&g_{1a}\\\
0&0&0&0&0&g_{1a}&g_{1b}&0&\widetilde{\omega}_{b}+\bar{\omega}_{1}^{ge}+\chi_{1b}&g_{2b}\\\
0&0&0&0&0&0&0&g_{1a}&g_{2b}&\bar{\omega}_{1}^{ge}+\omega_{2}^{ge}\\\
\end{pmatrix},$

where
$\displaystyle\textstyle\widetilde{\omega}_{a}=\omega_{a}-\mathrm{i}\kappa_{a}/2$
and
$\displaystyle\textstyle\widetilde{\omega}_{b}=\omega_{b}-\mathrm{i}\kappa_{b}/2$
with $\displaystyle\textstyle\kappa_{a}=2\pi f^{2}_{a}(\omega_{a})$ and
$\displaystyle\textstyle\kappa_{b}=2\pi f^{2}_{b}(\omega_{b})$ being the
photon leakage rate from resonators $\displaystyle\textstyle A$ and
$\displaystyle\textstyle B$ into waveguides $\displaystyle\textstyle\alpha$
and $\displaystyle\textstyle\beta$, correspondingly. The derivation of Eq.
(15) is given in Appendix C and follows the approach demonstrated in Refs.
[76] and [77].

For computing the probability amplitudes standing in Eq. (13), we use the
symbolic solution of Eq. (15) expressed as
$\displaystyle\textstyle\boldsymbol{\mu}(t)=\exp(-\mathrm{i}\boldsymbol{\Omega}t)\boldsymbol{\mu}(0)$.
The matrix exponential is evaluated numerically using the built-in function
MatrixExp of Mathematica.

Figure 3 provides the series of snapshots of the joint spectra of the emitted
photons evaluated for the specific moments of time and various parameters of
the system. The computations reveal that switching on the coupling between
resonator $\displaystyle\textstyle B$ and
$\displaystyle\textstyle|e\rangle_{1}\leftrightarrow|f\rangle_{1}$ transition
of the first transmon results in the pronounced frequency anti-correlation of
the emitted photons, which manifests itself as an accumulation of the joint
spectrum near the line given by
$\displaystyle\textstyle\nu+\nu^{\prime}=\omega_{a}+\omega_{b}$.

Figure 4: (Upper row) Joint spectra of the emitted photons
$\displaystyle\textstyle S_{\nu,\nu^{\prime}}(t)$ versus (bottom row) the
corresponding time-domain probability densities $\displaystyle\textstyle
A_{\tau,\tau^{\prime}}(t)$ at $\displaystyle\textstyle t=0.5\,\mathrm{\mu s}$
for different values of $\displaystyle\textstyle g_{1b}$: (left column)
$\displaystyle\textstyle g_{1b}=0$; (central column) $\displaystyle\textstyle
g_{1b}/2\pi=10\,\mathrm{MHz}$; (right column) $\displaystyle\textstyle
g_{1b}/2\pi=25\,\mathrm{MHz}$. The rest of the parameters are as described in
the caption of Fig. 3. The dashed lines correspond
$\displaystyle\textstyle\nu+\nu^{\prime}=\omega_{a}+\omega_{b}$ (upper row)
and $\displaystyle\textstyle\tau^{\prime}=\tau$ (bottom row). Dotted lines are
defined as $\displaystyle\textstyle\tau^{\prime}=\tau+\delta\tau$. Plots were
normalized on $\displaystyle\textstyle\max\\{S_{\nu,\nu^{\prime}}(t)\\}$
(upper row) and $\displaystyle\textstyle\max\\{A_{\tau,\tau^{\prime}}(t)\\}$
(bottom row) and then multiplied on the corresponding value of
$\displaystyle\textstyle p_{\alpha\beta}(t)$.

For the analysis of the temporal properties of the emitted photons, we
introduce their time-domain joint probability amplitude
$\displaystyle\textstyle\phi_{\tau,\tau^{\prime}}(t)$ expressed via the
Fourier transform of the joint spectral amplitude
$\displaystyle\textstyle\varPhi_{\nu,\nu^{\prime}}(t)$ as:

$\phi_{\tau,\tau^{\prime}}(t)=\frac{1}{2\pi}\int\mathrm{d}\nu\,\mathrm{e}^{\mathrm{i}\nu\tau}\int\mathrm{d}\nu^{\prime}\,\mathrm{e}^{\mathrm{i}\nu^{\prime}\tau^{\prime}}\,\varPhi_{\nu,\nu^{\prime}}(t).$
(16)

One can interpret the quantity $\displaystyle\textstyle
A_{\tau,\tau^{\prime}}(t)=|\phi_{\tau,\tau^{\prime}}(t)|^{2}$ as the
distribution of the joint probability density of finding one photon at point
$\displaystyle\textstyle x=\upsilon_{\alpha}t$ in waveguide
$\displaystyle\textstyle\alpha$ at instant $\displaystyle\textstyle\tau$ along
with finding another photon at point $\displaystyle\textstyle
x=\upsilon_{\beta}t$ in waveguide $\displaystyle\textstyle\beta$ at instant
$\displaystyle\textstyle\tau^{\prime}$, assuming that the coupling points of
the waveguides to the resonators are at $\displaystyle\textstyle x=0$. Here,
$\displaystyle\textstyle\upsilon_{\alpha}$ and
$\displaystyle\textstyle\upsilon_{\beta}$ stand for the photon group velocity
in waveguide $\displaystyle\textstyle\alpha$ and
$\displaystyle\textstyle\beta$, correspondingly. Figure 4 demonstrates the
joint spectra $\displaystyle\textstyle S_{\nu,\nu^{\prime}}(t)$ of the emitted
photons at $\displaystyle\textstyle t=0.5\,\mathrm{\mu s}$ versus the
respective time-domain probability densities $\displaystyle\textstyle
A_{\tau,\tau^{\prime}}(t)$ for different values of $\displaystyle\textstyle
g_{1b}$. One can notice that the frequency anti-correlation of the emitted
photons arising for $\displaystyle\textstyle g_{1b}\neq 0$ is accompanied by
their time correlation, which appears as a dilution of the time-domain
probability density for $\displaystyle\textstyle\tau^{\prime}\leq\tau$ and its
aggregation in the region $\displaystyle\textstyle\tau^{\prime}>\tau$ along
the line $\displaystyle\textstyle\tau^{\prime}=\tau+\delta\tau$. This result
indicates that when we switch on the interaction between the second transmon
and resonator $\displaystyle\textstyle B$, the emission of the photon into
waveguide $\displaystyle\textstyle\beta$ starts _after_ the emission of the
photon into waveguide $\displaystyle\textstyle\alpha$. The delay between the
photons is roughly estimated as
$\displaystyle\textstyle\delta\tau\approx(\kappa_{b}/2)^{-1}$. Such behavior
is consistent with the considerations presented in Sec. II.

Thus, the general idea discussed in Sec. II is supported by the numerical
results confirming that the proposed setup can emit the TF entangled photon
pairs. As a quantitative measure of photon entanglement, we use the von
Neumann entropy, whose definition, along with the extensive analysis of its
dependence of the parameters of the system, is given in Sec. VI.

Figure 5: Dynamics of the populations of the transmons eigenlevels (upper
row), resonators populations (middle row), and waveguides populations (bottom
row). Parameters of the system used for computations are the following: (left
column) $\displaystyle\textstyle g_{1b}=0$; (middle column)
$\displaystyle\textstyle g_{1b}/2\pi=10\,\mathrm{MHz}$; (right column)
$\displaystyle\textstyle g_{1b}/2\pi=25\,\mathrm{MHz}$. The rest of the
parameters are the same as in Fig. 4. For a better perception of the numerical
results, we supplement the plots for $\displaystyle\textstyle g_{1b}\neq 0$
(middle and right columns) with the corresponding curves (shown in lower
saturation) for the decoupled case $\displaystyle\textstyle g_{1b}=0$.

Using the time-dependent wavefunction of the system, one can compute various
observables such as populations of the transmon eigenlevels, the resonators,
and the waveguides at arbitrary moments of time. The population of
$\displaystyle|s\rangle_{j}$ eigenlevel of the $\displaystyle\textstyle j$-th
transmon is given by $\displaystyle\textstyle
P_{sj}(t)=\langle\varPsi(t)|\sigma^{ss}_{j}|\varPsi(t)\rangle$. With the
wavefunction expressed by Eq. (13), one obtains

$\begin{split}P_{e1}(t)&\,=|X(t)|^{2}+|Y_{b}(t)|^{2}+\int^{\infty}_{0}\mathrm{d}\nu|\varTheta^{\beta}_{\nu}(t)|^{2},\\\
P_{f1}(t)&\,=|Q(t)|^{2}\end{split}$ (17)

for the populations of $\displaystyle|e\rangle_{1}$ and
$\displaystyle|f\rangle_{1}$ eigenlevels of the first transmon, respectively.
The population of $\displaystyle\textstyle|e\rangle_{2}$ eigenlevel of the
second transmon is given by

$P_{e2}(t)=|X(t)|^{2}+|Y_{a}(t)|^{2}+\int^{\infty}_{0}\mathrm{d}\nu|\varTheta^{\alpha}_{\nu}(t)|^{2}.$
(18)

The populations of resonators $\displaystyle\textstyle A$ and
$\displaystyle\textstyle B$ are expressed in terms of the probability
amplitudes as

$\begin{split}P_{a}(t)&\,=\langle\varPsi(t)|a^{\dagger}a|\varPsi(t)\rangle\\\
&\,=|R(t)|^{2}+|Y_{a}(t)|^{2}+\int^{\infty}_{0}\mathrm{d}\nu|\varXi^{\beta}_{\nu}(t)|^{2},\end{split}$
(19a)

$\begin{split}P_{b}(t)=&\,\langle\varPsi(t)|b^{\dagger}b|\varPsi(t)\rangle\\\
&\,=|R(t)|^{2}+|Y_{b}(t)|^{2}+\int^{\infty}_{0}\mathrm{d}\nu|\varXi^{\alpha}_{\nu}(t)|^{2}.\end{split}$
(19b)

Finally, the waveguide populations are computed as

$\begin{split}P_{\alpha}(t)&\,=\int^{\infty}_{0}\mathrm{d}\nu\,\langle\varPsi(t)|\alpha^{\dagger}_{\nu}\alpha_{\nu}|\varPsi(t)\rangle\\\
&\,=p_{\alpha\beta}(t)+\int^{\infty}_{0}\mathrm{d}\nu\left[|\varXi^{\alpha}_{\nu}(t)|^{2}+|\varTheta^{\alpha}_{\nu}(t)|^{2}\right],\end{split}$
(20a)
$\begin{split}P_{\beta}(t)&\,=\int^{\infty}_{0}\mathrm{d}\nu\,\langle\varPsi(t)|\beta^{\dagger}_{\nu}\beta_{\nu}|\varPsi(t)\rangle\\\
&\,=p_{\alpha\beta}(t)+\int^{\infty}_{0}\mathrm{d}\nu\left[|\varXi^{\beta}_{\nu}(t)|^{2}+|\varTheta^{\beta}_{\nu}(t)|^{2}\right],\end{split}$
(20b)

where $\displaystyle\textstyle p_{\alpha\beta}(t)$ is given by Eq. (14).

Figure 5 demonstrates the dynamics of the populations of the transmons
eigenlevels as well as the populations of the resonators and the waveguides
for different values of coupling $\displaystyle\textstyle g_{1b}$ between
resonator $\displaystyle\textstyle B$ and
$\displaystyle\textstyle|e\rangle_{1}\leftrightarrow|f\rangle_{1}$ transition
of the first transmon. In agreement with the earlier qualitative
considerations and numerical results, we observe that switching on the
coupling between the first transmon and resonator $\displaystyle\textstyle B$
leads to pronounced modification of the emission dynamics compared to the
decoupled case.

Figure 6: Plots (a)–(f) show the dependence of the entanglement entropy
$\displaystyle\textstyle S_{\mathrm{ent}}$ on $\displaystyle\textstyle g_{1b}$
for different values of $\displaystyle\textstyle g_{1a}$ and
$\displaystyle\textstyle\kappa$ (encoded by the color gradient), where we set
$\displaystyle\textstyle(\kappa_{a}=\kappa_{b})=\kappa$. Parameters of the
system used for computations are the following: (a) $\displaystyle\textstyle
g_{1a}/2\pi=5\,\textrm{MHz}$, (b) $\displaystyle\textstyle
g_{1a}/2\pi=7\,\textrm{MHz}$, (c) $\displaystyle\textstyle
g_{1a}/2\pi=10\,\textrm{MHz}$, (d) $\displaystyle\textstyle
g_{1a}/2\pi=15\,\textrm{MHz}$, (e) $\displaystyle\textstyle
g_{1a}/2\pi=20\,\textrm{MHz}$, (f) $\displaystyle\textstyle
g_{1a}/2\pi=30\,\textrm{MHz}$. Plots (g) and (h) demonstrate the effect of
$\displaystyle\textstyle\kappa$ on the maximal value of the entanglement
entropy $\displaystyle\textstyle S^{\mathrm{max}}_{\mathrm{ent}}$, which can
be achieved for the given values of $\displaystyle\textstyle g_{1a}$ and
$\displaystyle\textstyle g_{2b}$. Parameters used in plot (g):
$\displaystyle\textstyle g_{1a}/2\pi=5\,\mathrm{MHz}$ (blue circles),
$\displaystyle\textstyle g_{1a}/2\pi=7\,\mathrm{MHz}$ (red triangles),
$\displaystyle\textstyle g_{1a}/2\pi=10\,\mathrm{MHz}$ (grey diamonds).
Parameters used in plot (h): $\displaystyle\textstyle
g_{1a}/2\pi=15\,\mathrm{MHz}$ (blue circles), $\displaystyle\textstyle
g_{1a}/2\pi=20\,\mathrm{MHz}$ (red triangles), $\displaystyle\textstyle
g_{1a}/2\pi=30\,\mathrm{MHz}$ (grey diamonds). For all plots, we set
$\displaystyle\textstyle g_{2b}/2\pi=10\,\textrm{MHz}$.

##  VI Photon entanglement

In the course of its evolution, the considered system eventually reaches the
state when all emission processes are finished, the transmons have relaxed to
their ground states, and the photons propagate in the waveguides as free
excitations. This final state of the system
$\displaystyle\textstyle|\Psi_{\mathrm{fi}}\rangle$ is given by

$|\Psi_{\mathrm{fi}}\rangle=\int\mathrm{d}\nu\int\mathrm{d}\nu^{\prime}\Phi^{\mathrm{out}}_{\nu,\nu^{\prime}}\alpha^{\dagger}_{\nu}\beta^{\dagger}_{\nu^{\prime}}|\varnothing\rangle,$
(21)

with $\displaystyle\textstyle\Phi^{\mathrm{out}}_{\nu,\nu^{\prime}}$ being the
joint spectral amplitude of the outgoing photon pair. If the emitted photons
are not entangled, their joint spectral amplitude can be factorized into a
product of single-photon amplitudes
$\displaystyle\textstyle\Phi^{\mathrm{out}}_{\nu,\nu^{\prime}}=\phi_{\alpha}(\nu)\phi_{\beta}(\nu^{\prime})$,
while the joint spectral amplitude of TF entangled photon pair is not
factorable
$\displaystyle\textstyle\Phi^{\mathrm{out}}_{\nu,\nu^{\prime}}\neq\phi_{\alpha}(\nu)\phi_{\beta}(\nu^{\prime})$.
For testing the factorizability of the joint spectral amplitude of the emitted
photons $\displaystyle\textstyle\Phi^{\mathrm{out}}_{\nu,\nu^{\prime}}$, we
perform its Schmidt decomposition [78, 79]:

$\Phi^{\mathrm{out}}_{\nu,\nu^{\prime}}=\sum_{j}\sqrt{\lambda_{j}}\,\varphi_{j,\nu}\vartheta_{j,\nu^{\prime}},$
(22)

where the weights $\displaystyle\textstyle\lambda_{j}\geq 0$ satisfying the
condition $\displaystyle\textstyle\sum_{j}\lambda_{j}=1$ are usually referred
to as the Schmidt coefficients, and the single-photon spectral amplitudes
$\displaystyle\textstyle\varphi_{j,\nu}$ and
$\displaystyle\textstyle\vartheta_{j,\nu}$ are called the Schmidt modes. The
latter constitute a complete set of orthonormal functions obeying the
relations
$\displaystyle\textstyle\sum_{j}\zeta^{*}_{j,\nu}\zeta_{j,\nu^{\prime}}=\delta(\nu-\nu^{\prime})$
and
$\displaystyle\textstyle\int\mathrm{d}\nu\,\zeta^{*}_{j,\nu}\zeta_{k,\nu}=\delta_{j,k}$,
where $\displaystyle\textstyle\zeta=\varphi$ or $\displaystyle\vartheta$.

The Schmidt coefficients $\displaystyle\textstyle\lambda_{j}$ are determined
by solution of the eigenvalue problem [78]:

$\begin{split}\int\mathrm{d}\nu^{\prime}\,\mathcal{K}^{\varphi}_{\nu,\nu^{\prime}}\,\varphi_{j,\nu^{\prime}}=\lambda_{j}\varphi_{j,\nu},\\\
\int\mathrm{d}\nu^{\prime}\,\mathcal{K}^{\vartheta}_{\nu,\nu^{\prime}}\,\vartheta_{j,\nu^{\prime}}=\lambda_{j}\vartheta_{j,\nu},\end{split}$
(23)

where the integral kernels
$\displaystyle\textstyle\mathcal{K}^{\varphi}_{\nu,\nu^{\prime}}$ and
$\displaystyle\textstyle\mathcal{K}^{\vartheta}_{\nu,\nu^{\prime}}$ are
defined as

$\begin{split}\mathcal{K}^{\varphi}_{\nu,\nu^{\prime}}=\int\mathrm{d}\varpi\,(\Phi^{\mathrm{out}}_{\nu,\varpi})^{*}\Phi^{\mathrm{out}}_{\nu^{\prime},\varpi},\\\
\mathcal{K}^{\vartheta}_{\nu,\nu^{\prime}}=\int\mathrm{d}\varpi\,(\Phi^{\mathrm{out}}_{\varpi,\nu})^{*}\Phi^{\mathrm{out}}_{\varpi,\nu^{\prime}}.\end{split}$
(24)

Figure 7: Dependence of $\displaystyle\textstyle
S^{\mathrm{max}}_{\mathrm{ent}}$ on $\displaystyle\textstyle\kappa_{b}$ for
different values of $\displaystyle\textstyle\kappa_{a}$:
$\displaystyle\textstyle\kappa_{a}/2\pi=10\,\mathrm{MHz}$ (blue circles),
$\displaystyle\textstyle\kappa_{a}/2\pi=25\,\mathrm{MHz}$ (red triangles), and
$\displaystyle\textstyle\kappa_{a}/2\pi=50\,\mathrm{MHz}$ (grey diamonds), and
$\displaystyle\textstyle g_{1a}$: (a) $\displaystyle\textstyle
g_{1a}/2\pi=5\,\mathrm{MHz}$, (b) $\displaystyle\textstyle
g_{1a}/2\pi=10\,\mathrm{MHz}$, (c) $\displaystyle\textstyle
g_{1a}/2\pi=20\,\mathrm{MHz}$. The rest of the system parameters are the same
as in Fig. 6.

The eigenvalue problem in Eq. (23) is treated numerically using the tools
provided by the Mathematica system. We transform the eigenvalue problem for
the integral operators into that for the matrices. For this purpose we
discretize the kernels
$\displaystyle\textstyle\mathcal{K}^{\varphi}_{\nu,\nu^{\prime}}$ and
$\displaystyle\textstyle\mathcal{K}^{\vartheta}_{\nu,\nu^{\prime}}$ into a
uniform $\displaystyle\textstyle\mathsf{N}\times\mathsf{N}$ grid on a square
domain spanning $\displaystyle\textstyle\pm\varDelta\omega$ around the
resonator frequencies $\displaystyle\textstyle\omega_{a}$ and
$\displaystyle\textstyle\omega_{b}$, where
$\displaystyle\textstyle\varDelta\omega$ satisfies the condition
$\displaystyle\textstyle\int^{\omega_{a}+\varDelta\omega}_{\omega_{a}-\varDelta\omega}\mathrm{d}\nu\int^{\omega_{b}+\varDelta\omega}_{\omega_{b}-\varDelta\omega}\mathrm{d}\nu^{\prime}|\Phi^{\mathrm{out}}_{\nu,\nu^{\prime}}|^{2}\geq\Theta$
with $\displaystyle\textstyle 0<\Theta<1$. The values of
$\displaystyle\textstyle\mathcal{K}^{\varphi}_{\nu,\nu^{\prime}}$ and
$\displaystyle\textstyle\mathcal{K}^{\vartheta}_{\nu,\nu^{\prime}}$ in the
grid nodes are computed by the numerical evaluation of the integrals in Eq.
(24) using the NIntegrate function. The joint spectral amplitude of the
outgoing photons is evaluated as
$\displaystyle\textstyle\Phi^{\mathrm{out}}_{\nu,\nu^{\prime}}=\varPhi_{\nu,\nu^{\prime}}(t_{\infty})$,
where the moment of time $\displaystyle\textstyle t_{\infty}$ is determined as
$\displaystyle\textstyle p_{\alpha\beta}(t_{\infty})>0.999$. Then, we use the
Eigenvalues function for determining the set of eigenvalues
$\displaystyle\textstyle\\{\Lambda_{j}\\}$ for the obtained matrices. A finite
discretization of the bounded domain of photon frequencies results in
$\displaystyle\textstyle(\mathcal{I}\equiv\sum^{\mathsf{N}}_{j=1}\Lambda_{j})<1$,
so we make a normalization
$\displaystyle\textstyle\\{\lambda_{j}\\}=\\{\Lambda_{j}\\}/\mathcal{I}$ to
ensure that $\displaystyle\textstyle\sum_{j}\lambda_{j}=1$. We set
$\displaystyle\textstyle\mathsf{N}=100$ and
$\displaystyle\textstyle\Theta=0.99$ for all computations. For the system
parameters we use for computations, the extension of the frequency domain (by
setting larger $\displaystyle\textstyle\Theta$) and using the finer grid (by
increasing $\displaystyle\textstyle\mathsf{N}$) has only a minor effect on the
evaluated values of the Schmidt coefficients and the entanglement entropies.

As a measure of entanglement of the emitted photons, we use the entanglement
(von Neumann) entropy $\displaystyle\textstyle S_{\mathrm{ent}}$, which is
expressed via the Schmidt coefficients as [80]:

$S_{\mathrm{ent}}=-\sum_{j}\lambda_{j}\log_{2}\lambda_{j}.$ (25)

The non-zero entanglement entropy, $\displaystyle\textstyle
S_{\mathrm{ent}}>0$, implies that emitted photons are entangled.

Figure 6 aggregates the results of computations demonstrating the dependence
of the entanglement entropy $\displaystyle\textstyle S_{\mathrm{ent}}$ of the
emitted photons on the interrelation between the transmon-resonator coupling
parameters and the photon leakage rate $\displaystyle\textstyle\kappa$ from
the resonators to the corresponding waveguides for the case
$\displaystyle\textstyle\kappa=(\kappa_{a}=\kappa_{b})$. Computations reveal
that for $\displaystyle\textstyle g_{1a}<g_{2b}$, the entanglement entropy
rapidly grows with the increase of $\displaystyle\textstyle g_{1b}$, reaching
its maximum $\displaystyle\textstyle S^{\mathrm{max}}_{\mathrm{ent}}$ for some
value of $\displaystyle\textstyle g_{1b}$, then slowly decreasing with the
further increase of $\displaystyle\textstyle g_{1b}$. In this regime, faster
photon leakage to the waveguides (i.e., shorter photon lifetimes inside the
resonators) leads to stronger photon entanglement, as illustrated in Figs.
6(g). The increase of the ratio $\displaystyle\textstyle g_{1a}/g_{2b}$
eventually leads to the opposite behavior, when the entanglement weakens with
the increase of the photon leakage rates, which is shown in Fig. 6(h). Figure
7 demonstrates the results of computations for
$\displaystyle\textstyle\kappa_{a}\neq\kappa_{b}$. Figure 7(a) shows that for
$\displaystyle\textstyle g_{1a}<g_{2b}$, stronger entanglement is achieved for
larger values of $\displaystyle\textstyle\kappa_{a}$. The crossover to the
regime of $\displaystyle\textstyle g_{1a}>g_{2b}$ results in the opposite
behavior of the entanglement entropy when stronger entanglement is achieved
for the lower photon leakage rate $\displaystyle\textstyle\kappa_{a}$, which
is shown in Fig. 7(c). Lower photon leakage rate
$\displaystyle\textstyle\kappa_{b}$ gives stronger entanglement, but the
dependence of the entanglement entropy on $\displaystyle\textstyle\kappa_{b}$
is rather weak, especially for $\displaystyle\textstyle g_{1a}>g_{2b}$. Thus,
the results of computations suggest that stronger entanglement of the emitted
photons is achieved by increasing the ratios $\displaystyle\textstyle
g_{2b}/g_{1a}$ and $\displaystyle\textstyle g_{1b}/g_{2b}$ in the regime of
$\displaystyle\textstyle\kappa_{a}\gtrsim g_{1a}$ and
$\displaystyle\textstyle\kappa_{b}\gtrsim g_{2b}$.

Let us elucidate some features of the obtained behavior of the entanglement
entropy using the general considerations discussed in Sec. II. The relaxation
of the first transmon not only delivers the photon via resonator
$\displaystyle\textstyle A$ to waveguide $\displaystyle\textstyle\alpha$ but
also triggers the relaxation of the second transmon leading to the emission of
the photon into waveguide $\displaystyle\textstyle\beta$. For gaining stronger
temporal correlations between the emitted photons, we need the photon to be
emitted into the waveguide $\displaystyle\textstyle\beta$ shortly after the
emission of the photon into waveguide $\displaystyle\textstyle\alpha$
triggered the relaxation of the second transmon. That is attained by
increasing the coupling between the second transmon and resonator
$\displaystyle\textstyle B$, i.e., by increasing the ratio
$\displaystyle\textstyle g_{2b}/g_{1a}$. The larger ratios
$\displaystyle\textstyle g_{1b}/g_{2b}$ are required for the efficient
inhibition of the second transmon relaxation until the moment the first
transmon decays, delivering a photon to waveguide
$\displaystyle\textstyle\alpha$. The mechanism of that inhibition is clarified
in the second paragraph of Sec. II.

##  VII Discussion and Summary

Having outlined the scheme and operational principle of the on-demand source
of microwave TF entangled photon pairs and investigated its performance, let
us now briefly discuss the general idea of how the proposed scheme can be
extended for the generation of multiphoton entangled states. For this purpose,
we consider the setup for the generation of the three-photon entangled states
illustrated in Fig. 8. Compared to the original setup for the generation of
photon pairs shown in Fig. 2, here we added one more resonator (marked as
$\displaystyle\textstyle C$) with frequency
$\displaystyle\textstyle\omega_{c}$ coupled to the output transmission line
and a third transmon coupled to resonators $\displaystyle\textstyle A$ and
$\displaystyle\textstyle C$. Following the logic of Sec. II, the frequencies
$\displaystyle\textstyle\omega^{ge}_{3}$ and
$\displaystyle\textstyle\omega^{ef}_{3}$ of
$\displaystyle|g\rangle_{3}\leftrightarrow|e\rangle_{3}$ and
$\displaystyle\textstyle|e\rangle_{3}\leftrightarrow|f\rangle_{3}$ transitions
of the third transmon are set that
$\displaystyle\textstyle\omega^{ge}_{3}=\omega_{c}$ and
$\displaystyle\textstyle\omega^{ef}_{3}=\omega_{a}$. The interrelation between
the frequencies of the resonators and the transition frequencies of the
transmons is schematically shown in Fig. 8(b). Such a choice of frequencies
enables the resonant excitation exchange between
$\displaystyle\textstyle|g\rangle_{3}\leftrightarrow|e\rangle_{3}$ transition
and resonator $\displaystyle\textstyle C$ and
$\displaystyle\textstyle|e\rangle_{3}\leftrightarrow|f\rangle_{3}$ transition
and resonator $\displaystyle\textstyle A$, while the excitation exchange
between $\displaystyle\textstyle|g\rangle_{3}\leftrightarrow|e\rangle_{3}$
transition and between resonator $\displaystyle\textstyle A$ and
$\displaystyle\textstyle|e\rangle_{3}\leftrightarrow|f\rangle_{3}$ transition
and resonator $\displaystyle\textstyle C$ is inhibited due to detuning. The
relaxation of the third transmon excited state
$\displaystyle\textstyle|e\rangle_{3}$ decouples
$\displaystyle\textstyle|e\rangle_{3}\leftrightarrow|f\rangle_{3}$ transition
from resonator $\displaystyle\textstyle A$, triggering the relaxation of the
first transmon excited state $\displaystyle\textstyle|e\rangle_{1}$, which, in
turn, triggers the relaxation of the second transmon. Thus, one may anticipate
that the photons are emitted into the corresponding waveguides in TF entangled
triples.

Figure 8: (a) The variant of extension of the proposed scheme of generation of
TF entangled photons: schematics of the potential circuit QED setup for
generation of three-photon entangled states. (b) Scheme illustrating the
relationships between the frequencies of the resonators and the transition
frequencies of the transmons in the setup shown in Fig. 8(a).

To summarize, we proposed a feasible scheme for efficient on-demand generation
of pairs of microwave TF entangled photons propagating in separate waveguides
and analyzed its performance. The design and working parameters of the
considered setup are readily accessible for the state-of-the-art
superconducting circuit QED technologies. Although we outlined a rather
concrete superconducting circuit QED implementation of an on-demand source of
microwave TF entangled photon pairs, we should stress that the general
principle of operation does not rely on the details of the circuit QED
realization and can be applied to other physical systems. We focused our
analysis on the pulsed on-demand regime of photon pairs generation. Therefore,
the consideration of the continuous regime of operation is of interest for
future work. Besides, a detailed investigation of the extended schemes (as
shown in Fig. 8) for the generation of multiphoton TF entangled states
constitutes a potential research direction as well.

###### Acknowledgements.

The author thanks Andrii Semenov and Andrii Sokolov for useful comments. This
work was supported by the National Academy of Sciences of Ukraine through the
Program of Postdoctoral Researches.

##  Appendix A Pair of 2LEs coupled to resonator

The paradigmatic system described in the first paragraph of Sec. II and
schematically shown in Fig. 1(a) is modeled by the Hamiltonian as follows

$\begin{split}{\mathcal{H}}=&\,\omega
c^{\dagger}c+\sum_{j=1}^{2}(\omega+\Delta_{j})\sigma^{+}_{j}\sigma^{-}_{j}+\sum_{j=1}^{2}g_{j}(c^{\dagger}\sigma^{-}_{j}+\sigma^{+}_{j}c)\\\
&\,+\int^{\infty}_{0}\mathrm{d}\nu\nu
A^{\dagger}_{\nu}A_{\nu}+\int\mathrm{d}\nu
f(\nu)(A^{\dagger}_{\nu}c+c^{\dagger}A_{\nu}).\end{split}$ (26)

The first three terms in the above Hamiltonian describe the single-mode
resonator with frequency $\displaystyle\textstyle\omega$ coupled to a pair of
2LEs, where $\displaystyle\textstyle\Delta_{j}$ stands for the detuning
between the frequencies of the $\displaystyle\textstyle j$-th 2LE and the
resonator, parameter $\displaystyle\textstyle g_{j}$ is the coupling strength
of the resonator to the $\displaystyle\textstyle j$-th 2LE
($\displaystyle\textstyle j\in\\{1,2\\}$). Operator $\displaystyle\textstyle
c$ ($\displaystyle\textstyle c^{\dagger}$) annihilates (creates) a photon in
the resonator, and $\displaystyle\textstyle\sigma^{+}_{j}$
($\displaystyle\textstyle\sigma^{-}_{j}$) rises (lowers) the state of the
$\displaystyle\textstyle j$-th 2LE. The last pair of terms in the Hamiltonian
in Eq. (26) describes the waveguide, represented by a bath of independent
bosonic modes, and its coupling to the resonator with strength
$\displaystyle\textstyle f(\nu)$. Operator $\displaystyle\textstyle A_{\nu}$
($\displaystyle\textstyle A^{\dagger}_{\nu}$) annihilates (creates) a photon
with frequency $\displaystyle\textstyle\nu$ propagating in the waveguide.

Now, let us demonstrate that in the single-excitation case, the paradigmatic
system composed of a pair of 2LEs coupled to the resonator, considered in Sec.
II and illustrated in Fig. 1(a), can be equivalently represented as a 2LE
coupled to the V-configuration 3LE. To proceed, we formally represent the
Hamiltonian $\displaystyle\textstyle{\mathcal{H}}$ given by Eq. (26) in the
eigenbasis of the 2LE-resonator Hamiltonian
$\displaystyle\textstyle{\mathcal{H}}_{\mathrm{JC}}$. The latter is
constituted by the terms describing the resonator, the first 2LE, and their
interaction in Eq. (26), which together form the Hamiltonian of the Jaynes-
Cummings (JC) model [81]. The JC system ground state
$\displaystyle|\mathcal{G}\rangle=|0\rangle_{\mathrm{r}}|g\rangle_{1}$ is a
state with a vacuum field in the resonator
$\displaystyle\textstyle|0\rangle_{\mathrm{r}}$ and the first 2LE in the
ground state $\displaystyle\textstyle|g\rangle_{1}$. The excited eigenstates
of $\displaystyle{\mathcal{H}}_{\mathrm{JC}}$ are the superpositions of the
state $\displaystyle\textstyle|n\rangle_{\mathrm{r}}|g\rangle_{1}$ containing
$\displaystyle\textstyle n$ photons in the resonator with the ground-state 2LE
and the state $\displaystyle\textstyle|n-1\rangle_{\mathrm{r}}|e\rangle_{1}$
with $\displaystyle\textstyle n-1$ photons in the resonator and the excited-
state 2LE [81, 68]:

$\left(\begin{array}[]{c}|\mathcal{E}_{n}^{-}\rangle\\\
|\mathcal{E}_{n}^{+}\rangle\end{array}\right)=\left(\begin{array}[]{cc}\cos\mu_{n}&-\sin\mu_{n}\\\
\sin\mu_{n}&\cos\mu_{n}\end{array}\right)\left(\begin{array}[]{c}|n\rangle_{\mathrm{r}}|g\rangle_{1}\\\
|n-1\rangle_{\mathrm{r}}|e\rangle_{1}\end{array}\right),$

where

$\tan\mu_{n}=\sqrt{\frac{\Lambda_{n}-\Delta_{1}}{\Lambda_{n}+\Delta_{1}}},\quad\Lambda_{n}=\sqrt{4ng^{2}+\Delta_{1}^{2}}.$

The JC eigenstates $\displaystyle\textstyle|\mathcal{E}^{\pm}_{n}\rangle$
correspond to the eigenfrequencies $\displaystyle\textstyle
E^{\pm}_{n}=n\omega+(\Delta_{1}\pm\Lambda_{n})/2$.

Restricting ourselves to the resonant regime of the 2LE-resonator coupling
($\displaystyle\textstyle\Delta_{1}=\Delta_{2}=0$) and the single-excitation
domain, as considered in Sec. II, we project the system Hamiltonian
$\displaystyle\textstyle{\mathcal{H}}$ on the lowest eigenstates
$\displaystyle\textstyle|\mathcal{G}\rangle$ and
$\displaystyle\textstyle|\mathcal{E}^{\pm}_{1}\rangle$ of the JC Hamiltonian.
Thus, one has
$\displaystyle\textstyle{\mathcal{H}}\rightarrow\varPi{\mathcal{H}}\varPi$,
where
$\displaystyle\textstyle\varPi=|\mathcal{G}\rangle\langle\mathcal{G}|+\sum_{\pm}|\mathcal{E}^{\pm}_{1}\rangle\langle\mathcal{E}^{\pm}_{1}|$,
which gives:

$\begin{split}{\mathcal{H}}=&\,\omega\sigma^{+}_{2}\sigma^{-}_{2}+\sum_{\pm}(\omega\pm
g_{1})|\mathcal{E}^{\pm}\rangle\langle\mathcal{E}^{\pm}|+\int^{\infty}_{0}\mathrm{d}\nu\nu
A^{\dagger}_{\nu}A_{\nu}\\\
&\,+\frac{g_{2}}{\sqrt{2}}\sum_{\pm}\left(\sigma^{+}_{2}|\mathcal{G}\rangle\langle\mathcal{E}^{\pm}|+|\mathcal{E}^{\pm}\rangle\langle\mathcal{G}|\sigma^{-}_{2}\right)\\\
&\,+\int^{\infty}_{0}\mathrm{d}\nu\,\frac{f(\nu)}{\sqrt{2}}\sum_{\pm}(A^{\dagger}_{\nu}|\mathcal{G}\rangle\langle\mathcal{E}^{\pm}|+|\mathcal{E}^{\pm}\rangle\langle\mathcal{G}|A_{\nu}),\end{split}$

where for brevity we dropped the subscripts indicating the photon number in
the notations of JC eigenstates. The above Hamiltonian describes the 3LE
coupled to the 2LE and the waveguide. The 3LE is constituted by the ground
state $\displaystyle|\mathcal{G}\rangle$ and a pair of excited states
$\displaystyle|\mathcal{E}^{\pm}\rangle$. The frequency of
$\displaystyle\textstyle|\mathcal{G}\rangle\leftrightarrow|\mathcal{E}^{\pm}\rangle$
transition is $\displaystyle\omega\pm g_{1}$. The levels of the 3LE are
arranged in the $\displaystyle\mathsf{V}$-type configuration. Each transition
of this $\displaystyle\mathsf{V}$3LE is coupled to the second 2LE with
strength $\displaystyle g_{2}/\sqrt{2}$ and to the waveguide with strength
$\displaystyle f(\nu)/\sqrt{2}$. This equivalent representation of the system
composed of the resonator resonantly coupled to a pair of 2LEs is illustrated
in Fig. 1(b).

### Dynamics of the resonator-2LEs system

Figure 9: Effect of the 2LEs couplings ratio $\displaystyle\textstyle
g_{1}/g_{2}$ on the evolution of the waveguide population
$\displaystyle\textstyle P_{\mathrm{w}}(t)$ in the system shown in Fig. 1(a).
The parameters used for calculations are as follows: $\displaystyle\textstyle
g_{2}/2\pi=10\,\mathrm{MHz}$ and
$\displaystyle\textstyle\kappa/2\pi=2\,\mathrm{MHz}$.

The waveguide population at time $\displaystyle\textstyle t$ is given by:

$P_{\mathrm{w}}(t)=\kappa\int^{t}_{0}\mathrm{d}\tau\langle
c^{\dagger}(\tau)c(\tau)\rangle,$ (27)

with $\displaystyle\textstyle\kappa=2\pi f^{2}(\omega)$ being the resonator
decay rate into the waveguide. The averaging in Eq. (27) goes over the initial
state of the system
$\displaystyle\textstyle|\psi_{0}\rangle=\sigma^{+}_{2}|0\rangle$ with
$\displaystyle\textstyle|0\rangle=|0\rangle_{\mathrm{r}}|g\rangle_{1}|g\rangle_{2}$
being the vacuum state of the resonator-2LEs system. In the single-excitation
case, the resonator population $\displaystyle\textstyle\langle
c^{\dagger}(t)c(t)\rangle$ can be evaluated using the relation
$\displaystyle\textstyle\langle
c^{\dagger}(t)c(t)\rangle=\langle\psi_{0}|c^{\dagger}(t)|0\rangle\langle
0|c(t)|\psi_{0}\rangle$. The evolution of the matrix element
$\displaystyle\textstyle\langle 0|c(t)|\psi_{0}\rangle$ is governed by the set
of equations of motion:

$\frac{\partial}{\partial t}\\!\left(\\!\begin{array}[]{c}\langle
0|c(t)|\psi_{0}\rangle\\\ \langle 0|\sigma^{-}_{1}(t)|\psi_{0}\rangle\\\
\langle
0|\sigma^{-}_{2}(t)|\psi_{0}\rangle\end{array}\\!\right)=-\mathrm{i}\begin{pmatrix}\widetilde{\omega}&g_{1}&g_{2}\\\
g_{1}&\omega&0\\\
g_{2}&0&\omega\end{pmatrix}\\!\\!\left(\\!\begin{array}[]{c}\langle
0|c(t)|\psi_{0}\rangle\\\ \langle 0|\sigma^{-}_{1}(t)|\psi_{0}\rangle\\\
\langle 0|\sigma^{-}_{2}(t)|\psi_{0}\rangle\end{array}\\!\right),$

where $\displaystyle\textstyle\widetilde{\omega}=\omega-\mathrm{i}\kappa/2$.

Using the solution of the above set of equations in Eq. (27), one obtains

$P_{\mathrm{w}}(t)=\frac{g_{2}^{2}}{g_{1}^{2}+g_{2}^{2}}\Bigg{[}1-F(t)\exp\left(-\frac{\kappa}{2}t\right)\Bigg{]},$
(28)

with $\displaystyle\textstyle F(t)$ given by

$F(t)=1+\kappa\sin(\varOmega t)\,\frac{\kappa\sin(\varOmega
t)+4\varOmega\cos(\varOmega t)}{8\varOmega^{2}},$

where
$\displaystyle\textstyle\varOmega=\sqrt{g_{1}^{2}+g_{2}^{2}-(\kappa/4)^{2}}$.

Figure 9 shows the dependence of the dynamics of the waveguide population
$\displaystyle\textstyle P_{\mathrm{w}}(t)$ on the ratio of the couplings
$\displaystyle\textstyle g_{1}/g_{2}$. Calculations demonstrate that the
increase of $\displaystyle\textstyle g_{1}/g_{2}$ leads to the inhibition of
the photon leakage into the waveguide. It follows from Eq. (28) that
$\displaystyle\textstyle P_{\mathrm{w}}(t\rightarrow\infty)<1$ for
$\displaystyle\textstyle g_{1}\neq 0$. For $\displaystyle\textstyle g_{1}\gg
g_{2}$, one has $\displaystyle\textstyle P_{\mathrm{w}}(t\rightarrow\infty)\ll
1$ implying that the leakage into the waveguide is suppressed and the photon
remains trapped within the resonator-2LEs system.

##  Appendix B Derivation of the effective Hamiltonian

###  B.1 Operators in the dressed basis

Using the definition of the unitary operator
$\displaystyle\textstyle{\mathcal{U}}$ given by Eq. (8) along with the Baker-
Campbell-Hausdorff formula, one obtains the following expressions for the
resonator photon annihilation operators in the dressed basis:

${\mathcal{U}}^{\dagger}a\,{\mathcal{U}}\approx\left(1+\frac{\lambda_{1a}^{2}}{2}\mathcal{Z}^{fe}_{1}\right)a+\lambda_{1a}\sigma_{1}^{ef}-\frac{\lambda_{1a}\lambda_{1b}}{2}b^{\dagger}\sigma_{1}^{gf},$
(29a)
${\mathcal{U}}^{\dagger}b\,{\mathcal{U}}\approx\left(1+\frac{\lambda_{1b}^{2}}{2}\mathcal{Z}^{eg}_{1}+\frac{\lambda_{2b}^{2}}{2}\mathcal{Z}^{fe}_{2}\right)b+\lambda_{1b}\sigma_{1}^{ge}+\lambda_{2b}\sigma_{2}^{ef}+\frac{\lambda_{1a}\lambda_{1b}}{2}a^{\dagger}\sigma_{1}^{gf},$
(29b)

where we keep the terms up to the second order in the small parameters
$\displaystyle\lambda_{1a}$, $\displaystyle\lambda_{1b}$, and
$\displaystyle\lambda_{2b}$. The resonator photon number operators
$\displaystyle\textstyle{\mathcal{N}}_{a}$ and
$\displaystyle\textstyle{\mathcal{N}}_{b}$ in the new basis read as

${\mathcal{U}}^{\dagger}{\mathcal{N}}_{a}\,{\mathcal{U}}\approx\left(1+\lambda_{1a}^{2}\mathcal{Z}^{fe}_{1}\right){\mathcal{N}}_{a}+\lambda_{1a}\left(a^{\dagger}\sigma_{1}^{ef}+\sigma_{1}^{fe}a\right)+\lambda^{2}_{1a}\sigma_{1}^{ff}-\frac{\lambda_{1a}\lambda_{1b}}{2}\left(a^{\dagger}b^{\dagger}\sigma_{1}^{gf}+\sigma_{1}^{fg}ab\right),$
(30a)
$\begin{split}{\mathcal{U}}^{\dagger}{\mathcal{N}}_{b}\,{\mathcal{U}}\approx&\,\left(1+\lambda_{1b}^{2}\mathcal{Z}^{eg}_{1}+\lambda_{2b}^{2}\mathcal{Z}^{fe}_{2}\right){\mathcal{N}}_{b}+\lambda_{1b}\left(b^{\dagger}\sigma_{1}^{ge}+\sigma_{1}^{eg}b\right)+\lambda_{2b}\left(b^{\dagger}\sigma_{2}^{ef}+\sigma_{2}^{fe}b\right)\\\
&\,+\lambda_{1b}^{2}\sigma_{1}^{ee}+\lambda_{2b}^{2}\sigma_{2}^{ff}+\frac{\lambda_{1a}\lambda_{1b}}{2}\left(a^{\dagger}b^{\dagger}\sigma_{1}^{gf}+\sigma_{1}^{fg}ab\right)+\lambda_{1b}\lambda_{2b}\left(\sigma_{2}^{fe}\sigma_{1}^{ge}+\sigma_{1}^{eg}\sigma_{2}^{ef}\right).\end{split}$
(30b)

In the dressed basis, the transmon ladder operators
$\displaystyle\textstyle\sigma^{ge}_{1,2}$ and
$\displaystyle\textstyle\sigma^{ef}_{1,2}$ acquire the form as follows

$\begin{split}{\mathcal{U}}^{\dagger}\sigma_{1}^{ge}{\mathcal{U}}\approx&\left[1-\frac{\lambda_{1a}^{2}}{2}{\mathcal{N}}_{a}-\frac{\lambda_{1b}^{2}}{2}(2{\mathcal{N}}_{b}+1)\right]\sigma_{1}^{ge}+\lambda_{1b}\mathcal{Z}^{eg}_{1}b+\lambda_{1a}a^{\dagger}\sigma_{1}^{gf}\\\
&\,-\lambda^{2}_{1b}\sigma_{1}^{eg}b^{2}+\lambda_{1a}\lambda_{1b}a^{\dagger}b\sigma_{1}^{ef}+\frac{\lambda_{1a}\lambda_{1b}}{2}\sigma_{1}^{fe}ab+\frac{\lambda_{1b}\lambda_{2b}}{2}\mathcal{Z}^{eg}_{1}\sigma_{2}^{ef},\end{split}$
(31a)
$\begin{split}{\mathcal{U}}^{\dagger}\sigma_{1}^{ef}{\mathcal{U}}\approx&\left[1-\frac{\lambda_{1a}^{2}}{2}(2{\mathcal{N}}_{a}+1)-\frac{\lambda^{2}_{1b}}{2}({\mathcal{N}}_{b}+1)\right]\sigma_{1}^{ef}+\lambda_{1a}\mathcal{Z}^{fe}_{1}a-\lambda_{1b}b^{\dagger}\sigma_{1}^{gf}\\\
&\,+\lambda_{1a}\lambda_{1b}\,b^{\dagger}a\,\sigma_{1}^{ge}+\frac{\lambda_{1a}\lambda_{1b}}{2}\,ab\,\sigma_{1}^{eg}-\lambda_{1a}^{2}\sigma_{1}^{fe}a^{2}-\frac{\lambda_{1b}\lambda_{2b}}{2}\sigma_{1}^{gf}\sigma_{2}^{fe},\end{split}$
(31b)
${\mathcal{U}}^{\dagger}\sigma_{2}^{ge}{\mathcal{U}}\approx\left(1-\frac{\lambda_{2b}^{2}}{2}{\mathcal{N}}_{b}\right)\sigma_{2}^{ge}+\lambda_{2b}b^{\dagger}\sigma_{2}^{gf}+\frac{\lambda_{1b}\lambda_{2b}}{2}\sigma_{1}^{eg}\sigma_{2}^{gf}$
(31c)
${\mathcal{U}}^{\dagger}\sigma_{2}^{ef}{\mathcal{U}}\approx\left[1-\frac{\lambda_{2b}^{2}}{2}(2{\mathcal{N}}_{b}+1)\right]\sigma_{2}^{ef}+\lambda_{2b}\mathcal{Z}^{fe}_{2}b-\lambda_{2b}^{2}\sigma_{2}^{fe}b^{2}+\frac{\lambda_{1b}\lambda_{2b}}{2}\mathcal{Z}^{fe}_{2}\sigma_{1}^{ge}$
(31d)

For the projection operators on the transmons eigenstates
$\displaystyle\textstyle\sigma_{1,2}^{ee}$ and
$\displaystyle\textstyle\sigma_{1,2}^{ff}$, one has:

$\begin{split}{\mathcal{U}}^{\dagger}\sigma_{1}^{ee}{\mathcal{U}}\approx&\left(1-\lambda_{1b}^{2}\right)\sigma_{1}^{ee}+\lambda_{1a}\left(a^{\dagger}\sigma_{1}^{ef}+\sigma_{1}^{fe}a\right)-\lambda_{1b}\left(b^{\dagger}\sigma_{1}^{ge}+\sigma_{1}^{eg}b\right)+\lambda_{1a}^{2}\sigma_{1}^{ff}\\\
&\,+\lambda_{1a}^{2}\mathcal{Z}^{fe}_{1}{\mathcal{N}}_{a}-\lambda_{1b}^{2}\mathcal{Z}^{eg}_{1}{\mathcal{N}}_{b}-\lambda_{1a}\lambda_{1b}\left(a^{\dagger}b^{\dagger}\sigma_{1}^{gf}+\sigma_{1}^{fg}ab\right)-\frac{\lambda_{1b}\lambda_{2b}}{2}\left(\sigma_{2}^{fe}\sigma_{1}^{ge}+\sigma_{1}^{eg}\sigma_{2}^{ef}\right),\end{split}$
(32a)
$\begin{split}{\mathcal{U}}^{\dagger}\sigma_{1}^{ff}{\mathcal{U}}\approx\left(1-\lambda_{1a}^{2}\right)\sigma_{1}^{ff}-\lambda_{1a}\left(a^{\dagger}\sigma_{1}^{ef}+\sigma_{1}^{fe}a\right)-\lambda_{1a}^{2}\mathcal{Z}^{fe}_{1}{\mathcal{N}}_{a}+\frac{\lambda_{1a}\lambda_{1b}}{2}\left(a^{\dagger}b^{\dagger}\sigma_{1}^{gf}+\sigma_{1}^{fg}ab\right),\end{split}$
(32b)
${\mathcal{U}}^{\dagger}\sigma_{2}^{ee}{\mathcal{U}}\approx\sigma_{2}^{ee}+\lambda_{2b}\left[b^{\dagger}\sigma_{2}^{ef}+\sigma_{2}^{fe}b\right]+\lambda_{2b}^{2}\mathcal{Z}^{fe}_{2}{\mathcal{N}}_{b}+\lambda_{2b}^{2}\sigma_{2}^{ff}+\frac{\lambda_{1b}\lambda_{2b}}{2}\left(\sigma_{2}^{fe}\sigma_{1}^{ge}+\sigma_{1}^{eg}\sigma_{2}^{ef}\right),$
(32c)
${\mathcal{U}}^{\dagger}\sigma_{2}^{ff}{\mathcal{U}}\approx(1-\lambda_{2b}^{2})\sigma_{2}^{ff}-\lambda_{2b}\left(\sigma_{2}^{fe}b+b^{\dagger}\sigma_{2}^{ef}\right)-\lambda_{2b}^{2}\mathcal{Z}^{fe}_{2}{\mathcal{N}}_{b}-\frac{\lambda_{1b}\lambda_{2b}}{2}\left(\sigma_{2}^{fe}\sigma_{1}^{ge}+\sigma_{1}^{eg}\sigma_{2}^{ef}\right).$
(32d)

Using Eqs. (29)–(32) and keeping the terms up to the first order in small
parameters $\displaystyle\textstyle\lambda_{1a}$,
$\displaystyle\textstyle\lambda_{1b}$, and
$\displaystyle\textstyle\lambda_{2b}$, one arrives at the Hamiltonian given by
Eq. (9).

###  B.2 Role of Purcell filters

Note that in the dressed basis the Hamiltonian of the resonator-waveguide
couplings $\displaystyle\textstyle{\mathcal{H}}_{\mathrm{r-w}}$ acquires the
form:
$\displaystyle\textstyle{\mathcal{U}}^{\dagger}{\mathcal{H}}_{\mathrm{r}-\mathrm{w}}{\mathcal{U}}\approx{\mathcal{H}}_{\mathrm{r}-\mathrm{w}}+{\mathcal{H}}_{\mathrm{t}-\mathrm{w}}$,
where $\displaystyle\textstyle{\mathcal{H}}_{\mathrm{r}-\mathrm{w}}$ is given
by Eq. (7), and the term
$\displaystyle\textstyle{\mathcal{H}}_{\mathrm{t}-\mathrm{w}}$ reads as

${\mathcal{H}}_{\mathrm{t}-\mathrm{w}}=\lambda_{1a}\int^{\infty}_{0}\mathrm{d}\nu
f_{a}(\nu)\big{(}\alpha^{\dagger}_{\nu}\sigma_{1}^{ef}+\sigma_{1}^{fe}\alpha_{\nu}\big{)}+\lambda_{1b}\int^{\infty}_{0}\mathrm{d}\nu
f_{b}(\nu)\big{(}\beta^{\dagger}_{\nu}\sigma_{1}^{ge}+\sigma_{1}^{eg}\beta_{\nu}\big{)}+\lambda_{2b}\int^{\infty}_{0}\mathrm{d}\nu
f_{b}(\nu)\big{(}\beta^{\dagger}_{\nu}\sigma_{2}^{ef}+\sigma_{2}^{fe}\beta_{\nu}\big{)}.$
(33)

The first pair of terms in
$\displaystyle\textstyle{\mathcal{H}}_{\mathrm{t}-\mathrm{w}}$ describes the
coupling between
$\displaystyle\textstyle|e\rangle_{1}\leftrightarrow|f\rangle_{1}$ and
$\displaystyle\textstyle|g\rangle_{1}\leftrightarrow|e\rangle_{1}$ transitions
of the first transmon and waveguides $\displaystyle\textstyle\alpha$ and
$\displaystyle\textstyle\beta$, respectively. The third term in Eq. (33)
describes the direct coupling between
$\displaystyle\textstyle|e\rangle_{2}\leftrightarrow|f\rangle_{2}$ transition
of the second transmon and waveguide $\displaystyle\textstyle\beta$ leading to
the relaxation of the state $\displaystyle\textstyle|f\rangle_{2}$. These
direct couplings of the transmon transitions to the waveguides lead to that
both excitations can be emitted into one waveguide instead of being emitted
into the separate waveguides. For mitigating this superfluous emission
processes, we use the Purcell filters. In the proposed setup (see Fig. 2), the
latter are represented by the resonators side-coupled to the output
waveguides. The frequencies of these resonators correspond to the dressed
frequencies $\displaystyle\textstyle\bar{\omega}_{1}^{ge}$ and
$\displaystyle\textstyle\bar{\omega}_{1}^{ef}$ of
$\displaystyle\textstyle|g\rangle_{1}\leftrightarrow|e\rangle_{1}$ and
$\displaystyle\textstyle|e\rangle_{1}\leftrightarrow|f\rangle_{1}$
transitions. The filter resonators reject the radiation in the narrow bands of
frequencies in the vicinity of $\displaystyle\textstyle\bar{\omega}_{1}^{ge}$
and $\displaystyle\textstyle\bar{\omega}_{1}^{ef}$ and transmit the radiation
with frequencies outside these bands [63]. Thus, the Purcell filters inhibit
the unwanted emission channels of the first transmon described by the first
two terms in the Hamiltonian (33). As we demonstrate in Sec. V, no more than
one photon could reside in either of the resonators, so
$\displaystyle\textstyle|e\rangle_{2}\rightarrow|f\rangle_{2}$ transition of
the second transmon is not involved in the dynamics of the system, and
$\displaystyle\textstyle|f\rangle_{2}$ level is not excited. Thus, the Purcell
filter cutting the emission from this transition into waveguide
$\displaystyle\textstyle\beta$ is not required. Since the processes of direct
relaxation of transmons into the waveguides are essentially suppressed by the
Purcell filters and occur on the timescales much longer than the photon
emission times, for simplicity, in the Hamiltonian
$\displaystyle\textstyle\mathcal{H}_{\mathrm{eff}}$, we drop the term
$\displaystyle\textstyle{\mathcal{H}}_{\mathrm{t}-\mathrm{w}}$ describing
these processes.

##  Appendix C Derivation of evolution equations

For deriving Eq. (15) governing the evolution of the probability amplitudes,
we follow the lines of Appendix B in Ref. [76]. For this purpose, we start
with the derivation of the equations of motion for the operators.

The effective Hamiltonian (9) generates the equations of motion for the
waveguide variables $\displaystyle\textstyle\alpha_{\nu}(t)$ and
$\displaystyle\textstyle\beta_{\nu}(t)$ as follows:

$\frac{\partial}{\partial
t}\alpha_{\nu}(t)=-\mathrm{i}\nu\alpha_{\nu}(t)-\mathrm{i}f_{a}(\nu)a(t),$
(34a) $\frac{\partial}{\partial
t}\beta_{\nu}(t)=-\mathrm{i}\nu\beta_{\nu}(t)-\mathrm{i}f_{b}(\nu)b(t),$ (34b)

with the formal solutions written as

$\alpha_{\nu}(t)=\tilde{\alpha}_{\nu}(t)-\mathrm{i}f_{a}(\nu)\int^{t}_{0}\mathrm{d}\tau\mathrm{e}^{-\mathrm{i}\nu(t-\tau)}a(\tau),$
(35a)
$\beta_{\nu}(t)=\tilde{\beta}_{\nu}(t)-\mathrm{i}f_{b}(\nu)\int^{t}_{0}\mathrm{d}\tau\mathrm{e}^{-\mathrm{i}\nu(t-\tau)}b(\tau),$
(35b)

where
$\displaystyle\textstyle\tilde{\alpha}_{\nu}(t)=\alpha_{\nu}(0)\mathrm{e}^{-\mathrm{i}\nu
t}$ and
$\displaystyle\textstyle\tilde{\beta}_{\nu}(t)=\beta_{\nu}(0)\mathrm{e}^{-\mathrm{i}\nu
t}$ stand for the annihilation operators of a photon propagating as a free
excitation in waveguide $\displaystyle\textstyle\alpha$ and
$\displaystyle\textstyle\beta$, respectively.

The equation of motion for the annihilation operator $\displaystyle\textstyle
a(t)$ of the photon in resonator $\displaystyle\textstyle A$ reads as

$\frac{\partial}{\partial
t}a(t)=-\mathrm{i}\left[\omega_{a}+\chi_{1a}\mathcal{Z}_{1}^{fe}(t)\right]a(t)-\mathrm{i}g_{1a}\sigma_{1}^{ge}(t)-\mathrm{i}\varUpsilon
b^{\dagger}(t)\sigma_{1}^{gf}(t)-\mathrm{i}\int^{\infty}_{0}\mathrm{d}\nu
f_{a}(\nu)\alpha_{\nu}(t).$ (36)

Plugging Eq. (35a) into the last term on the right-hand side of the above
equation gives

$\begin{split}\int^{\infty}_{0}\mathrm{d}\nu
f_{a}(\nu)\alpha_{\nu}(t)&=\alpha_{\mathrm{in}}(t)-\mathrm{i}\int^{t}_{0}\mathrm{d}\tau\int^{\infty}_{0}\mathrm{d}\nu
f^{2}_{a}(\nu)\mathrm{e}^{-\mathrm{i}\nu(t-\tau)}a(\tau),\end{split}$ (37)

where
$\displaystyle\textstyle\alpha_{\mathrm{in}}(t)=\int^{\infty}_{0}\mathrm{d}\nu\,f_{a}(\nu)\widetilde{\alpha}_{\nu}(t)$.
Using Eq. (36) and that $\displaystyle\textstyle|\chi_{1a}|\ll\omega_{a}$, one
can represent the operator $\displaystyle\textstyle a(t)$ as
$\displaystyle\textstyle
a(t)=\mathfrak{a}(t)\mathrm{e}^{-\mathrm{i}\omega_{a}t}$, where
$\displaystyle\textstyle\mathfrak{a}(t)$ stands for the slowly-varying part of
$\displaystyle a(t)$. With such a representation, one can notice that only the
frequencies in the vicinity of $\displaystyle\textstyle\omega_{a}$ contribute
significantly to the integral over $\displaystyle\textstyle\tau$ on the right-
hand side of Eq. (37). Thus, one can extend the lower bound of integration
over $\displaystyle\textstyle\nu$ to $\displaystyle\textstyle-\infty$, and
neglect the frequency dependence of the coupling parameter
$\displaystyle\textstyle f_{a}(\nu)$ assuming that $\displaystyle\textstyle
f_{a}(\nu)\approx f_{a}(\omega_{a})$. Using these approximations, one obtains

$\int^{\infty}_{0}\mathrm{d}\nu
f_{a}(\nu)\alpha_{\nu}(t)\approx\alpha_{\mathrm{in}}(t)-\mathrm{i}\frac{\kappa_{a}}{2}a(t),\quad\kappa_{a}=2\pi
f^{2}_{a}(\omega_{a}).$ (38)

Substituting this result into Eq. (36) yields

$\frac{\partial}{\partial
t}a(t)=-\mathrm{i}\left[\omega_{a}-\mathrm{i}\frac{\kappa_{a}}{2}+\chi_{1a}\mathcal{Z}_{1}^{fe}(t)\right]a(t)-\mathrm{i}g_{1a}\sigma_{1}^{ge}(t)-\mathrm{i}\varUpsilon
b^{\dagger}\sigma_{1}^{gf}(t)-\mathrm{i}\alpha_{\mathrm{in}}(t).$ (39)

The equation of motion for the operator $\displaystyle b(t)$ annihilating a
photon in the resonator $\displaystyle B$ is obtained similarly to Eq. (39)
and reads as

$\frac{\partial}{\partial
t}b(t)=-\mathrm{i}\left[\omega_{b}-\mathrm{i}\frac{\kappa_{b}}{2}+\chi_{1b}\mathcal{Z}_{1}^{eg}(t)+\chi_{2b}\mathcal{Z}_{2}^{fe}(t)\right]b(t)-\mathrm{i}g_{1b}\sigma_{1}^{ef}(t)-\mathrm{i}g_{2b}\sigma_{2}^{ge}(t)-\mathrm{i}\varUpsilon
a^{\dagger}(t)\sigma_{1}^{gf}(t)-\mathrm{i}\beta_{\mathrm{in}}(t),$ (40)

where
$\displaystyle\textstyle\beta_{\mathrm{in}}(t)=\int^{\infty}_{0}\mathrm{d}\nu\,f_{b}(\nu)\tilde{\beta}_{\nu}(t)$.
For the derivation of Eq. (40), we used the approximate relation

$\int^{\infty}_{0}\mathrm{d}\nu
f_{b}(\nu)\beta_{\nu}(t)\approx\beta_{\mathrm{in}}(t)-\mathrm{i}\frac{\kappa_{b}}{2}b(t),\quad\kappa_{b}=2\pi
f^{2}_{b}(\omega_{b}),$ (41)

which is obtained by analogy to Eq. (38).

The equations of motion for the required ladder operators of the transmons
read as follows:

$\frac{\partial}{\partial
t}\sigma_{1}^{ge}(t)=-\mathrm{i}\left[\omega_{1}^{ge}+\chi_{1b}-\chi_{1a}{\mathcal{N}}_{a}(t)+2\chi_{1b}{\mathcal{N}}_{b}(t)\right]\sigma_{1}^{ge}(t)+\mathrm{i}g_{1a}\mathcal{Z}_{1}^{eg}(t)a(t)-\mathrm{i}g_{1b}b^{\dagger}(t)\sigma_{1}^{gf}(t)+\mathrm{i}\varUpsilon\sigma_{1}^{fe}(t)a(t)b(t),$
(42a) $\frac{\partial}{\partial
t}\sigma_{1}^{gf}(t)=-\mathrm{i}\left[\omega_{1}^{ge}+\omega_{1}^{ef}+\chi_{1a}\right]\sigma_{1}^{gf}(t)-\mathrm{i}g_{1a}\sigma_{1}^{ef}(t)a-\mathrm{i}g_{1b}\sigma_{1}^{ge}(t)b(t)+\mathrm{i}\varUpsilon\mathcal{Z}_{1}^{fg}(t)a(t)b(t),$
(42b) $\frac{\partial}{\partial
t}\sigma_{2}^{ge}(t)=-\mathrm{i}\left[\omega_{2}^{ge}-\chi_{2b}{\mathcal{N}}_{b}(t)\right]\sigma_{2}^{ge}(t)+\mathrm{i}g_{2b}\mathcal{Z}_{2}^{eg}(t)b(t).$
(42c)

Using Eq. (13), one can express
$\displaystyle\textstyle\varPhi_{\nu,\nu^{\prime}}(t)$ as
$\displaystyle\textstyle\varPhi_{\nu,\nu^{\prime}}(t)=\langle\varnothing|\alpha_{\nu}\beta_{\nu^{\prime}}|\varPsi(t)\rangle$.
Then, we use the standard quantum-mechanical relations
$\displaystyle\textstyle|\varPsi(t)\rangle=\mathrm{e}^{-\mathrm{i}\mathcal{H}_{\mathrm{eff}}t}|\varPsi(0)\rangle$
and
$\displaystyle\textstyle\mathcal{O}(t)=\mathrm{e}^{\mathrm{i}\mathcal{H}_{\mathrm{eff}}t}\mathcal{O}\mathrm{e}^{-\mathrm{i}\mathcal{H}_{\mathrm{eff}}t}$,
where $\displaystyle\textstyle\mathcal{O}$ and
$\displaystyle\textstyle\mathcal{O}(t)$ denote the operators in the
Schrödinger and Heisenberg representations, correspondingly. This leads us to
the expression
$\displaystyle\textstyle\varPhi_{\nu,\nu^{\prime}}(t)=\langle\varnothing|\alpha_{\nu}(t)\beta_{\nu^{\prime}}(t)|\Psi_{\mathrm{in}}\rangle$,
where we used that
$\displaystyle\textstyle|\varPsi(0)\rangle\equiv|\Psi_{\mathrm{in}}\rangle$
and
$\displaystyle\textstyle\mathrm{e}^{-\mathrm{i}\mathcal{H}_{\mathrm{eff}}t}|\varnothing\rangle=|\varnothing\rangle$.
Using such an approach for the rest of the probability amplitudes entering the
state vector given by Eq. (13), one obtains:

$\boldsymbol{\mu}(t)=\left[\begin{array}[]{c}\varPhi_{\nu,\nu^{\prime}}(t)\\\
\varXi^{\alpha}_{\nu}(t)\\\ \varXi^{\beta}_{\nu}(t)\\\
\varTheta^{\alpha}_{\nu}(t)\\\ \varTheta^{\beta}_{\nu}(t)\\\ R(t)\\\ S(t)\\\
W^{a}(t)\\\ W^{b}(t)\\\
X(t)\end{array}\right]=\big{\langle}\varnothing\big{|}\\!\left[\\!\begin{array}[]{c}\alpha_{\nu}(t)\beta_{\nu^{\prime}}(t)\\\
{\alpha}_{\nu}(t){b}(t)\\\ {\beta}_{\nu}(t){a}(t)\\\
{\alpha}_{\nu}(t){\sigma}_{2}^{ge}(t)\\\
{\beta}_{\nu}(t){\sigma}_{2}^{ge}(t)\\\ {a}(t){b}(t)\\\
{\sigma}_{1}^{gf}(t)\\\ {a}(t){\sigma}_{2}^{ge}(t)\\\
{b}(t){\sigma}_{1}^{ge}(t)\\\
{\sigma}_{1}^{ge}(t){\sigma}_{2}^{ge}(t)\end{array}\\!\right]\\!\big{|}\Psi_{\mathrm{in}}\big{\rangle}.$
(43)

Then, using the equations of motion for the corresponding operators derived
above, we derive the set of evolution equations for the matrix elements
standing on the right-hand side of Eq. (43). This leads us to the set of
evolution equations for the probability amplitudes given by Eq. (15), where we
used the narrowband approximation assuming that $\displaystyle\textstyle
f_{a}(\nu)\approx f_{a}(\omega_{a})$ and $\displaystyle\textstyle
f_{b}(\nu)\approx f_{b}(\omega_{b})$.

## References

  * Kurpiers  _et al._ [2019] P. Kurpiers, M. Pechal, B. Royer, P. Magnard, T. Walter, J. Heinsoo, Y. Salathé, A. Akin, S. Storz, J.-C. Besse, S. Gasparinetti, A. Blais, and A. Wallraff, Quantum communication with time-bin encoded microwave photons, Phys. Rev. Applied 12, 044067 (2019).
  * Lukens and Lougovski [2017] J. M. Lukens and P. Lougovski, Frequency-encoded photonic qubits for scalable quantum information processing, Optica 4, 8 (2017).
  * Brecht  _et al._ [2015] B. Brecht, D. V. Reddy, C. Silberhorn, and M. G. Raymer, Photon temporal modes: A complete framework for quantum information science, Phys. Rev. X 5, 041017 (2015).
  * Ansari  _et al._ [2018] V. Ansari, J. M. Donohue, B. Brecht, and C. Silberhorn, Tailoring nonlinear processes for quantum optics with pulsed temporal-mode encodings, Optica 5, 534 (2018).
  * Schlawin  _et al._ [2018] F. Schlawin, K. E. Dorfman, and S. Mukamel, Entangled two-photon absorption spectroscopy, Accounts of Chemical Research 51, 2207 (2018).
  * Gea-Banacloche [1989] J. Gea-Banacloche, Two-photon absorption of nonclassical light, Phys. Rev. Lett. 62, 1603 (1989).
  * Javanainen and Gould [1990] J. Javanainen and P. L. Gould, Linear intensity dependence of a two-photon transition rate, Phys. Rev. A 41, 5088 (1990).
  * Tabakaev  _et al._ [2021] D. Tabakaev, M. Montagnese, G. Haack, L. Bonacina, J.-P. Wolf, H. Zbinden, and R. T. Thew, Energy-time-entangled two-photon molecular absorption, Phys. Rev. A 103, 033701 (2021).
  * Oka [2018] H. Oka, Two-photon absorption by spectrally shaped entangled photons, Phys. Rev. A 97, 033814 (2018).
  * Oka [2020] H. Oka, Entangled two-photon absorption spectroscopy for optically forbidden transition detection, The Journal of Chemical Physics 152, 044106 (2020).
  * Zhang  _et al._ [2022] Z. Zhang, T. Peng, X. Nie, G. S. Agarwal, and M. O. Scully, Entangled photons enabled time-frequency-resolved coherent Raman spectroscopy and applications to electronic coherences at femtosecond scale, Light: Science & Applications 11, 274 (2022).
  * Lloyd [2008] S. Lloyd, Enhanced sensitivity of photodetection via quantum illumination, Science 321, 1463 (2008).
  * Barzanjeh  _et al._ [2015] S. Barzanjeh, S. Guha, C. Weedbrook, D. Vitali, J. H. Shapiro, and S. Pirandola, Microwave quantum illumination, Phys. Rev. Lett. 114, 080503 (2015).
  * Barzanjeh  _et al._ [2020] S. Barzanjeh, S. Pirandola, D. Vitali, and J. M. Fink, Microwave quantum illumination using a digital receiver, Science Advances 6, eabb0451 (2020).
  * Tittel  _et al._ [2000] W. Tittel, J. Brendel, H. Zbinden, and N. Gisin, Quantum cryptography using entangled photons in energy-time Bell states, Phys. Rev. Lett. 84, 4737 (2000).
  * Ali-Khan  _et al._ [2007] I. Ali-Khan, C. J. Broadbent, and J. C. Howell, Large-alphabet quantum key distribution using energy-time entangled bipartite states, Phys. Rev. Lett. 98, 060503 (2007).
  * Nunn  _et al._ [2013] J. Nunn, L. J. Wright, C. Söller, L. Zhang, I. A. Walmsley, and B. J. Smith, Large-alphabet time-frequency entangled quantum key distribution by means of time-to-frequency conversion, Opt. Express 21, 15959 (2013).
  * Franson [1989] J. D. Franson, Bell inequality for position and time, Phys. Rev. Lett. 62, 2205 (1989).
  * Tittel  _et al._ [1998] W. Tittel, J. Brendel, H. Zbinden, and N. Gisin, Violation of Bell inequalities by photons more than 10 km apart, Phys. Rev. Lett. 81, 3563 (1998).
  * Cabello  _et al._ [2009] A. Cabello, A. Rossi, G. Vallone, F. De Martini, and P. Mataloni, Proposed Bell experiment with genuine energy-time entanglement, Phys. Rev. Lett. 102, 040401 (2009).
  * D’Angelo  _et al._ [2001] M. D’Angelo, M. V. Chekhova, and Y. Shih, Two-photon diffraction and quantum lithography, Phys. Rev. Lett. 87, 013602 (2001).
  * Giovannetti  _et al._ [2001] V. Giovannetti, S. Lloyd, and L. Maccone, Quantum-enhanced positioning and clock synchronization, Nature 412, 417 (2001).
  * Varnavski  _et al._ [2022] O. Varnavski, C. Gunthardt, A. Rehman, G. D. Luker, and T. Goodson, Quantum light-enhanced two-photon imaging of breast cancer cells, The Journal of Physical Chemistry Letters 13, 2772 (2022).
  * Gu  _et al._ [2017] X. Gu, A. F. Kockum, A. Miranowicz, Y. xi Liu, and F. Nori, Microwave photonics with superconducting quantum circuits, Physics Reports 718-719, 1 (2017).
  * Wendin [2017] G. Wendin, Quantum information processing with superconducting circuits: a review, Reports on Progress in Physics 80, 106001 (2017).
  * Blais  _et al._ [2020] A. Blais, S. M. Girvin, and W. D. Oliver, Quantum information processing and quantum optics with circuit quantum electrodynamics, Nature Physics 16, 247 (2020).
  * Marquardt [2007] F. Marquardt, Efficient on-chip source of microwave photon pairs in superconducting circuit QED, Phys. Rev. B 76, 205416 (2007).
  * Chang  _et al._ [2016] Y. Chang, A. González-Tudela, C. Sánchez Muñoz, C. Navarrete-Benlloch, and T. Shi, Deterministic down-converter and continuous photon-pair source within the bad-cavity limit, Phys. Rev. Lett. 117, 203602 (2016).
  * Sánchez-Burillo  _et al._ [2016] E. Sánchez-Burillo, L. Martín-Moreno, J. J. García-Ripoll, and D. Zueco, Full two-photon down-conversion of a single photon, Phys. Rev. A 94, 053814 (2016).
  * Sathyamoorthy  _et al._ [2016] S. R. Sathyamoorthy, A. Bengtsson, S. Bens, M. Simoen, P. Delsing, and G. Johansson, Simple, robust, and on-demand generation of single and correlated photons, Phys. Rev. A 93, 063823 (2016).
  * Leppäkangas  _et al._ [2013] J. Leppäkangas, G. Johansson, M. Marthaler, and M. Fogelström, Nonclassical photon pair production in a voltage-biased Josephson junction, Phys. Rev. Lett. 110, 267004 (2013).
  * Westig  _et al._ [2017] M. Westig, B. Kubala, O. Parlavecchio, Y. Mukharsky, C. Altimiras, P. Joyez, D. Vion, P. Roche, D. Esteve, M. Hofheinz, M. Trif, P. Simon, J. Ankerhold, and F. Portier, Emission of nonclassical radiation by inelastic Cooper pair tunneling, Phys. Rev. Lett. 119, 137001 (2017).
  * Peugeot  _et al._ [2021] A. Peugeot, G. Ménard, S. Dambach, M. Westig, B. Kubala, Y. Mukharsky, C. Altimiras, P. Joyez, D. Vion, P. Roche, D. Esteve, P. Milman, J. Leppäkangas, G. Johansson, M. Hofheinz, J. Ankerhold, and F. Portier, Generating two continuous entangled microwave beams using a dc-biased Josephson junction, Phys. Rev. X 11, 031008 (2021).
  * Flurin  _et al._ [2012] E. Flurin, N. Roch, F. Mallet, M. H. Devoret, and B. Huard, Generating entangled microwave radiation over two transmission lines, Phys. Rev. Lett. 109, 183901 (2012).
  * Gasparinetti  _et al._ [2017] S. Gasparinetti, M. Pechal, J.-C. Besse, M. Mondal, C. Eichler, and A. Wallraff, Correlations and entanglement of microwave photons emitted in a cascade decay, Phys. Rev. Lett. 119, 140504 (2017).
  * Ren  _et al._ [2020] W. Ren, W. Liu, C. Song, H. Li, Q. Guo, Z. Wang, D. Zheng, G. S. Agarwal, M. O. Scully, S.-Y. Zhu, H. Wang, and D.-W. Wang, Simultaneous excitation of two noninteracting atoms with time-frequency correlated photon pairs in a superconducting circuit, Phys. Rev. Lett. 125, 133601 (2020).
  * Perelshtein  _et al._ [2022] M. R. Perelshtein, K. V. Petrovnin, V. Vesterinen, S. Hamedani Raja, I. Lilja, M. Will, A. Savin, S. Simbierowicz, R. Jabdaraghi, J. Lehtinen, L. Grönberg, J. Hassel, M. Prunnila, J. Govenius, G. Paraoanu, and P. Hakonen, Broadband continuous-variable entanglement generation using a Kerr-free Josephson metamaterial, Phys. Rev. Applied 18, 024063 (2022).
  * Kubala  _et al._ [2015] B. Kubala, V. Gramich, and J. Ankerhold, Non-classical light from superconducting resonators coupled to voltage-biased Josephson junctions, Physica Scripta T165, 014029 (2015).
  * Rolland  _et al._ [2019] C. Rolland, A. Peugeot, S. Dambach, M. Westig, B. Kubala, Y. Mukharsky, C. Altimiras, H. le Sueur, P. Joyez, D. Vion, P. Roche, D. Esteve, J. Ankerhold, and F. Portier, Antibunched photons emitted by a dc-biased Josephson junction, Phys. Rev. Lett. 122, 186804 (2019).
  * Grimm  _et al._ [2019] A. Grimm, F. Blanchet, R. Albert, J. Leppäkangas, S. Jebari, D. Hazra, F. Gustavo, J.-L. Thomassin, E. Dupont-Ferrier, F. Portier, and M. Hofheinz, Bright on-demand source of antibunched microwave photons based on inelastic cooper pair tunneling, Phys. Rev. X 9, 021016 (2019).
  * Ma  _et al._ [2021] S.-l. Ma, X.-k. Li, Y.-l. Ren, J.-k. Xie, and F.-l. Li, Antibunched $\displaystyle n$-photon bundles emitted by a Josephson photonic device, Phys. Rev. Research 3, 043020 (2021).
  * li Ma  _et al._ [2022] S. li Ma, J. kun Xie, Y. long Ren, X. ke Li, and F. li Li, Photon-pair blockade in a Josephson-photonics circuit with two nondegenerate microwave resonators, New Journal of Physics 24, 053001 (2022).
  * Roch  _et al._ [2012] N. Roch, E. Flurin, F. Nguyen, P. Morfin, P. Campagne-Ibarcq, M. H. Devoret, and B. Huard, Widely tunable, nondegenerate three-wave mixing microwave device operating near the quantum limit, Phys. Rev. Lett. 108, 147701 (2012).
  * Bergeal  _et al._ [2012] N. Bergeal, F. Schackert, L. Frunzio, and M. H. Devoret, Two-mode correlation of microwave quantum noise generated by parametric down-conversion, Phys. Rev. Lett. 108, 123902 (2012).
  * Clarke and Wilhelm [2008] J. Clarke and F. K. Wilhelm, Superconducting quantum bits, Nature 453, 1031 (2008).
  * Kjaergaard  _et al._ [2020] M. Kjaergaard, M. E. Schwartz, J. Braumüller, P. Krantz, J. I.-J. Wang, S. Gustavsson, and W. D. Oliver, Superconducting qubits: Current state of play, Annual Review of Condensed Matter Physics 11, 369 (2020).
  * Kyriienko and Sørensen [2016] O. Kyriienko and A. S. Sørensen, Continuous-wave single-photon transistor based on a superconducting circuit, Phys. Rev. Lett. 117, 140503 (2016).
  * Stolyarov [2020] E. V. Stolyarov, Single-photon switch controlled by a qubit embedded in an engineered electromagnetic environment, Phys. Rev. A 102, 063709 (2020).
  * Barends  _et al._ [2013] R. Barends, J. Kelly, A. Megrant, D. Sank, E. Jeffrey, Y. Chen, Y. Yin, B. Chiaro, J. Mutus, C. Neill, P. O’Malley, P. Roushan, J. Wenner, T. C. White, A. N. Cleland, and J. M. Martinis, Coherent Josephson qubit suitable for scalable quantum integrated circuits, Phys. Rev. Lett. 111, 080502 (2013).
  * Koch  _et al._ [2007] J. Koch, T. M. Yu, J. Gambetta, A. A. Houck, D. I. Schuster, J. Majer, A. Blais, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Charge-insensitive qubit design derived from the Cooper pair box, Phys. Rev. A 76, 042319 (2007).
  * Houck  _et al._ [2009] A. A. Houck, J. Koch, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Life after charge noise: recent results with transmon qubits, Quantum Inf. Proc. 8, 105 (2009).
  * Chen  _et al._ [2014] Y. Chen, C. Neill, P. Roushan, N. Leung, M. Fang, R. Barends, J. Kelly, B. Campbell, Z. Chen, B. Chiaro, A. Dunsworth, E. Jeffrey, A. Megrant, J. Y. Mutus, P. J. J. O’Malley, C. M. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. C. White, M. R. Geller, A. N. Cleland, and J. M. Martinis, Qubit architecture with high coherence and fast tunable coupling, Phys. Rev. Lett. 113, 220502 (2014).
  * Barends  _et al._ [2014] R. Barends, J. Kelly, A. Megrant, A. Veitia, D. Sank, E. Jeffrey, T. C. White, J. Mutus, A. G. Fowler, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, C. Neill, P. O’Malley, P. Roushan, A. Vainsencher, J. Wenner, A. N. Korotkov, A. N. Cleland, and J. M. Martinis, Superconducting quantum circuits at the surface code threshold for fault tolerance, Nature 508, 500 (2014).
  * Kelly  _et al._ [2015] J. Kelly, R. Barends, A. G. Fowler, A. Megrant, E. Jeffrey, T. C. White, D. Sank, J. Y. Mutus, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, I.-C. Hoi, C. Neill, P. J. J. O’Malley, C. Quintana, P. Roushan, A. Vainsencher, J. Wenner, A. N. Cleland, and J. M. Martinis, State preservation by repetitive error detection in a superconducting quantum circuit, Nature 519, 66 (2015).
  * Ristè  _et al._ [2017] D. Ristè, M. P. da Silva, C. A. Ryan, A. W. Cross, A. D. Córcoles, J. A. Smolin, J. M. Gambetta, J. M. Chow, and B. R. Johnson, Demonstration of quantum advantage in machine learning, npj Quantum Information 3, 16 (2017).
  * Place  _et al._ [2021] A. P. M. Place, L. V. H. Rodgers, P. Mundada, B. M. Smitham, M. Fitzpatrick, Z. Leng, A. Premkumar, J. Bryon, A. Vrajitoarea, S. Sussman, G. Cheng, T. Madhavan, H. K. Babla, X. H. Le, Y. Gang, B. Jäck, A. Gyenis, N. Yao, R. J. Cava, N. P. de Leon, and A. A. Houck, New material platform for superconducting transmon qubits with coherence times exceeding 0.3 milliseconds, Nature Communications 12, 1779 (2021).
  * Wang  _et al._ [2022] C. Wang, X. Li, H. Xu, Z. Li, J. Wang, Z. Yang, Z. Mi, X. Liang, T. Su, C. Yang, G. Wang, W. Wang, Y. Li, M. Chen, C. Li, K. Linghu, J. Han, Y. Zhang, Y. Feng, Y. Song, T. Ma, J. Zhang, R. Wang, P. Zhao, W. Liu, G. Xue, Y. Jin, and H. Yu, Towards practical quantum computers: transmon qubit with a lifetime approaching 0.5 milliseconds, npj Quantum Information 8, 3 (2022).
  * Krantz  _et al._ [2019] P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando, S. Gustavsson, and W. D. Oliver, A quantum engineer’s guide to superconducting qubits, Applied Physics Reviews 6, 021318 (2019).
  * Motzoi  _et al._ [2009] F. Motzoi, J. M. Gambetta, P. Rebentrost, and F. K. Wilhelm, Simple pulses for elimination of leakage in weakly nonlinear qubits, Phys. Rev. Lett. 103, 110501 (2009).
  * Chow  _et al._ [2010] J. M. Chow, L. DiCarlo, J. M. Gambetta, F. Motzoi, L. Frunzio, S. M. Girvin, and R. J. Schoelkopf, Optimized driving of superconducting artificial atoms for improved single-qubit gates, Phys. Rev. A 82, 040305(R) (2010).
  * Machnes  _et al._ [2018] S. Machnes, E. Assémat, D. Tannor, and F. K. Wilhelm, Tunable, flexible, and efficient optimization of control pulses for practical qubits, Phys. Rev. Lett. 120, 150401 (2018).
  * Leonard  _et al._ [2019] E. Leonard, M. A. Beck, J. Nelson, B. G. Christensen, T. Thorbeck, C. Howington, A. Opremcak, I. V. Pechenezhskiy, K. Dodge, N. P. Dupuis, M. D. Hutchings, J. Ku, F. Schlenker, J. Suttle, C. Wilen, S. Zhu, M. G. Vavilov, B. L. T. Plourde, and R. McDermott, Digital coherent control of a superconducting qubit, Phys. Rev. Applied 11, 014009 (2019).
  * Reed  _et al._ [2010] M. D. Reed, B. R. Johnson, A. A. Houck, L. DiCarlo, J. M. Chow, D. I. Schuster, L. Frunzio, and R. J. Schoelkopf, Fast reset and suppressing spontaneous emission of a superconducting qubit, Applied Physics Letters 96, 203110 (2010).
  * Peropadre  _et al._ [2013] B. Peropadre, D. Zueco, F. Wulschner, F. Deppe, A. Marx, R. Gross, and J. J. García-Ripoll, Tunable coupling engineering between superconducting resonators: From sidebands to effective gauge fields, Phys. Rev. B 87, 134504 (2013).
  * Kounalakis  _et al._ [2018] M. Kounalakis, C. Dickel, A. Bruno, N. K. Langford, and G. A. Steele, Tuneable hopping and nonlinear cross-Kerr interactions in a high-coherence superconducting circuit, npj Quantum Information 4, 38 (2018).
  * Wu  _et al._ [2018] Y. Wu, L.-P. Yang, M. Gong, Y. Zheng, H. Deng, Z. Yan, Y. Zhao, K. Huang, A. D. Castellano, W. J. Munro, K. Nemoto, D.-N. Zheng, C. P. Sun, Y.-x. Liu, X. Zhu, and L. Lu, An efficient and compact switch for quantum circuits, npj Quantum Information 4, 50 (2018).
  * Klimov and Sanchez-Soto [2000] A. B. Klimov and L. L. Sanchez-Soto, Method of small rotations and effective Hamiltonians in nonlinear quantum optics, Phys. Rev. A 61, 063802 (2000).
  * Blais  _et al._ [2004] A. Blais, R.-S. Huang, A. Wallraff, S. M. Girvin, and R. J. Schoelkopf, Cavity quantum electrodynamics for superconducting electrical circuits: An architecture for quantum computation, Phys. Rev. A 69, 062320 (2004).
  * Megrant  _et al._ [2012] A. Megrant, C. Neill, R. Barends, B. Chiaro, Y. Chen, L. Feigl, J. Kelly, E. Lucero, M. Mariantoni, P. J. J. O’Malley, D. Sank, A. Vainsencher, J. Wenner, T. C. White, Y. Yin, J. Zhao, C. J. Palmstrøm, J. M. Martinis, and A. N. Cleland, Planar superconducting resonators with internal quality factors above one million, Applied Physics Letters 100, 113510 (2012).
  * Bruno  _et al._ [2015] A. Bruno, G. de Lange, S. Asaad, K. L. van der Enden, N. K. Langford, and L. DiCarlo, Reducing intrinsic loss in superconducting resonators by surface treatment and deep etching of silicon substrates, Applied Physics Letters 106, 182601 (2015).
  * Reagor  _et al._ [2016] M. Reagor, W. Pfaff, C. Axline, R. W. Heeres, N. Ofek, K. Sliwa, E. Holland, C. Wang, J. Blumoff, K. Chou, M. J. Hatridge, L. Frunzio, M. H. Devoret, L. Jiang, and R. J. Schoelkopf, Quantum memory with millisecond coherence in circuit QED, Phys. Rev. B 94, 014506 (2016).
  * Heidler  _et al._ [2021] P. Heidler, C. M. F. Schneider, K. Kustura, C. Gonzalez-Ballestero, O. Romero-Isart, and G. Kirchmair, Non-Markovian effects of two-level systems in a niobium coaxial resonator with a single-photon lifetime of 10 milliseconds, Phys. Rev. Applied 16, 034024 (2021).
  * Paik  _et al._ [2011] H. Paik, D. I. Schuster, L. S. Bishop, G. Kirchmair, G. Catelani, A. P. Sears, B. R. Johnson, M. J. Reagor, L. Frunzio, L. I. Glazman, S. M. Girvin, M. H. Devoret, and R. J. Schoelkopf, Observation of high coherence in Josephson junction qubits measured in a three-dimensional circuit QED architecture, Phys. Rev. Lett. 107, 240501 (2011).
  * Reagor  _et al._ [2013] M. Reagor, H. Paik, G. Catelani, L. Sun, C. Axline, E. Holland, I. M. Pop, N. A. Masluk, T. Brecht, L. Frunzio, M. H. Devoret, L. Glazman, and R. J. Schoelkopf, Reaching 10 ms single photon lifetimes for superconducting aluminum cavities, Applied Physics Letters 102, 192604 (2013).
  * Flurin  _et al._ [2015] E. Flurin, N. Roch, J. D. Pillet, F. Mallet, and B. Huard, Superconducting quantum node for entanglement and storage of microwave radiation, Phys. Rev. Lett. 114, 090503 (2015).
  * Stolyarov [2019] E. V. Stolyarov, Few-photon Fock-state wave packet interacting with a cavity-atom system in a waveguide: Exact quantum state dynamics, Phys. Rev. A 99, 023857 (2019).
  * Sokolov and Stolyarov [2020] A. M. Sokolov and E. V. Stolyarov, Single-photon limit of dispersive readout of a qubit with a photodetector, Phys. Rev. A 101, 042306 (2020).
  * Law  _et al._ [2000] C. K. Law, I. A. Walmsley, and J. H. Eberly, Continuous frequency entanglement: Effective finite Hilbert space and entropy control, Phys. Rev. Lett. 84, 5304 (2000).
  * Eberly [2006] J. H. Eberly, Schmidt analysis of pure-state entanglement, Laser Physics 16, 921 (2006).
  * Bennett  _et al._ [1996] C. H. Bennett, H. J. Bernstein, S. Popescu, and B. Schumacher, Concentrating partial entanglement by local operations, Phys. Rev. A 53, 2046 (1996).
  * Shore and Knight [1993] B. W. Shore and P. L. Knight, The Jaynes-Cummings model, J. Mod. Opt. 40, 1195 (1993).