# Deterministic single-photon source in the ultrastrong coupling regime

Jie Peng<EMAIL_ADDRESS>Hunan Key Laboratory for Micro-Nano Energy
Materials and Devices
and School of Physics and Optoelectronics, Xiangtan University, Hunan 411105,
China  Jianing Tang  Hunan Key Laboratory for Micro-Nano Energy Materials and
Devices
and School of Physics and Optoelectronics, Xiangtan University, Hunan 411105,
China  Pinghua Tang  Hunan Key Laboratory for Micro-Nano Energy Materials and
Devices
and School of Physics and Optoelectronics, Xiangtan University, Hunan 411105,
China  Zhongzhou Ren  School of Physics Science and Engineering, Tongji
University, Shanghai 200092, China  Junlong Tian  College of Big Data and
Information Engineering, Guizhou University, Guiyang 550025, China  Nancy
Barraza  International Center of Quantum Artificial Intelligence for Science
and Technology (QuArtist) and Department of Physics, Shanghai University,
200444 Shanghai, China  Gabriel Alvarado Barrios  Kipu Quantum,
Greifswalderstrasse 226, 10405 Berlin, Germany  Lucas Lamata  Departamento de
Física Atómica, Molecular y Nuclear, Universidad de Sevilla, 41080 Sevilla,
Spain  Instituto Carlos I de Física Teórica y Computacional, 18071 Granada,
Spain  Enrique Solano<EMAIL_ADDRESS>International Center of Quantum
Artificial Intelligence for Science and Technology (QuArtist) and Department
of Physics, Shanghai University, 200444 Shanghai, China  Kipu Quantum,
Greifswalderstrasse 226, 10405 Berlin, Germany  IKERBASQUE, Basque Foundation
for Science, Plaza Euskadi 5, 48009 Bilbao, Spain  F. Albarrán-Arriagada
Departamento de Física, Universidad de Santiago de Chile (USACH), Avenida
Víctor Jara 3493, 9170124, Estación Central, Chile  Center for the Development
of Nanoscience and Nanotechnology, 9170124, Estación Central, Chile

###### Abstract

Deterministic single-photon sources are important and ubiquitous in quantum
information protocols. However, to the best of our knowledge, none of them
work in the ultrastrong light-matter coupling regime, and each excitation
process can only emit one photon. We propose a deterministic single-photon
source in circuit QED which can work in the ultrastrong coupling regime. Here,
two qubits are excited simultaneously in one process and two deterministic
single photons can be sequentially emitted with an arbitrary time separation.
This happens through two consecutive adiabatic transfers along the one-photon
solutions of the two-qubit Rabi and Jaynes-Cummings model, which has constant
eigenenergy in the whole coupling regime. Unlike the stimulated Raman
adiabatic passage, the system goes back to the initial state of another period
automatically after photon emission. Our scheme can approach unity single-
photon efficiency, indistinguishability, and purity simultaneously. With the
assistance of the Stark shift, a deterministic single photon can be generated
within a time proportional to the inverse of the resonator frequency.

_Introduction.–_ Single photon sources are fundamental building blocks in
quantum information, with applications ranging from quantum computation klm ;
kok  to quantum communication cirac ; pan1 ; pan2  and sensing degen . The
different technologies for single-photon sources can be classified into two
families thomas . The first one considers a nonlinear material process such as
spontaneous parametric downconversion lou ; burn , where a pump laser
illuminates a material with a $\chi(2)$ optical nonlinearity, creating two
photons, one of which can be used to herald the creation of the other.
However, this process is probabilistic and inefficient bon ; thomas . The
second approach is based on single quantum emitters, which deterministically
emit one photon at a time. This has been demonstrated in atoms Darquie ,
molecules Christian , ions Maunz , color centers in diamonds Kurtsiefer , and
quantum dots Charles . However, the spontaneous emission in all directions
makes photon collection difficult. One way to solve this problem is to couple
the emitters to a cavity to enhance the radiation into the cavity mode, so
that the photon can be emitted from the cavity through a certain direction ion
. Meanwhile, the emission rate can be increased by Purcell effect circuit ,
and the pure dephasing is reduced due to the coupling to a fixed cavity mode,
making the photons more indistinguishable tunable . Using a quantum dot
coupled to cavities, near-unity indistinguishability and purity are realized
simultaneously, with an extraction efficiency of $66\%$ dx , a polarized
single-photon efficiency of $60\%$ wang  and an overall efficiency of $57\%$
tomm , respectively. Besides high efficiency tunable , purity Bozyigit  and
indistinguishability lang , a microwave single photon source realized in
circuit QED can also be tunable tunable ; peng .

The Purcell effect, which enhances the emission rate, is proportional to the
qubit-cavity coupling strength Kaer . When we enter the strong coupling
regime, the swap between the qubit excited state and the single photon state
in the cavity will be accelerated. Since now the ultrastrong ultra  and even
deep-strong coupling deep  have been realized in circuit QED, it is natural to
consider implementing a single photon source using these stronger coupling
strengths, which most likely will accelerate the photon generation speed.
However, to the best of our knowledge, this has not been realized yet. Because
the ultrastrong coupling will invoke counter-rotating terms, which excite the
qubit and photon simultaneously, all photon number states become connected,
and the system, described by the quantum Rabi model braak ; chen , normally
involve an infinite number of photons zhong ; xie ; yu .

Recently, we have found special dark states of the two-qubit quantum Rabi
model pj1 ; pj2 ; pj3 , which contain at most one photon, and have constant
eigenenergy in the whole coupling regime. Here, we propose to implement a
deterministic single photon source that takes advantage of the ultrastrong
coupling using these dark states. Our scheme begins with qubits and the cavity
ground state, and then the qubits are excited simultaneously. Next, the system
undergoes an adiabatic evolution along the aforementioned dark state, and end
up with a product state of a single photon state and a two-qubit singlet
state. Because the reaching of ultrastrong coupling and peculiarity of the
dark state, the target state can be generated with fidelity $99.2\%$ in a time
of $68\omega^{-1}$, where $\omega$ is the cavity frequency. Furthermore, with
the addition of Stark shift terms Grimsmo1 ; Grimsmo2 ; cong , the minimum
energy gap between the dark state and its closest eigenstates can be increased
to $\approx 0.54\omega$, which accelerates the adiabatic evolution to
$12\omega^{-1}$ with fidelity reaching $99\%$, showing a sign of ultrafast
state generation ultrafast . Moreover, after the photon is emitted, the system
will be in an eigenstate of the two-qubit JC model with the excitation number
$C=1$, which has also constant eigenenergy in the whole coupling regime. So
through another adiabatic process, a single photon can be generated and
emitted through dissipation, leaving the system in its ground state, which is
also the initial state of a period. Unlike the STIRAP, no recycling of the
qubit between photon generations is required, which will increase the
repetition rate Wilk . Here one excitation process emits two single photons
with efficiencies over 99% with an arbitrary time separation, and their purity
and indistiguishability both approach unity. Finally, we propose an
experimental realization in circuit QED.

_Scheme and circuit QED implementation.–_ Our scheme of the deterministic
single-photon source is based on the one-photon solutions to the two-qubit
Rabi and JC models ($\hbar=1$)

$\displaystyle H_{Rabi}$ $\displaystyle=$ $\displaystyle\omega
a^{{\dagger}}a+g_{1}{\sigma_{1x}}(a+{a^{\dagger}})+g_{2}{\sigma_{2x}}(a+a^{{\dagger}})$
(1) $\displaystyle+\Delta_{1}\sigma_{1z}+\Delta_{2}{\sigma_{2z}},$
$\displaystyle H_{JC}$ $\displaystyle=$ $\displaystyle\omega
a^{\dagger}a+g_{1}(\sigma_{1}^{\dagger}a+\sigma_{1}a^{\dagger})+g_{2}(\sigma_{2}^{{\dagger}}a+\sigma_{2}a^{\dagger})$
(2) $\displaystyle+\Delta_{1}\sigma_{1z}+\Delta_{2}\sigma_{2z},$

where $a^{\dagger}$ and $a$ are the photon creation and annihilation operators
with frequency $\omega$, respectively. Also, $\sigma_{j\alpha}(\alpha=x,y,z)$
are the Pauli matrices corresponding to the $j$-th qubit. $2\Delta_{j}$ is the
energy level splitting of the $j$-th qubit, and $g_{j}$ is the qubit-photon
coupling parameter between the resonator and $j$-th qubit. Normally, there is
no solution with finite photon numbers to $H_{Rabi}$ since all the photon
number states are connected. However, we have found special dark states with
at most one-photon pj1

$\displaystyle|\psi_{R}\rangle=\frac{1}{{\cal
N}}[(\Delta_{1}-\Delta_{2})|0\uparrow\uparrow\rangle+g|1\rangle(\downarrow\uparrow\rangle-\uparrow\downarrow\rangle)],$
(3)

where $\Delta_{1}+\Delta_{2}=\omega$ and $g_{1}=g_{2}=g$, corresponding to the
horizontal line $E=\omega$ in Fig. 1(a). Surprisingly, horizontal lines
$E/\omega=N~{}(N=0,1,2,\ldots)$ emerge in the energy spectrum of the two-qubit
JC model under the same conditions, as shown in Fig. 1(b).

Figure 1: (a) The numerical spectrum of two-qubit quantum Rabi model with
$\Delta_{1}=0.8\omega$, $\Delta_{2}=0.2\omega$, $g=g_{1}=g_{2}$. The red lines
have even parity while the blue lines have odd parity. (b) The spectrum of
two-qubit JC model with the same parameters. (c) Relevant energy levels and
transitions of our scheme. (d) Setup for the scheme: Two superconducting
qubits are coupled to one resonator, whose photon emission rate into the TL is
controlled by a variable coupler C. The lower part is a superconducting
circuit design for the two-qubit Rabi and JC models with tunable couplings and
qubit frequencies sl .

The corresponding eigenstates read pj4

$\displaystyle|\psi_{JC}\rangle=\frac{1}{{\cal
N^{\prime}}}[(\Delta_{1}-\Delta_{2})|1\downarrow\downarrow\rangle+g(|0\downarrow\uparrow\rangle-|0\uparrow\downarrow\rangle)]$
(4)

for $N=0$ pj4 . We will show that $|\psi_{R}\rangle$ and $|\psi_{JC}\rangle$
can be used to produce a special single-photon source through two consecutive
adiabatic transfers.

Our scheme is depicted in Fig. 1 (c) and (d). Two qubits are coupled to one
resonator, which is connected to a transmission line (TL) through a variable
coupler C, so that its dissipation rate $\kappa_{c}$ is tunable. The qubits
and the resonator are cooled down to the ground state
$|0\downarrow\downarrow\rangle$ initially. Then pumping pulses,

$H=\frac{\Omega}{2}(\sigma_{1}^{\dagger}{e^{-i\omega_{q1}t}}+\sigma_{1}{e^{i\omega_{q1}t}}+\sigma_{2}^{\dagger}{e^{-i\omega_{q2}t}}+\sigma_{2}{e^{i\omega_{q2}t}}),$
(5)

are applied to excite the qubits, where $\omega_{q1}$, $\omega_{q2}$ are the
frequencies of the two qubits, respectively. Initially, we set
$\omega_{q1}=2\Delta_{1}\approx 2\omega$, $\omega_{q2}=2\Delta_{2}\approx 0$,
so the crosstalk can be neglected and the coupling $g=g_{1}=g_{2}$ is set to
zero. After excitation, $|0\uparrow\uparrow\rangle$ just corresponds to
$|\psi_{R}\rangle$ of Eq. (3) since $g=0$ and $\Delta_{1}\neq\Delta_{2}$.
Next, we increase $g$ to a nonzero value and decrease $\Delta_{1}-\Delta_{2}$
to $0$, so that the state $|0\uparrow\uparrow\rangle$ evolves adiabatically to
$|1\downarrow\uparrow\rangle-|1\uparrow\downarrow\rangle$ through
$|\psi_{R}\rangle$. Then, we activate C to increase the dissipation rate into
the TL to a very large value, which can be achieved in a few

Figure 2: (a) Shows the evolution of system states during one period, where
$|\psi_{B}\rangle=|\downarrow\uparrow-\uparrow\downarrow\rangle$. (b) The
corresponding evolution of the parameters. (c) The waveforms of the emitted
photons. The solid line corresponds to the numerical simulation, fitted by an
exponential decaying function $\exp(-\kappa t)$ and a Gaussian function,
represented by the dashed line. (d) A sequence of single photons are generated
when we repeat the above processes.

nanoseconds catch . The single photon can be emitted into the TL and the
system state becomes
$|0\downarrow\uparrow\rangle-|0\uparrow\downarrow\rangle$, corresponding to
the one-photon solution $|\psi_{JC}\rangle$ of Eq. (4) for the two-qubit JC
model. After an arbitrary time separation, we begin another adiabatic
evolution along $|\psi_{JC}\rangle$ by decreasing $g$ to zero and increasing
$\Delta_{1}-\Delta_{2}$ to a nonzero value, which generates
$|1\downarrow\downarrow\rangle$, and the single photon can be emitted during
this process. After that, the system returns to the initial state
$|0\downarrow\downarrow\rangle$ which can be utilized for the next period
automatically. Furthermore, we propose a superconducting circuit design for
the two-qubit Rabi and JC models with tunable couplings and qubit frequencies,
shown in the lower part of Fig. 1 (d), with a detailed demonstration in sl .

Usually, one excitation process can only emit one photon for deterministic
single-photon sources, and the system has to be repumped afterwards, so there
is always an unavoidable time separation between the photons emitted. However,
with our protocol one pumping process can emit two single photons, since the
system reaches the “excited state” of the second emission process
automatically after the first photon is emitted, which also provides an
arbitrarily controllable time separation of the photons.

We simulate the process described above using the Lindblad master equation,
with the resonator intrinsic dissipation rate $\kappa_{in}=10^{-4}\omega$, the
qubit damping rate $\gamma=10^{-5}\omega$ and dephasing rate
$\gamma_{\phi}=2\times 10^{-5}\omega$ sl . The evolution of different states
in a period is shown in Fig. 2 (a), with corresponding parameters depicted in
Fig. 2 (b). The absolute values of the detunnings of both qubits with respect
to the resonator are always the same during the evolution, which implies that
the magnetic flux used to tune the qubit frequencies can also be used to tune
the coupling strength, maintaining $g_{1}=g_{2}$ sl . This will make our
implementation easier to realize. Note that $g$ can be arbitrary values during
$100\omega^{-1}-150\omega^{-1}$, because the qubit singlet state is decoupled
from the cavity when emitting the first photon. Here, we choose $g=0$ to have
a better single photon purity. Although the fast adiabatic evolution along
$|\psi_{R}\rangle$ may seem difficult according to the adiabatic theorem,
which states that the system will evolve along $|\psi_{R}\rangle$ if

$\left|\frac{\langle
E_{n}(t)|\dot{H}_{R}|\psi_{R}(t)\rangle}{(E_{n}-E_{R})^{2}}\right|\ll 1,$ (6)

where there is an energy level $|E_{n}\rangle$ very close to
$|\psi_{R}\rangle$, as shown in Fig. 1 (a). However,
$|1(\downarrow\uparrow-\uparrow\downarrow)\rangle$ can be fast generated from
$|\uparrow\uparrow\rangle$ in $68\omega^{-1}$ with a fidelity of $99.2\%$
through the adiabatic evolution along $|\psi_{R}\rangle$. This is mainly
because the peculiarity of the dark state and reaching of ultrastrong coupling
regime sl . The photon emission rate into the TL is shown in Fig. 2 (c), where
the time separation between two emitted photons is zero. The first photon has
a typical single-sided exponential waveform decaying with the lifetime of the
resonator, while the second photon has a Gaussian shape which is optimal for
tolerance to mode mismatch Sonia . We use an exponential function
$\exp(-\kappa t)$ and a Gaussian function to fit the simulated data and find
good consistency, as can be seen in Fig. 2 (c). The system is reset to
$|0\downarrow\downarrow\rangle$ automatically after the second photon
emission, so we can repeat the above process to obtain a sequence of single
photons, as shown in Fig. 2 (d).

_Figures of merit–_ There are commonly three important figures of merit for a
single photon source Reimer : efficiency, purity and indistinguishability.
First, the generation and collection efficiency are defined as the photon
generated and collected for one excitation pulse, respectively, which is equal
to one for a perfect single photon source. Here, two single photons are
generated in one excitation process, with emission probabilities
(efficiencies), $\kappa_{c}\int a^{{\dagger}}(t)a(t){\rm d}t$, both larger
than $99\%$. The generation and collection efficiencies both approach unity,
since the qubits are strongly coupled to the resonator, and their interaction
is much faster than the decoherence rate of the qubit. The single photon is
generated in the resonator and almost fully collected by the TL through the
variable coupler because $\kappa_{c}$ can be tuned to be $1000$ times the
intrinsic decay rate in a few nanoseconds catch .

Figure 3: (a)The simulated second-order correlation function in HBT experiment
for the single-photon source proposed above, which characterizes its purity.
(b) Our scheme to collect two single photons separately. (c) The second-order
correlation function in HBT experiment for the first single photon source. (d)
The second-order correlation function in HBT experiment for the second single
photon source. (e) The second-order correlation function in HOM experiment for
the first single photon source, which characterizes its indistinguishability.
(f) The second-order correlation function in HOM experiment for the second
single photon source. All these simulated correlation functions are normalized
to the highest peak.

The second requirement for a single photon source is high purity, which means
exactly one photon is emitted at a time. It can be characterized by the
second-order correlation functions in the Hanbury-Brown-Twiss (HBT) experiment

$\displaystyle G^{(2)}(t,\tau)$ $\displaystyle=$
$\displaystyle\langle{a^{\dagger}}(t){a^{\dagger}}(t+\tau)a(t+\tau)a(t)\rangle,$
$\displaystyle G^{(2)}(\tau)$ $\displaystyle=$ $\displaystyle\int
G^{(2)}(t,\tau)dt,$ (7)

where the former represents the probability of detecting a photon at time $t$
and another one at time $t+\tau$, and the later is its integration over $t$,
as shown in Fig. 3 (a). The data has been normalized to the highest peak.
There is a strongly suppressed $G^{(2)}(0)\approx 0.025$, giving a clear
evidence of the single-photon emission. $G^{(2)}(\tau)$ has the same period as
the single-photon source, and two different kinds of peaks other than
$\tau=0$, because the waveforms of photons in a period are not identical. The
first peak arises at $\tau\approx 88\omega^{-1}$, corresponding to the
coincidence detection of the first and the second photons, and the second peak
arises at $\tau\approx 173\omega^{-1}$, which has the same height as the first
peak caused by the coincidence detection of the second and the third photons.
The highest peak arises clearly at $\tau\approx 260\omega^{-1}=T_{period}$
because of the periodicity.

The third requirement is high indistinguishability, which means the photons
must be identical in all degrees of freedom, to ensure the successful
implementation of the two-photon gate through interference. The Hong, Ou and
Mandel (HOM) experiment hom  is used to characterize this notion, where two
photons are interfered on a 50/50 beam splitter. If they are completely
indistinguishable, they will coalesce and exit through the same beam splitter
output. Here two single photons are emitted in one period with different
temporal waveforms which may be made identical with the photon waveform
reshaping technology Kie  to implement a two-photon CNOT gate. Since the time
separation between these two single photons is arbitrary, we have much freedom
to choose the delay time between the control and target photons.

To make better use of this single photon source, we can add another output
channel to the resonator with a variable coupler to collect these two single
photons separately, as shown in Fig. 3 (b) , so that the photons in each
channel are indistinguishable, and with only one excitation process in a
period, we efficiently generate two sequences of indistinguishable single
photons. The first one has a typical single-sided exponential waveform. Its
simulated $G^{(2)}(\tau)$ is shown in Fig. 3 (c), with $G^{(2)}(0)\approx
0.026$. The second one has a Gaussian shape, whose $G^{(2)}(\tau)$ is shown in
Fig. 3 (d), with $G^{(2)}(0)\approx 0.02$. They can be applied in different
situations. We also simulated the normalized HOM experiment counts
$G^{(2)}_{HOM}(\tau)$ of photon detection in different outputs at time delay
$\tau$. There is a clear vanishing second order correlation function
$G^{(2)}_{HOM}(0)$ for the first photon shown in Fig. 3 (e) because of the
indistinguishability, while typical small peaks around $\tau=0$ arise for the
second photon in Fig. 4 (f), because the decoherence of the qubits cause
random photon frequency differences. However, the first photon has already
been generated before emission, and its spatio-temporal coherence is governed
by resonator decay, showing no significant additional dephasing lang . The
indistinguishability is defined by dividing the area of the peaks around
$\tau=0$ by that of the uncorrelated peak around $\tau=T$, and subtracting
this number from unity. It can be calculated as Kaer

$I=\frac{\int_{0}^{\infty}{\rm d}t\int_{0}^{\infty}{\rm d}\tau|\langle
a^{{\dagger}}(t+\tau)a(t)\rangle|^{2}}{\int_{0}^{\infty}{\rm
d}t\int_{0}^{\infty}{\rm d}\tau\langle
a^{\dagger}(t+\tau)a(t+\tau)\rangle\langle a(t)^{{\dagger}}a(t)\rangle}.$ (8)

Because of periodicity, we only need to consider the integration over one
period quantum-dot . The indistinguishability, $I$, reaches $99.99\%$ for the
first photon and $96.1\%$ for the second photon.

Another way to generate indistinguishable photons is to turn off the variable
coupler $C$ during the adiabatic evolution from
$|0\rangle(|\downarrow\uparrow\rangle-|\uparrow\downarrow)$ to
$|1\downarrow\downarrow\rangle$, so the second photon will be generated inside
the resonator and then emitted into the TL when $C$ is turned on, just like
the first photon. These two single photons will clearly be indistinguishable,
with $I\approx 1$.

_Ultrafast generation of the single photon with the assistance of the Stark
shift.–_ Another requirement for the single photon source is high speed,
meanwhile, the most prominent advantage of the ultratrong coupling is the
ultrafast state generation, which can be realized here by adding Stark shift
terms to the two-qubit Rabi model Eckle ; yf ; li1 ; li2

Figure 4: (a) The numerical spectrum of $H_{RS}$ with even parity when
parameters evolve according to (b). (b) The time-dependent parameters during
the adiabatic evolution from $|0\uparrow\uparrow\rangle$ to
$|1\psi_{B}\rangle$ along $|\psi_{RS}\rangle$, where
$U=U_{1}=U_{2}=0.5\omega$. (c) The numerical spectrum of $H_{RS}$ of even
parity with the same parameters as (a) except for $U_{1}=U_{2}=0.01\omega$.
(d) Ultrafast generation of $|1\psi_{B}\rangle$ through adiabatic evolution
along $|\psi_{RS}\rangle$, with parameters evolving according to (b).

$\displaystyle H_{RS}$ $\displaystyle=$ $\displaystyle\omega
a^{{\dagger}}a+g_{1}{\sigma_{1x}}(a+{a^{\dagger}})+g_{2}{\sigma_{2x}}(a+a^{{\dagger}})$
(9)
$\displaystyle+\Delta_{1}\sigma_{1z}+\Delta_{2}{\sigma_{2z}}+U_{1}a^{{\dagger}}a\sigma_{1z}+U_{2}a^{{\dagger}}a\sigma_{2z},$

where $U_{1}$ and $U_{2}$ are couplings of the stark terms. The photon
frequency is shifted by $U_{1}\sigma_{1z}+U_{2}\sigma_{2z}$, so the stability
of the system requires $U_{1}+U_{2}\leq\omega$. We find a similar special dark
state

$\displaystyle|\psi_{RS}\rangle$ $\displaystyle=$ $\displaystyle\frac{1}{{\cal
N}}[(\Delta_{1}-\Delta_{2}+U_{1}-U_{2})|0\uparrow\uparrow\rangle$ (10)
$\displaystyle+g(|1\downarrow\uparrow\rangle-|1\uparrow\downarrow\rangle)]$

with constant energy $E=\omega$ in the whole coupling regime, under the same
condition $\Delta_{1}+\Delta_{2}=\omega$ and $g=g_{1}=g_{2}$ as for $H_{Rabi}$
sl . Interestingly, the energy gap between $|\psi_{RS}\rangle$ and its closest
eigenstates is enlarged by the stark terms. E. g.,
$|\psi_{RS}\rangle=|0\uparrow\uparrow\rangle$ and
$|2\downarrow\downarrow\rangle$ are degenerate for $H_{Rabi}$ at $g=0$, but
the stark terms will reduce the energy of the latter by $U_{1}+U_{2}$, so
these two energy levels are separated in the spectrum. This phenomena is quite
significant and interesting when $U_{1}+U_{2}=\omega$, where all
$|n\downarrow\downarrow\rangle$ become degenerate ground states at $g=0$ with
$E=-\omega$, as shown in Fig. 4 (a). These infinite many energy levels lie
below $E=-\omega$ in the spectrum, which should have appeared in the upper
part if $U_{1}+U_{2}<\omega$, as can be seen from Fig. 4 (c). All these states
become quite far from $|\psi_{RS}\rangle$, and the minimum gap between
$|\psi_{RS}\rangle$ and other states $|E_{n}\rangle$ is about $0.54~{}\omega$,
as can be seen in Fig. 4 (a). Therefore, according to the adiabatic theorem
Eq. (6), the adiabatic process should be much faster than for cases where the
minimum gap is much smaller than $\omega$. With the parameters changing as
Fig. 4 (b), $|1\psi_{B}\rangle$ can be generated from
$|0\uparrow\uparrow\rangle$ through an adiabatic evolution along
$|\psi_{RS}\rangle$ with fidelity $99.2\%$, operating at a time proportional
to the inverse of the resonator frequency $f$
($12\omega^{-1}=\frac{6}{\pi}f^{-1}$), which is a sign of ultrafast quantum-
state generation ultrafast , as shown in Fig. 4 (d). Note that many other
choices of $U$ also significantly accelerate the adiabatic evolution. E. g.,
$|1\psi_{B}\rangle$ can be generated in $12\omega^{-1}$ with fidelity $99\%$
for $U_{1}=U_{2}=0.45\omega$. Following the same procedure detailed above, a
faster single-photon source can be implemented.

_Conclusion.–_ We have proposed a scheme for a deterministic single-photon
source that can work in the ultrastrong coupling regime. That is through two
consecutive adiabatic transfers along the one-photon solutions of the two-
qubit Rabi model and JC model, respectively. An important advantage of our
scheme is that one pumping process can emit two deterministic single photons
with an arbitrary time separation that is easily controlled. Furthermore, with
our protocol, the system evolves naturally to the initial state of the next
period after photon emission. We characterize our single photon source by
calculating efficiency, purity, and indistinguishability, showing that all of
them can reach near unity values simultaneously. Moreover, by introducing
Stark shift terms, we can accelerate the speed of the single-photon generation
to a degree proportional to the inverse of the resonator frequency. Our scheme
paves the way for the application of ultrastrong coupling in fast computation
and deterministic state generation.

_Acknowledgements.–_ We acknowledge Guillermo Romero for helpful discussions.
This work was supported by the Natural Science Foundation of Hunan Province,
China (2022JJ30556), the National Natural Science Foundation of China (Grant
No. 11874316), the National Basic Research Program of China (2015CB921103),
the Program for Changjiang Scholars and Innovative Research Team in University
(No. IRT13093), Junta de Andalucía (P20_00617 and US-1380840), ANID Subvención
a la Instalación en la Academia SA77210018, ANID Proyecto Basal AFB 180001.

## References

  * (1) E. Knill, R. Laflamme, and G. J. Milburn, A scheme for efficient quantum computation with linear optics, Nature 409, 46 (2001).
  * (2) P. Kok, W. J. Munro, K. Nemoto, T. C. Ralph, J. P. Dowling, and G. J. Milburn, Linear optical quantum computing with photonic qubits, Rev. Mod. Phys. 79, 135 (2007).
  * (3) I. Cirac, P. Zoller, J. Kimble, and H. Mabuchi, Quantum State Transfer and Entanglement Distribution among Distant Nodes in a Quantum Network, Phys. Rev. Lett. 78, 3221 (1997).
  * (4) J. W. Pan, Z. B. Chen, C. Y. Lu, H. Weinfurter, A. Zeilinger, and M. Żukowski, Multiphoton entanglement and interferometry. Rev. Mod. Phys. 84, 777 (2012).
  * (5) J. Yin et al., Satellite-based entanglement distribution over 1,200 kilometers.  Science 356, 1140–1144 (2017).
  * (6) C. L. Degen, F. Reinhard, and P. Cappellaro, Quantum sensing, Rev. Mod. Phys. 89, 035002 (2017).
  * (7) S. Thomas and P. Senellart, The race for the ideal single-photon source is on, Nat. Nanotechnol. 16, 367 (2021).
  * (8) W. H. Louisell, A. Yariv, and A. E. Siegman, Quantum Fluctuations and Noise in Parametric Processes. I., Phys. Rev. 124, 1646 (1961).
  * (9) D. C. Burnham and D. L. Weinberg, Observation of Simultaneity in Parametric Production of Optical Photon Pairs, Phys. Rev. Lett. 25, 84 (1970).
  * (10) D. Bonneau, G. J. Mendoza, J. L. O’Brien and M. G. Thompso, Effect of loss on multiplexed single-photon sources, New J. Phys. 17, 043057 (2015).
  * (11) B. Darquie, M. P. A. Jones, J. Dingjan, J. Beugnon, S. Bergamini, Y. Sortais, G. Messin, A. Browaeys and P. Grangier, Controlled Single-Photon Emission from a Single Trapped Two-Level Atom, Science 309, 454 (2005).
  * (12) C. Brunel, B. Lounis, P. Tamarat, and M. Orrit, Triggered Source of Single Photons based on Controlled Single Molecule Fluorescence, Phys. Rev. Lett. 83, 2722 (1999).
  * (13) P. Maunz, D. L. Moehring, S. Olmschenk, K. C. Younge, D. N. Matsukevich, and C. Monroe, Nat. Phys. 3, 538 (2007).
  * (14) C. Kurtsiefer, S. Mayer, P. Zarda, and H. Weinfurter, Stable Solid-State Source of Single Photons, Phys. Rev. Lett. 85, 290 (2000).
  * (15) C. Santori, M. Pelton, G. Solomon, Y. Dale, and Y. Yamamoto, Triggered Single Photons from a Quantum Dot, Phys. Rev. Lett. 86, 1502 (2001).
  * (16) H. G. Barros, A. Stute, T. E. Northup, C. Russo1, P. O. Schmidt and R. Blatt, Deterministic single-photon source from a single ion, New J. Phys. 11, 103004 (2009).
  * (17) A. A. Houck, D. I. Schuster, J. M. Gambetta, J. A. Schreier, B. R. Johnson, J. M. Chow, L. Frunzio, J. Majer, M. H. Devoret, S. M. Girvin, and R. J. Schoelkopf, Generating single microwave photons in a circuit, Nature 449, 328 (2007).
  * (18) Y. Zhou, Z. H. Peng, Y. Horiuchi, O. V. Astafiev, and J.S. Tsai, Tunable Microwave Single-Photon Source Based on Transmon Qubit with High Efficiency, Phys. Rev. Applied 13, 034007 (2020).
  * (19) X. Ding, Y. He, Z. C. Duan, N. Gregersen, M. C. Chen, S. Unsleber, S. Maier, C. Schneider, M. Kamp, S. Höfling, C. Y. Lu, and J. W. Pan, On-Demand Single Photons with High Extraction Efficiency and Near-Unity Indistinguishability from a Resonantly Driven Quantum Dot in a Micropillar, Phys. Rev. Lett. 116, 020401 (2016).
  * (20) H. Wang, Y. M. He, T. H. Chung et al., Towards optimal single-photon sources from polarized microcavities, Nat. Photonics 13, 770 (2019).
  * (21) N. Tomm, A. Javadi, N. O. Antoniadis et al., A bright and fast source of coherent single photons, Nat. Nanotechnol. 16, 399 (2021).
  * (22) D. Bozyigit, C. Lang, L. Steffen, J. M. Fink, C. Eichler, M. Baur, R. Bianchetti, P. J. Leek, S. Filipp, M. P. da Silva , A. Blais, and A. Wallraff, Antibunching of microwave-frequency photons observed in correlation measurements using linear detectors, Nature Phys. 7, 154 (2011).
  * (23) C. Lang, C. Eichler, L. Steffen, J. M. Fink, M. J. Woolley, A. Blais, and A. Wallraff, Correlations, indistinguishability and entanglement in Hong–Ou–Mandel experiments at microwave frequencies, Nature Phys. 9, 345 (2013).
  * (24) Z. H. Peng, S. E. de Graaf, J. S. Tsai, and O. V. Astafiev, Tuneable on-demand single-photon source in the microwave range, Nat. Commun. 7, 12588 (2016).
  * (25) P. Kaer, N. Gregersen, and J. Mork, The role of phonon scattering in the indistinguishability of photons emitted from semiconductor cavity QED systems, New J. Phys. 15, 035027 (2013).
  * (26) T. Niemczyk, F. Deppe, H. Hueb, E. P. Menzel, F. Hocke, M. J. Schwarz, J. J. Garcia-Ripoll, D. Zueco, T. Hümmer, E. Solano, A. Marx, and R. Gross, Circuit quantum electrodynamics in the ultrastrong-coupling regime,  Nature Phys. 6, 772 (2010).
  * (27) F. Yoshihara, T. Fuse, S. Ashhab, K. Kakuyanagi, S. Saito, and K. Semba, Superconducting qubit–oscillator circuit beyond the ultrastrong-coupling regime, Nature Phys. 13, 44 (2017).
  * (28) D. Braak, Integrability of the Rabi model, Phys. Rev. Lett. 107, 100401 (2011).
  * (29) Q. H. Chen, C. Wang, S. He, T. Liu, and K. L. Wang, Exact solvability of the quantum Rabi model using Bogoliubov operators, Phys. Rev. A 86, 023822 (2012).
  * (30) H. H. Zhong, Q. T. Xie, M. T. Batchelor, and C. H. Lee, Analytical eigenstates for the quantum Rabi model, J. Phys. A: Math. Theor. 46, 415302 (2013).
  * (31) Q. T. Xie, H. H. Zhong, M. T. Batchelor, and C. H. Lee, The quantum Rabi model: solution and dynamics,  J. Phys. A: Math. Theor. 50, 113001 (2017).
  * (32) Y. Y. Zhang, Q. H. Chen, and Y. Zhao, Generalized rotating-wave approximation to biased qubit-oscillator systems, Phys. Rev. A 87, 033827 (2013).
  * (33) J. Peng, Z. Z. Ren, D. Braak, G. J. Guo, G. X. Ju, X. Zhang, and X. Y. Guo, Solution of the two-qubit quantum Rabi model and its exceptional eigenstates,  J. Phys. A: Math. Theor. 47, 265303 (2014).
  * (34) J. Peng, C. X. Zheng, G. J. Guo, X. Y. Guo, X. Zhang, C. S. Deng, G. X. Ju, Z. Z. Ren, L. Lamata, and E. Solano, Dark-like states for the multi-qubit and multi-photon Rabi models, J. Phys. A: Math. Theor. 50, 174003 (2017).
  * (35) J. Peng, J. C. Zheng, J. Yu, P. H. Tang, G. A. Barrios, J. X. Zhong, E. Solano, F. Albarrán-Arriagada, and L. Lamata, One-Photon Solutions to the Multiqubit Multimode Quantum Rabi Model for Fast $W$-State Generation, Phys. Rev. Lett. 127, 043604 (2021).
  * (36) A. L. Grimsmo and S. Parkins, Cavity-QED simulation of qubit-oscillator dynamics in the ultrastrong-coupling regime, Phys. Rev. A 87, 033814 (2013).
  * (37) A. L. Grimsmo and Scott Parkins, Open Rabi model with ultrastrong coupling plus large dispersive-type nonlinearity: Nonclassical light via a tailored degeneracy, Phys. Rev. A 89, 033802 (2014).
  * (38) L. Cong, S. Felicetti, J. Casanova, L. Lamata, E. Solano, and I. Arrazola, Selective interactions in the quantum Rabi model, Phys. Rev. A 101, 032350 (2020).
  * (39) G. Romero, D. Ballester, Y. M. Wang, V. Scarani, and E. Solano, Ultrafast Quantum Gates in Circuit QED, Phys. Rev. Lett. 108, 120501 (2012).
  * (40) T. Wilk, S. C. Webster, H. P. Specht, G. Rempe, and A. Kuhn, Polarization-Controlled Single Photons, Phys. Rev. Lett. 98, 063601 (2007).
  * (41) J. Peng, Z. Z. Ren, H. T. Yang, G. J. Guo, X. Zhang, G. X. Ju, X. Y. Guo, C. S. Deng, and Guolin Hao, Algebraic structure of the two-qubit quantum Rabi model and its solvability using Bogoliubov operators, J. Phys. A: Math. Theor. 48, 285301 (2015).
  * (42) Y. Yin et al., Catch and Release of Microwave Photon States, Phys. Rev. Lett. 110, 107001 (2013).
  * (43) See Supplemental Material for the demonstration of the circuit design for the implementation of the two-qubit Rabi and JC model, the Lindblad master equation used for numerical simulation, the explanation of valid adiabatic evolution along $|\psi\rangle$, and the method to obtain the one-photon special dark state solution to the Rabi-Stark model.
  * (44) S. Buckley, K. Rivoire1, and J. Vučković, Engineered quantum dot single-photon sources, Rep. Prog. Phys. 75, 126503 (2012)
  * (45) M. E. Reimer and C. Cher, The quest for a perfect single-photon source, Nat. Photon. 13, 734 (2019).
  * (46) C. K. Hong, Z. Y. Ou, and L. Mandel, Measurement of subpicosecond time intervals between two photons by interference, Phys. Rev. Lett. 59, 2044 (1987).
  * (47) D. Kielpinski, J. F. Corney, and H. M. Wiseman, Quantum Optical Waveform Conversion, Phys. Rev. Lett. 106, 130501 (2011).
  * (48) A. Kiraz, M. Atatüre, and A. Imamoğlu, Quantum-dot single-photon sources: Prospects for applications in linear optics quantum-information processing, Phys. Rev. A 69, 032305 (2004).
  * (49) H. P. Eckle and H. Johannesson, A generalization of the quantum Rabi model: exact solution and spectral structure,  J. Phys. A: Math. Theor. 50, 294004 (2017).
  * (50) Y. F. Xie, L. W. Duan, and Q. H. Chen, Quantum Rabi–Stark model: solutions and exotic energy spectra, J. Phys. A: Math. Theor. 52, 245304 (2019).
  * (51) Z. M. Li and M. T. Batchelor, Hidden symmetry and tunneling dynamics in asymmetric quantum Rabi models, Phys. Rev. A 103, 023719 (2021).
  * (52) X. L. Lu, Z. M. Li, V. V. Mangazeev, and Murray T Batchelor, Hidden symmetry operators for asymmetric generalized quantum Rabi models, Chinese Phys. B 31, 014210 (2022).

###

Supplementary Material: Deterministic single photon source in the ultrastrong
coupling regime

This supplementary material contains four parts: (1) Demonstration of the
circuit design for the implementation of the two-qubit quantum Rabi and
Jaynes-Cummings (JC) model with independently tunable couplings and qubit
frequencies. (2) The Lindblad master equation used for numerical simulation.
(3) Peculiarities of the special dark state $|\psi\rangle$ for the valid and
fast adiabatic evolution along it. (4) Method to obtain the one-photon special
dark-sate solution to the Rabi-Stark model.

##  S1 Circuit design for the two-qubit quantum Rabi and JC models

Figure S1: Circuit diagram for the implementation of the two-qubit quantum
Rabi model and JC models with tunable coupling and qubit frequency.

The experimental implementation of our protocol for a single-photon source
implies the implementation of the two-qubit quantum Rabi model with tunable
coupling and tunable qubit frequency. One of the platforms that allow us to
implement tunable coupling between qubits and resonators is superconducting
circuits, where also the qubits have a tunable energy gap, and time-dependent
coupling strength can be engeniered Blais2020NatPhys ; Devoret2013Science ;
Blais2021RevModPhys . In particular, we propose the circuit design given in
Fig. S1 which is the single oscillator version of the circuit design presented
in Ref. Peng2021PhysRevLett .

This circuit is described by the Lagrangian

$\displaystyle\mathcal{L}=$
$\displaystyle\sum_{j=1}^{2}\bigg{[}\frac{C_{g}}{2}(\dot{\Phi}_{q,j}-V_{g,j})^{2}+\frac{C_{q}}{2}\dot{\Phi}_{q,j}^{2}+E_{J,j}\cos{\varphi_{q,j}}+C_{s}\dot{\Phi}_{s,j}^{2}+2E_{s}\cos{(\varphi^{(j)}_{ext})}\cos{\varphi_{s,j}}\bigg{]}+\bigg{[}\frac{C_{R}}{2}\dot{\Phi}_{C}^{2}-\frac{\Phi_{C}^{2}}{2L_{R}}\bigg{]}$
$\displaystyle+\sum_{j=1}^{2}\left[\frac{C_{c}}{2}(\dot{\Phi}_{q,j}-\dot{\Phi}_{s,j})^{2}+\frac{C_{c}}{2}(\dot{\Phi}_{C}-\dot{\Phi}_{s,j})^{2}\right]$
$\displaystyle=$
$\displaystyle\frac{1}{2}\left[\vec{\dot{\Phi}}^{\intercal}\hat{C}\vec{\dot{\Phi}}-\vec{\dot{\Phi}}^{\intercal}\hat{C}_{g}\vec{V}_{g}-\vec{V}_{g}^{\intercal}\hat{C}_{g}\vec{\dot{\Phi}}\right]-U(\vec{\Phi})$
(S1)

where

$\displaystyle\vec{\Phi}=\begin{pmatrix}\Phi_{q,1}\\\ \Phi_{s,1}\\\
\Phi_{c}\\\ \Phi_{s,2}\\\
\Phi_{q,2}\end{pmatrix},\quad\vec{V}_{g}=\begin{pmatrix}V_{g,1}\\\ 0\\\ 0\\\
0\\\ V_{g,2}\end{pmatrix},\quad\hat{C}=\begin{pmatrix}C_{qbt}&-C_{c}&0&0&0\\\
-C_{c}&C_{sqd}&-C_{c}&0&0\\\ 0&-C_{c}&C_{cav}&-C_{c}&0\\\
0&0&-C_{c}&C_{sqd}&-C_{c}\\\
0&0&0&-C_{c}&C_{qbt}\end{pmatrix},\quad\hat{C}_{g}=\begin{pmatrix}C_{g}&0&0&0&0\\\
0&0&0&0&0\\\ 0&0&0&0&0\\\ 0&0&0&0&0\\\ 0&0&0&0&C_{g}\end{pmatrix}$
$\displaystyle
U(\vec{\Phi})=\frac{\Phi_{C}^{2}}{2L_{R}}-\sum_{j=1}^{2}\left[E_{J,j}\cos{\varphi_{q,j}}+2E_{s}\cos{(\varphi^{(j)}_{ext})}\cos{\varphi_{s,j}}\right]$
(S2)

with $C_{qbt}=C_{g}+C_{q}+C_{c}$, $C_{sqd}=2(C_{s}+C_{c})$,
$C_{cav}=C_{R}+2C_{c}$, and $\varphi_{a,j}=2\pi\Phi_{a,j}/\Phi_{0}$, where
$\Phi_{0}=h/2e$ is the superconducting flux quantum and $2e$ is the electrical
charge of a Cooper pair.

To calculate the Hamiltonian of our system, we start calculating the conjugate
momenta $Q_{k}=\partial\mathcal{L}/\partial{\dot{\Phi_{k}}}$, and obtaining

$\displaystyle\vec{Q}=\hat{C}\vec{\dot{\Phi}}-\hat{C}_{g}\vec{V}_{g}$ (S3)

where

$\displaystyle\vec{Q}=\begin{pmatrix}Q_{q,1}\\\ Q_{s,1}\\\ Q_{c}\\\ Q_{s,2}\\\
Q_{q,2}\end{pmatrix}$ (S4)

Considering high plasma frequency approximation, we have that
$\dot{\Phi}_{s,j}\ll\dot{\Phi}_{q,j},\dot{\Phi}_{C}$ and
$\ddot{\Phi}_{s,j}~{}\ll~{}\ddot{\Phi}_{q,j},\ddot{\Phi}_{C}$, then the
conjugate momenta takes the form

$\displaystyle Q_{q,j}$ $\displaystyle=$ $\displaystyle
C_{qbt}\dot{\Phi}_{q,j}-Q_{g,j}\Rightarrow\dot{\Phi}_{q,j}=\frac{Q_{q,j}-Q_{g,j}}{C_{qbt}},$
$\displaystyle Q_{C}$ $\displaystyle=$ $\displaystyle
C_{cav}\dot{\Phi}_{C}\Rightarrow\dot{\Phi}_{C}=\frac{Q_{C}}{C_{cav}},$
$\displaystyle Q_{s,j}$ $\displaystyle=$ $\displaystyle-
C_{c}(\dot{\Phi}_{q,j}+\dot{\Phi}_{C})\Rightarrow
Q_{s,j}=-C_{c}\left(\frac{Q_{q,j}-Q_{g,j}}{C_{qbt}}+\frac{Q_{C}}{C_{cav}}\right).$
(S5)

Also, according to the Euler-Lagrange equations
$\frac{\partial\mathcal{L}}{\partial\Phi_{k}}-\frac{d}{dt}\frac{\partial\mathcal{L}}{\partial\dot{\Phi}_{k}}=0$,
we obtain

$\displaystyle\ddot{\Phi}_{q,j}$ $\displaystyle=$ $\displaystyle-\frac{2\pi
E_{J,j}}{C_{qbt}\phi_{0}}\sin(\varphi_{q,j})$ $\displaystyle\ddot{\Phi}_{C}$
$\displaystyle=$ $\displaystyle-\frac{1}{L_{R}C_{cav}}\Phi_{C}$
$\displaystyle\sin(\varphi_{s,j})$ $\displaystyle=$
$\displaystyle\frac{C_{c}\phi_{0}}{4\pi
E_{s}\cos(\varphi_{ext}^{(j)})}(\ddot{\Phi}_{q,j}+\ddot{\Phi}_{C})$ (S6)

where we have used the high-plasma frequency approximation. Also, we can
consider low-impedance approximation, which means that
$\sin(\varphi_{s,j})\approx\varphi_{s,j}$, then

$\varphi_{s,j}=-\frac{C_{c}\phi_{0}}{2\pi E_{s,j}^{eff}}\left[\frac{2\pi
E_{J,j}}{C_{qbt}\phi_{0}}\sin(\varphi_{q,j})+\frac{1}{L_{R}C_{cav}}\Phi_{C}\right],$
(S7)

where $E_{s,j}^{eff}=2E_{s}\cos(\varphi_{ext}^{(j)})$. By applying the
Legendre transformation
$\mathcal{H}(\vec{\Phi},\vec{Q})=\vec{Q}^{\intercal}\hat{C}^{-1}(\vec{Q}+\hat{C}_{g}\vec{V}_{g})-\mathcal{L}$,
we obtain
$\mathcal{H}(\vec{\Phi},\vec{Q})=\frac{1}{2}(\vec{Q}+\hat{C}_{g}\vec{V}_{g})^{\intercal}\hat{C}^{-1}(\vec{Q}+\hat{C}_{g}\vec{V}_{g})+U(\vec{\Phi})$,
so

$\displaystyle\mathcal{H}(\vec{\Phi},\vec{Q})=$
$\displaystyle\sum_{j=1}^{2}\left[\frac{1}{2C_{qbt}}(Q_{q,j}-2e\bar{n}_{g,j})^{2}-E_{J,j}\cos\varphi_{q,j}+\gamma_{q,j}\sin^{2}\varphi_{q,j}\right]+\frac{1}{2C_{cav}}Q_{c}^{2}+\frac{\Phi_{c}^{2}}{2L_{R}}+\gamma_{R}\Phi_{c}^{2}$
(S8)
$\displaystyle-\sum_{j=1}^{2}\left[\bar{g}_{j}\sin\varphi_{q,j}\Phi_{c}\right]$

where we have used the approximations given by Eq. (S5) and Eq. (S6), and we
define

$\displaystyle\gamma_{q,j}=\frac{C_{c}^{2}E_{J,j}^{2}}{2C_{qbt}^{2}E_{s,j}^{eff}},\quad\gamma_{R}=\frac{C_{c}^{2}\Phi_{0}^{2}}{8\pi^{2}C_{cav}^{2}L_{R}^{2}(E_{s,1}^{eff}+E_{s,2}^{eff})},\quad\bar{g}_{j}=\frac{C_{c}^{2}\Phi_{0}E_{J,j}}{2\pi
C_{cav}C_{qbt}L_{R}E_{s,j}^{eff}}$ (S9)

and also, we replace $C_{g}V_{g,j}=-2e\bar{n}_{g,j}$. We quantize the
Hamiltonian by promoting the classical variables to quantum operators,
$Q_{\alpha}\rightarrow\hat{Q}_{q,j}=2e\hat{n}_{q,j}$,
$\varphi_{q,j}\rightarrow\hat{\varphi}_{j}$, with the commutation relation
$[e^{i\hat{\varphi}_{q,j}},\hat{n}_{q,j}]=e^{i\hat{\varphi}_{q,j}}$. And for
the resonator
$Q_{c}\rightarrow\hat{Q}_{c}=\sqrt{\hbar\omega_{c}C_{cav}/2}(a^{\dagger}+a)$,
$\Phi_{c}\rightarrow\hat{\Phi}_{c}=i\sqrt{\hbar\omega_{c}L_{R}/2}(a-a^{\dagger})$.
Then our quantum Hamiltonian reads

$H=\sum_{j=1}^{2}\left[4E_{c}(\hat{n}_{q,j}-\bar{n}_{g,j})^{2}-E_{J,j}\cos\varphi_{q,j}+\gamma_{q,j}\sin^{2}\varphi_{q,j}\right]+\hbar\omega_{c}a^{\dagger}a-i\sum_{j=1}^{2}G_{j}\sin\hat{\varphi}_{q,j}(a-a^{\dagger})$
(S10)

with $E_{c}=e^{2}/2C_{qbt}$ and
$G_{j}=\bar{g}_{j}\sqrt{\frac{\hbar\omega_{c}L_{R}}{2}}$.

Using the charge number basis, i.e.,
$\hat{n}_{q,j}=\sum_{k}k|k_{j}\rangle\langle k_{j}|$ we have

$\cos\hat{\varphi}_{q,j}=\frac{1}{2}\left(\sum_{k}|k_{j}\rangle\langle
k_{j}+1|+\textrm{H.c}\right),\quad\sin\hat{\varphi}_{q,j}=-\frac{i}{2}\left(\sum_{k}|k_{j}\rangle\langle
k_{j}+1|-\textrm{H.c}\right).$ (S11)

Due to the anharmonicity of the system, we can use the two-level approximation
to obtain

$H=\sum_{j=1}^{2}\omega_{q_{j}}\sigma_{j}^{z}+\hbar\omega_{c}a^{\dagger}a-i\sum_{j=1}^{2}G_{j}\sigma_{j}^{y}(a-a^{\dagger}),$
(S12)

where $\omega_{q,j}\approx E_{J,j}$ Peng2021PhysRevLett . We use an external
flux of the form $\varphi_{ext}^{(j)}=A_{0}^{(j)}+\varphi_{AC}^{(j)}(t)$, with
$\varphi_{AC}^{(j)}(t)=A_{1}^{(j)}\cos(\mu_{1}^{(j)}t+\phi_{1}^{(j)})+A_{2}^{(j)}\cos(\mu_{2}^{(j)}t+\phi_{2}^{(j)})$.
If $|A_{1}^{(j)}|,~{}|A_{2}^{(j)}|\ll|A_{0}|$, then

$\displaystyle\frac{2}{E_{s,j}^{eff}}=\frac{1}{E_{s}\cos\varphi_{ext}^{(j)}}\approx\frac{1}{E_{s}\cos(A_{0}^{(j)})}\left[1+\frac{\sin(A_{0}^{(j)})}{\cos(A_{0}^{(j)})}\varphi_{AC}^{(j)}(t)\right]$
$\displaystyle\Rightarrow G_{j}=G_{0}^{(j)}+G_{1}^{(j)}(t),\quad
G_{0}^{(j)}=\frac{C_{c}^{2}\Phi_{0}E_{J,j}}{4\pi
C_{cav}C_{qbt}E_{s}\cos(A_{0}^{(j)})}\sqrt{\frac{\hbar\omega_{c}}{2L_{R}}},\quad
G_{1}^{(j)}(t)=G_{0}^{(j)}\frac{\sin(A_{0}^{(j)})}{\cos(A_{0}^{(j)})}\varphi_{AC}^{(j)}(t).$
(S13)

Now the Hamiltonian (S12) in the interaction picture is given by

$\displaystyle
H_{I}=\sum_{j=1}^{2}G_{0}^{(j)}\left\\{\sigma_{j}^{+}a^{\dagger}e^{i\Delta_{j}t}\left[1+\frac{\sin(A_{0}^{(j)})}{\cos(A_{0}^{(j)})}\varphi_{AC}^{(j)}(t)\right]-\sigma_{j}^{+}ae^{i\delta_{j}t}\left[1+\frac{\sin(A_{0}^{(j)})}{\cos(A_{0}^{(j)})}\varphi_{AC}^{(j)}(t)\right]\right\\}+\textrm{H.c}$
(S14)

with $\Delta_{j}=\omega_{q,j}+\omega_{c}$, and
$\delta_{j}=\omega_{q,j}-\omega_{c}$. If the qubits and the cavity are off-
resonant, then by adjusting the frequencies in $\varphi_{AC}(t)^{(j)}$, we can
activate and deactivate the different terms. In particular, choosing
$\mu_{1}^{(j)}=\Delta_{j}$ and $\mu_{2}^{(j)}=\delta_{j}$, we obtain

$\displaystyle
H_{I}\approx\sum_{j=1}^{2}\bar{G}_{0}^{(j)}\left(\sigma_{j}^{+}a^{\dagger}e^{i\phi_{1}}-\sigma_{j}^{+}ae^{i\phi_{2}}\right)+\textrm{H.c}$
(S15)

after rotating wave approximation, where
$\bar{G}_{0}^{(j)}=\frac{G_{0}^{(j)}\sin(A_{0}^{(j)})}{2\cos(A_{0}^{(j)})}$.
Finally, fixing $\phi_{1}^{(j)}=0$ and $\phi_{2}^{(j)}=\pi$, we obtain

$\displaystyle
H_{I}\approx\sum_{j=1}^{2}\bar{G}_{0}^{(j)}\sigma_{j}^{x}(a+a^{\dagger})$
(S16)

that correspond to the interaction term in the two-qubit Rabi model. We note
that $\bar{G}_{0}^{(j)}$ depends on the constant component in the external
flux, so it can be tuned independent of the qubit frequencies. The circuit
design allows us to activate and deactivate the rotating and counter-rotating
terms independently. Then, we can generate the two-qubit Rabi model activating
both terms at the same time, as well as the two-qubit JC model by the
activation only of the rotating terms.

##  S2 Lindblad master equation for numerical simulation

We use the following Lindblad form master equation

$\displaystyle\dot{\rho}=$
$\displaystyle-i[H_{pq},\rho]+\frac{\kappa}{2}(2a\rho
a^{\dagger}-a^{\dagger}a\rho-\rho
a_{i}^{\dagger}a_{i})+\sum_{j=1}^{2}\frac{\gamma_{j}}{2}(2\sigma_{j}\rho\sigma_{j}^{\dagger}-\sigma_{j}^{\dagger}\sigma_{j}\rho-\rho\sigma_{j}^{\dagger}\sigma_{j})$
(S17)
$\displaystyle+\sum_{j=1}^{2}\frac{\gamma_{j\phi}}{2}(\sigma_{jz}\rho\sigma_{jz}-\rho)$

to carry out the numerical simulation. Here, $\kappa$ is the photon decay rate
of the $i$th resonator, consisting of the intrinsic part $\kappa_{in}$ and
coupling part $\kappa_{c}$ with respect to the TL. $\gamma_{j}$ and
$\gamma_{j\phi}$ are the energy relaxation rate and the dephasing rate of the
$j$th qubit, respectively. We choose $\kappa_{in}=10^{-4}\omega$,
$\gamma_{j}=10^{-5}\omega$, $\gamma_{j\phi}=2\times 10^{-5}\omega$. Although
the ultrastrong coupling regime is reached, we use this Lindblad from master
equation because the decoherences are extremely small there. Meanwhile, the
qubit singlet state, which is decoupled from the resonator, has been generated
when $\kappa_{c}=10^{-1}\omega$ is turned on, and after the first photon is
emitted, we will enter the JC coupling regime.

##  S3 Peculiarities of the special dark state $|\psi_{R}\rangle$

There are two peculiarities of $|\psi_{R}\rangle$: 1. The energy level
$|E\rangle$ very close to $|\psi_{R}\rangle$ at Jaynes-Cummings (JC) coupling
regime in the spectrum (Fig. 1(a) in the main text) becomes degenerate with
$|\psi_{R}\rangle$ when the rotating wave approximation is applied algebra .
2. $\langle\psi_{E=\omega}|\dot{H}_{R}|\psi_{R}\rangle=0$, no matter how fast
the parameters changes, which also holds when $\dot{H}_{R}$ replaced by
$\dot{H}_{JC}$. In consequence, the adiabatic evolution along $|\psi\rangle$
is not restricted by the adiabatic theorem

$\left|\frac{\langle
E_{m}(t)|\dot{H}|E_{n}(t)\rangle}{(E_{m}-E_{n})^{2}}\right|\ll 1~{}~{}~{}m\neq
n~{}~{}~{}t\in[0,T]$ (S18)

at the degeneracy point. Numerical results show the adiabatic evolution can be
done in $68\omega^{-1}$ with fidelity $99\%$.

Here we prove $\langle\psi_{E=\omega}|\dot{H}_{R}|\psi_{R}\rangle=0$. In the
adiabatic evolution of $H_{R}$ as shown in Fig. 2(b) in the main text with
$\Delta_{1}+\Delta_{2}=\omega$, and $g_{1}=g_{2}=g$,

$\dot{H}_{R}=\dot{\Delta}_{1}\sigma_{1z}-\dot{\Delta}_{1}\sigma_{2z}+\dot{g}(a+a^{\dagger})\sigma_{1x}+\dot{g}(a+a^{\dagger})\sigma_{2x}.$
(S19)

So it is easy to find

$\dot{H}_{R}|\psi_{R}\rangle=\frac{1}{{\cal
N}}(|\downarrow,\uparrow\rangle+|\uparrow,\downarrow\rangle)\left[(\dot{g}(\Delta_{1}-\Delta_{2})-2\dot{\Delta}_{1}g)|1\rangle\right].$
(S20)

Substituting $E=\omega$ into $(H_{R}-E)|\psi_{E=\omega}\rangle=0$ and
projecting it into
$|1,\downarrow,\uparrow\rangle-|1,\uparrow,\downarrow\rangle$, we have

$\displaystyle(\Delta_{2}-\Delta_{1})\langle
1,\downarrow,\uparrow|\psi_{E=\omega}\rangle=(\Delta_{1}-\Delta_{2})\langle
1\uparrow,\downarrow|\psi_{E=\omega}\rangle.$ (S21)

Therefore,

$\displaystyle\langle 1,\downarrow,\uparrow|\psi_{E=\omega}\rangle=-\langle
1,\uparrow,\downarrow|\psi_{E=\omega}\rangle,$ (S22)

such that $\langle\psi_{E=\omega}|\dot{H}_{R}|\psi_{R}\rangle=0$ considering
Eq. (S20), no matter how fast the parameters change. This result can be easily
extended to the two-qubit JC model $H_{JC}$.

##  S4 One-photon special dark-state solution to the two-qubit quantum Rabi-
Stark model

We emphasize the deduction to obtain the one-photon dark-state solution to the
two-qubit quantum Rabi-Stark model here. A similar deduction can be found in
gaoxun . The two-qubit Rabi-Stark model reads

$\displaystyle H_{RS}=\omega
a^{{\dagger}}a+g_{1}{\sigma_{1x}}(a+{a^{\dagger}})+g_{2}{\sigma_{2x}}(a+a^{{\dagger}})+\Delta_{1}\sigma_{1z}+\Delta_{2}{\sigma_{2z}}+U_{1}a^{{\dagger}}a\sigma_{1z}+U_{2}a^{{\dagger}}a\sigma_{2z},$
(S23)

where $U_{1}$ and $U_{2}$ are couplings of the stark terms. The parity
$P=\sigma_{1z}\sigma_{2z}\exp(i\pi a^{\dagger}a)$ is still conserved here.
Supposing there is a solution with at most one photon
$|\psi_{RS}\rangle=c_{0}|0,\uparrow,\uparrow\rangle+c_{1}|0,\downarrow,\downarrow\rangle+c_{2}|1,\downarrow,\uparrow\rangle+c_{3}|1,\uparrow,\downarrow\rangle$
in the even parity subspace, then the eigenvalue equation
$(H_{RS}-E)|\psi_{RS}\rangle=0$ reads

$\displaystyle\small\left(\begin{array}[]{cccccc}\Delta_{1}+\Delta_{2}-E&0&g_{1}&g_{2}\\\
0&-\Delta_{1}-\Delta_{2}-E&g_{2}&g_{1}\\\
g_{1}&g_{2}&\omega-\Delta_{1}+\Delta_{2}-U_{1}+U_{2}-E&0\\\
g_{2}&g_{1}&0&\omega+\Delta_{1}-\Delta_{2}+U_{1}-U_{2}-E\\\
0&0&\sqrt{2}g_{1}&\sqrt{2}g_{2}\\\ 0&0&\sqrt{2}g_{2}&\sqrt{2}g_{1}\\\
\end{array}\right)\left(\begin{array}[]{c}c_{0}\\\ c_{1}\\\ c_{2}\\\ c_{3}\\\
\end{array}\right)=0.$ (S34)

There are more equations than variables in this linear equation set, which
normally does not have nontrivial solutions. However, quasi-exact solutions
exist when parameters themselves satisfy certain relations. We apply
elementary row transformations to the $6\times 4$ matrix and find the nonzero
rows can be less than the columns when $\Delta_{1}+\Delta_{2}=\omega=E$ and
$g_{1}=g_{2}=g$, with the matrix transformed into

$\begin{pmatrix}g&0&-\Delta_{1}+\Delta_{2}-U_{1}+U_{2}&0\\\ 0&1&0&0\\\
0&0&1&1\\\ 0&0&0&0\\\ 0&0&0&0\\\ 0&0&0&0\\\ \end{pmatrix},$ (S35)

which means nontrival solution

$|\psi_{RS}\rangle=(\Delta_{1}-\Delta_{2}+U_{1}-U_{2})|0\uparrow\uparrow\rangle+g|1\rangle(|\downarrow\uparrow-\uparrow\downarrow\rangle)$
(S36)

exist for arbitrary $g$ with constant energy, corresponding to a horizontal
line in the spectrum, which is a special dark state. $|\psi_{RS}\rangle$
reduces to $|\psi_{R}\rangle$ when $U_{1}=U_{2}$.

## References

  * (1) A. Blais, S. M. Girvin, and W. D. Oliver, Quantum information processing and quantum optics with circuit quantum electrodynamics, Nat. Phys. 16, 247 (2020).
  * (2) M. H. Devoret and R. J. Schoelkopf, Superconducting Circuits for Quantum Information: An Outlook, Science 339, 1169 (2013).
  * (3) A. Blais, A.L. Grimsmo, S. M. Girvin, and A. Wallraff, Circuit quantum electrodynamics, Rev. Mod. Phys. 93, 025005 (2021).
  * (4) J. Peng, J. Zheng, J. Yu, P. Tang, G. A. Barrios, J. Zhong, E. Solano, F. Albarrán-Arriagada, and L. Lamata, One-Photon Solutions to the Multiqubit Multimode Quantum Rabi Model for Fast $W$-State Generation, Phys. Rev. Lett. 127, 043604 (2021).
  * (5) J. Peng, Z. Z. Ren, H. T. Yang, G. J. Guo, X. Zhang, G. X. Ju, X. Y. Guo, C. S. Deng, and Guolin Hao, Algebraic structure of the two-qubit quantum Rabi model and its solvability using Bogoliubov operators, J. Phys. A: Math. Theor. 48, 285301 (2015).
  * (6) X. Gao, L. W. Duan, and J. Peng, Solutions to the multimode quantum Rabi model for ultrafast $W$-state generation, arXiv:2207.00775