One-photon Solutions to Multiqubit Multimode quantum Rabi model
Jie Peng, Juncong Zheng, Jing Yu, Pinghua Tang, G. Alvarado Barrios, Jianxin Zhong, Enrique Solano, F. Albarran-Arriagada, Lucas Lamata
aa r X i v : . [ qu a n t - ph ] F e b One-photon Solutions to Multiqubit Multimode quantum Rabi model
Jie Peng, ∗ Juncong Zheng, Jing Yu, Pinghua Tang, G. Alvarado Barrios, JianxinZhong, Enrique Solano,
2, 3, 4, 5, † F. Albarr´an-Arriagada, and Lucas Lamata ‡ Hunan Key Laboratory for Micro-Nano Energy Materials and Devices and Schoolof Physics and Optoelectronics, Xiangtan University, Hunan 411105, China International Center of Quantum Artificial Intelligence for Science and Technology (QuArtist)and Physics Department, Shanghai University, 200444 Shanghai, China Department of Physical Chemistry, University of the Basque Country UPV/EHU, Apartado 644, 48080 Bilbao, Spain IKERBASQUE, Basque Foundation for Science, Plaza Euskadi 5, 48009, Spain IQM, Nymphenburgerstr. 86, 80636 Munich, Germany Departamento de F´ısica At´omica, Molecular y Nuclear, Universidad de Sevilla, 41080 Sevilla, Spain
General solutions to the quantum Rabi model involve subspaces with unbounded number of pho-tons. However, for the multiqubit multimode case, we find special solutions with at most onephoton for arbitrary number of qubits and photon modes. Unlike the Juddian solution, ours existsfor arbitrary single qubit-photon coupling strength with constant eigenenergy. This corresponds toa horizontal line in the spectrum, while still being a qubit-photon entangled state. As a possibleapplication, we propose an adiabatic scheme for the fast generation of arbitrary single-photon multi-mode W states with nonadiabatic error less than 1%. Finally, we propose a superconducting circuitdesign, showing the experimental feasibility of the multimode multiqubit Rabi model.
Introduction.–
The quantum Rabi model [1, 2] de-scribes the interaction between a two-level system anda single photonic mode at the most fundamental level.Despite its simple form, the exact solution of the modelwas not found until 2011 [3]. Since the quantum Rabimodel involves both rotating and counter-rotating inter-action terms, all Fock states are connected and there is noclosed subspace, turning the Hamiltonian hard to solve.The quantum Rabi model plays an important role inquantum optics [4–7], molecular physics [8], condensedmatter physics [9] and quantum information [10, 11].However, in most applications we need to consider morethan one qubit and/or more than one mode. For in-stance, to perform a controlled gate [12], essential foruniversal quantum computing [13, 14], while multimodemodels are useful for the generation of the multipartiteentangled states [15, 16]. Therefore, a mathematical de-scription with physical implications for models with morethan one qubit and one mode is essential for the devel-opment of scalable and efficient protocols, suitable forcurrent technology demands.In this article, we find special solutions with at mostone photon to the multiqubit multimode quantum Rabimodel (MMQRM) for arbitrary number of qubits andmodes, although the interaction terms still connect allphoton number states. Unlike the well known Jud-dian solution [17] of the quantum Rabi model, thesequasi-exact solutions exist for arbitrary single qubit-photon coupling strength with constant eigenenergy, cor-responding to a horizontal line in the spectrum, whilestill being a qubit-photon entangled state. This co-herent superposition is what makes the photon popu- ∗ [email protected] † [email protected] ‡ [email protected] lation trapped in zero and one, and we call it a spe-cial dark state [18]. Furthermore, we use such solu-tions to propose a fast entangled state generation pro-tocol, where we simultaneously obtain a two-qubit Bellstate and an arbitrary single-photon M − mode W state | W i M = P Mi =1 g i | · · · i i +1 · · · M i [19–21] throughadiabatic passage. Here, g i just corresponds to the cou-pling strength between the qubits and the i -th photonmode. It is known that W states are robust under par-ticle loss [22] and a central resource in several quantuminformation processing protocols [23–27]. Consequently,various schemes have been presented to generate them[28–36]. Due to the reach of ultrastrong coupling and pe-culiarities of the special dark states, the most interestingadvantage of using the MMQRM is the fast generation(less than 70 ω − ) and high fidelity (larger than 99%),outperforming the fastest two-qubit CPHASE gate (30-45 ns) to date [37] for ω = 3 GHz. The generation timewill not change with the number of modes M and noexternal laser is needed. Since g i is adjustable, we cangenerate any W state in a unified and convenient way,such as the perfect W states which are useful in quantumteleportation [23]. Finally, we propose a superconduct-ing circuit design to show the experimental feasibility ofcatch and release of these W states. This result paves theway to the implementation of fast protocols in quantuminformation using the MMQRM. Special quasi-exact solutions to the multiqubit multi-mode quantum Rabi model.–
We present our method toobtain the quasi-exact solution with at most one photonof the multiqubit and multimode quantum Rabi model H pq = M X i =1 ω i a † i a i + M X i =1 N X j =1 g ij σ jx ( a i + a † i ) + N X j =1 ∆ j σ jz , (1)where a † i and a i are the i -th photon mode creation andannihilation operators with frequency ω i , respectively.Also, σ jα ( α = x, y, z ) are the Pauli matrices correspond-ing to the j -th qubit, 2∆ j is the energy level splittingof the j -th qubit, and g ij is the qubit-photon couplingparameter between the i -th mode and j -th qubit.Since Hamiltonian (1) breaks the U (1) symmetry, thereis no closed subspace consisting of finite photon numberstates. However, it has a Z symmetry with generator R = exp[ iπ P Mi =1 a † i a i ]Π j σ jz . Accordingly, we categorizeall N-qubit states {| ψ Nq i} into two sets correspondingto the eigenvuales of Π j σ jz , being 1 and −
1, and denotethem 2 N − dimensional row vectors, | ψ Nq + i and | ψ Nq − i ,respectively. We also denote all k-photon states with Mmodes by | k M i . Hence, there are two invariant subspaces | M , ψ Nq + i ↔ | M , ψ Nq − i ↔ | M , ψ Nq + i ↔ · · · (2) | M , ψ Nq − i ↔ | M , ψ Nq + i ↔ | M , ψ Nq − i ↔ · · · (3)with positive and negative parity, respectively. So H pq will take the following form in the ± parity subspace H ± pq = D ± O ± . . .O ± D ± O ± . . . O ± D ± O ± . . .. . . . . . . . . . . . . . . . . . , (4)where D ± k and O ± k ( k = 0 , , , , . . . ) are2 N − C kM + k − × N − C kM + k − matrices [38]. Althoughall Fock states are connected by the interaction terms,it is possible to find some special solutions with finitephoton number. Indeed, there are special dark statesin the single-mode multiqubit Rabi model consisting offinite Fock states [39–41].We search now for solutions with at most L pho-tons taking the form | ψ ± i = c ± M | M , ψ Nq, ± i + c ± M | M , ψ Nq, ± i + . . . + c ± LM | L M , ψ Nq, ± i in the MMQRMby solving the energy eigenequation( H ± pq − E ± ) | ψ ± i = 0 . (5)Substituting H ± pq from Eq. (4) into Eq. (5) and solv-ing for the coefficients c ± k M , we see there are more equa-tions than variables ( c ± k M ) provided we do not trun-cate H pq . Nevertheless, solutions could exist when theHamiltonian parameters themselves meet certain condi-tion to satisfy Eq. (5). Indeed, we have proved that for L > f ( g ij , ω i , ∆ j ) = 0 [38], just like the Juddian solution forthe single-qubit quantum Rabi model. However, we finda remarkable case L = 1, where the condition reduces to f ( ω i , ∆ j ) = 0 and f ( g ij ) = 0, so the solutions changefrom single points into horizontal lines in the spectrawhen ω i and ∆ j are fixed. The latter may have im-portant applications in quantum information, as will bediscussed below. Therefore, we will focus on the case L = 1, where Eq. (5) reduces to D ± − E ± O ± D ± − E ± O (cid:18) c ± M c ± M (cid:19) = 0 (6) Figure
1. (a) Spectrum of the two-qubit two-mode quantumRabi model with ω = ω = ω , ∆ = 0 . ω , ∆ = 0 . ω , g = g = g = g = g . (b) Spectrum of the three-qubittwo-mode quantum Rabi model with ∆ = ∆ = ∆ = ω = ω = ω , g = g + g = g , g = g = g = g . Redlines correspond to even parity while blue lines to odd parity. after elementary row matrix transformation. For the sin-gle qubit and single mode case, O is just a c-number, sothere is no nontrivial solution for O ± c ± = 0. However,for the multiqubit and multimode case, this equation canbe satisfied if matrix O ± has eigenvalue 0 with c ± M be-ing its corresponding eigenvector. Then, we make theelementary row transformation to the matrix in Eq. (6),and the solution is obtained if there are more columnsthan nonzero rows. At the same time, the condition forthe parameters is calculated.For the two-qubit and M-mode case, we find the specialsolution for even parity to be (see [38] for details) | ψ i = 1 N [(∆ − ∆ ) | M , ↑ , ↑i + | W M i ( | ↓ , ↑i − | ↑ , ↓i )](7)where | W M i = g | , , , . . . , i + g | , , , . . . , i + · · · + g M | , , , . . . , i with the condition ω i = ω for all i , g ij = g i for all j and ∆ + ∆ = ω = E + . This spe-cial solution has some novel properties: (1) It exists forarbitrary g i , with constant eigenenergy E + = ω , corre-sponding to a horizontal line in the spectrum, as shown inFig. 1 (a). (2) It is a special dark state because althoughthe interaction terms connect all Fock states, the popula-tion is trapped in vacuum and single photon multimodeW states [18]. There are two other similar dark-state so-lutions corresponding to odd parity which are describedin the supplementary material [38].For the three-qubit and M-mode case, the quasi-exactsolution for odd parity reads | ψ − i = | W M i ( | ↑ , ↓ , ↓i − | ↓ , ↑ , ↓i − | ↓ , ↓ , ↑i + | ↑ , ↑ , ↑i )+ ωg g | M , ↑ , ↑ , ↓i + ωg g | M , ↑ , ↓ , ↑i− ωg g g | M , ↓ , ↑ , ↑i , (8)where | W M i = g | , , , . . . , i + g | , , , . . . , i + · · · + g M | , , , . . . , i with the conditions ∆ j = ω i = ω = E − , g i = g i + g i . This eigenstate corresponds to thehorizontal line E/ω = 1 in Fig. 1 (b).As the qubit number grows, the existence conditionwill be harder to satisfy, but there are indeed quasi-exactsolutions for the N-qubit M-mode Rabi model formed bythe product of the two-qubit singlet Bell state | ψ B i = √ ( | ↓↑i − | ↑↓i ) and | ψ i or | ψ i , | ψ N i = | ψ i ⊗ ( | ψ B i ) ( N − / , (9) | ψ N i = | ψ i ⊗ ( | ψ B i ) ( N − / , (10)for even and odd N , respectively. Fast Generation of the arbitrary single-photon multi-mode W state.–
Here, the special dark state solution forthe two-qubit and M-mode quantum Rabi model | ψ i inEq. (7) is an excellent candidate for generating arbitraryW states | W M i through adiabatic passage: (1) The en-ergy gap limiting the adiabatic speed can be much largerthan normal cases due to its peculiarities, which will bedetailed later. (2) It consists of only | M i and | W M i forthe photon part, and the coefficient of the former is pro-portional to ∆ − ∆ , so once the qubit frequencies aretuned to be equal, we immediately arrive at | W M i forany nonzero g i , robust against its fluctuation. (3) Thecoefficient of | · · · i i +1 · · · i is just g i in | W M i , soit is very convenient to construct any W state with arbi-trary modes by adjusting g i directly, and the generationtime will stay the same. (4) At the beginning of theadiabatic passage, the initial state | M , ↑ , ↑i is easy toprepare, while in the end, we arrive at a W state | W M i and a qubit Bell state | ψ B i simultaneously, both veryuseful in quantum information processing, while no ex-ternal laser pulse is needed. Thereafter, the W state isautomatically stored in the resonator because the photonfield and qubit Bell state are decoupled.Our scheme is as follows. First, two qubits are excitedby pumping pulses and coupled to M resonators in vac-uum states with initial coupling strength g i = g i = g i = 0. The qubit frequencies are non-identical and al-ways satisfy ∆ +∆ = ω i = ω . Then, we slowly decrease | ∆ − ∆ | to 0 while increase g i to a nonzero value, so thatthe target state | W M ψ B i is obtained. As an example, thenumerical simulation for the adiabatic evolution of thetwo-qubit two-mode case is shown in Fig. 2 (c) with theadiabatic trajectory used to vary the parameters shownin Fig. 2 (a). The evolution time is just T = 100 ω − and the fidelity F = |h ψ ( T ) | W ψ B i| reaches 99 . Figure
2. Adiabatic evolution of the special dark state | ψ i of Eq. (7) from | M ↑↑i to | W M ψ B i for the two-qubit M-mode quantum Rabi model, where | ψ B i = √ ( | ↓↑ − ↑↓i ),and | W M i is the prototype W state with all g i ’s equal. (a)Adiabatic trajectory used to vary the parameters for M = 2. g = g = g . (b) Fidelity F M = |h ψ ( T ) | W M i| when T is fixedto 100 ω − for mode number M . (c) Population of differentstates during the adiabatic process for M = 2. (d) Interactiontime T M to reach F M > .
99 for mode M in unit of ω − . the evolution time is fixed to 100 ω − , the fidelities F M for the M-mode case are shown in Fig. 2 (b), which arealmost equal and higher than 99 . ω − .According to the current available circuit QED tech-nology, the transmon frequency ∆ /π can be tuned from0 to 6 GHz [42, 43], hence the resonator frequency ω/ π is chosen to be 3 GHz [44] to satisfy ∆ ± ∆ = ω at g/ π = 0. The maximum value of g/ π can be tunedto 0.7535 GHz [44], that is, 0 . ω/ π . Therefore, ourscheme shown in Fig. 2 (a) is within experimental reachand the adiabatic evolution takes only 33.3 ns with nona-diabatic error 0 . | W M ψ B i from | M ↑↑i for M = 2 , , . . . ,
10 will be 21.9 ns. Note that we have cho-sen the simplest linear adiabatic path shown in Fig. 2 (a),but we can also consider a “faster adiabatic” trajectoryto reduce the time such as in Refs. [45, 47, 48].Although we did not optimize the adiabatic path, thesimulated adiabatic evolution speed had already beenfaster than or similar to the optimized ones [45, 47, 48].This is due to the reach of ultrastrong coupling andthe peculiarities of the special dark state | ψ i [38].(1) There are C M +1 + 1 degenerate eigenstates | ψ E = ω i including | ψ i at Jaynes-Cummings coupling regimewhere the rotating wave approximation is applied. (2) h ψ E = ω | ˙ H | ψ i = 0, no matter how fast the parameterschange. As can be seen in the spectrum in Fig. 1 (a), Figure
3. (a) Schematic setup for the generation and re-lease of the W state: Two SQs are capacitively coupled to MCWRs. Each CWR is connected to a TL through a variablecoupler C, such that the photon emission rate into the TLis controllable. (b) A superconducting circuit design for thetwo-qubit two-mode Rabi model. there are three energy levels very close to the horizon-tal line dark state | ψ i at E = ω when g is small.Actually, they correspond to degenerate eigenstates of H pq with E = ω when rotating wave approximation isapplied [38], which is valid for small g , and we proved h ψ E = ω | ˙ H | ψ i = 0 [38], such that the actual energy gaplimiting the adiabatic speed according to the adiabatictheorem [49, 50] (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) h E m ( t ) | ˙ H | E n ( t ) i ( E m − E n ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≪ , m = n, t ∈ [0 , T ] . (11)is about | ∆ − ∆ | at Jaynes-Cummings coupling regime,which can be tuned to be ω , much larger than normalcases. For the same reason, the adiabatic speed is notlimited by the vanishing energy gap at the degeneracypoints around g ≈ .
357 and 0 . = ∆ and g = 0, we find the adiabatic evolution will be fasterif g is increased to the ultrastrong coupling regime. Catch and release of the W state. – On the other hand,just generating a W state inside resonators is not con-venient for its transport and detection, such that wepropose a scheme to store or extract them out on de-mand for practical usage in quantum teleportation andrelated tasks. Our scheme is depicted in Fig. 3 (a),where two superconducting qubits (SQs) in the center ofthe devices are capacitively coupled to N coplanar waveg-uide resonators (CWRs), which can be described by thetwo-qubit M-mode Rabi model. A detailed descriptionof the superconducting circuit design for the two-qubittwo-mode Rabi model (see Fig. 3 (b)) is shown in [38].Besides, there is an externally variable coupler to modu-late the decay rate κ c of each CWR through that coupler,such that its photon emission into the connected trans-mission line (TL) is controllable. In current experimentalsetups, κ c can be tuned to be 1000 times the CWR intrin- Figure
4. Numerical simulation for the catch and releaseof the prototype three-mode W state | W i = √ ( | i + | i + | i ) (left panel) and the four-mode W state | W ′ i = √ ( | i + √ | i + √ | i + 2 | i ) (right panel).(a) and (b) Population of different states inside resonators.(c) and (d) Emission rates into transmission lines. sic decay rate κ in [51], such that we can catch (generateand store) and release the W state on demand.At the first step, an arbitrary single photon M-modeW state and a two-qubit Bell state can be generated si-multaneously using the scheme discussed above. There-after, the qubit Bell state is decoupled from the CWRs,hence | W M i is naturally stored in resonators, and robustagainst any fluctuation on g i . After a desired time τ ,we turn on the coupling κ c and therefore the W state isreleased into the transmission lines. For realistic experi-mental considerations, we choose the intrinsic dissipationrate for each resonator κ in = 10 − ω and κ c = 10 − ω .The energy relaxation rate for each qubit γ j = 10 − ω and the dephasing rate γ jφ = 10 − ω . Using a Lind-blad equation [38], we present the numerical simulationin Fig. 4. The time cost for generating | W i and | W ′ i is still 100 ω − , but the fidelity is reduced to 98% due todissipation, damping and dephasing.However, the prototype W state | W i = | i + | i + | i generated in Fig. 4 (c) cannot be used to per-form quantum teleportation, while the so-called per-fect W state | W ′ i = ( | i + | i + √ | i ) can[23]. This state can be easily generated just by tuning g = √ g = √ g , as shown in Fig. 5 (a). The first twoqubits are kept by Alice while the third one is sent toBob to carry out the remote communication. If we knowthe distance L between them, we can delay the emissionof the W state into the first two transmission lines by L/C to assure they receive the qubits at the same time,which will increase the communication security. The cor-responding numerical simulation is shown in Fig. 5 (b).
Conclusions.–
All Fock states are excited by the dipoleinteraction in the ultrastrong coupling regime of thequantum Rabi model, where the rotating-wave approx-imation is not valid. Therefore, although the operation
Figure
5. Numerical simulation for the generation and con-trolled release of the perfect three-mode W state | W ′ i = ( | i + | i + √ | i ) for quantum teleportation. (a)Population of different states inside resonators. (b) Emissionrates into transmission lines. can be faster, it seems impossible to construct any kindof single photon state. However, we find that for themultiqubit multimode quantum Rabi model, there ex-ist special dark eigenstates consisting of only vacuumand single photon multimode W states for the photonpart in the whole coupling regime with constant energy.Accordingly, we propose a unified scheme to adiabati-cally generate arbitrary W states using the special dark-state solution to the two-qubit M-mode quantum Rabi model, being able to take advantages of the ultrastrongcoupling and avoid its dynamical complexities. Due tothe peculiarities of these states, the energy gap limitingthe speed can be tuned to be much larger than normalcases. Hence, the time cost according to the current cir-cuit QED technology (33ns with nonadiabatic error 0 . | , , . . . , i , . . . i in the W state is just proportional tothe coupling strength between the qubits and the i -thresonator, such that any W state can be obtained justby tuning g i . Moreover, it can also be released into thetransmission lines on demand, which are illustrated to beuseful in quantum information processing. Similar usesof other special dark states still need to be explored. Acknowledgements.–
This work was supported bythe National Natural Science Foundation of China(11704320), Natural Science Foundation of HunanProvince, China ( 2018JJ3482), the National BasicResearch Program of China (2015CB921103), the Pro-gram for Changjiang Scholars and Innovative ResearchTeam in University (No. IRT13093), the fundingfrom PGC2018-095113-B-I00, PID2019-104002GB-C21, Spanish Government PID2019-104002GB-C22(MCIU/AEI/FEDER, UE), Spanish GovernmentPGC2018-095113-B-I00 (MCIU/AEI/FEDER, UE),Basque Government IT986-16, as well as from QMiCS(820505) and OpenSuperQ (820363) of the EU Flagshipon Quantum Technologies, EU FET Open Grant Quro-morphic (828826), EPIQUS (899368) and Shang- haiSTCSM (Grant No. 2019SHZDZX01-ZX04). [1] I. I. Rabi, Phys. Rev. , 324 (1936); , 652 (1937).[2] E. T. Jaynes and F.W. Cummings, Proc. IEEE , 89(1963).[3] D. Braak, Phys. Rev. Lett. , 100401 (2011).[4] T. Werlang, A. V. Dodonov, E. I. Duzzioni, and C. J.Villas-Bˆoas, Phys. Rev. A , 053805 (2008).[5] J. Casanova, G. Romero, I. Lizuain, J. J. Garc´ıa-Ripoll,and E. Solano, Phys. Rev. Lett. , 263603 (2010).[6] J. S. Pedernales, I. Lizuain, S. Felicetti, G. Romero, L.Lamata, and E. Solano, Sci. Rep. , 15472 (2015).[7] L. Lamata, Sci. Rep. , 43768 (2017).[8] I. Thanopulos, E. Paspalakis, and Z. Kis, Chem. Phys.Lett. , 228 (2004).[9] E. K. Irish, Phys. Rev. Lett. , 173601 (2007).[10] A. Blais, R.-S. Huang, A. Wallraff, S. M. Girvin and R.J. Schoelkopf, Phys. Rev. A , 062320 (2004).[11] X.-Y. L¨u, G.-L. Zhu, L.-L. Zheng, and Y. Wu, Phys. Rev.A , 033807 (2018).[12] G. Romero, D. Ballester, Y. M. Wang, V. Scarani, andE. Solano, Phys. Rev. Lett. , 120501 (2012).[13] A. Barenco et al. , Phys. Rev. A , 3457 (1995).[14] M. Hua, M.-J. Tao, and F.-G. Deng, Phys. Rev. A , 012328 (2014).[15] Y.-H. Kang, Y.-H. Chen, Q.-C. Wu, B.-H. Huang, J.Song, and Y. Xia, Sci. Rep. , 36737 (2016).[16] M. Lu, Y. Xia, J. Song, and N. Ba An, J. Opt. Soc. Am.B , 2142 (2013).[17] B. R. Judd, J. Phys. C , 1685 (1979).[18] M. O. Scully and M. S. Zubairy, Quantum Optics(Cambridge University Press, Cambridge, England, 1997),Chap. 7.[19] L. Zhou, Y.-B. Sheng, W.-W. Cheng, L.-Y. Gong, andS.-M. Zhao, J. Opt. Soc. Am. B , 71 (2013).[20] Y.-B. Sheng, Y. Ou-Yang, L. Zhou and L. Wang, Quan-tum Inf. Process. , 1595 (2014).[21] A. Sharma, and A. A. Tulapurkar, Phys. Rev. A ,062330 (2020).[22] N. Kiesel, M. Bourennane, C. Kurtsiefer, H. Weinfurter,D. Kaszlikowski, W. Laskowski and M. Zukowski, J. Mod.Opt. , 1131 (2003).[23] P. Agrawal, A. Pati, Phys. Rev. A , 062320 (2006).[24] S.-B. Zheng, Phys. Rev. A , 054303 (2006).[25] L. Li and D. Qiu, J. Phys. A: Math. Theor. , 10871(2007). [26] A. K. Ekert, Phys. Rev. Lett. , 661 (1991).[27] L.-H. Yan, Y.-F. Gao, and J.-G. Zhao, Int. J. Theor.Phys. , 2445 (2009).[28] M. Eibl, N. Kiesel, M. Bourennane, C. Kurtsiefer, andH. Weinfurter, Phys. Rev. Lett. , 077901 (2004).[29] M. Menotti, L. Maccone, J. E. Sipe, and M. Liscidini,Phys. Rev. A , 013845 (2016).[30] B. Fang, M. Menotti, M. Liscidini, J. E. Sipe, and V. O.Lorenz, Phys. Rev. Lett. , 070508 (2019).[31] X. B. Zou, K. Pahlke, and W. Mathis, Phys. Rev. A ,044302 (2002).[32] Guo-Ping Guo, Chuan-Feng Li, Jian Li, Guang-Can Guo,Phys. Rev. A , 042102 (2002).[33] Z. J. Deng, M. Feng, and K. L. Gao, Phys. Rev. A ,014302 (2006).[34] Guang-Can Guo and Yong-Sheng Zhang, Phys. Rev. A , 042302 (2002).[35] Y.-H. Kang, Y.-H. Chen, Q.-C. Wu, B.-H. Huang, J.Song, and Y. Xia, Sci. Rep. , 36737 (2016).[36] V. M. Stojanovi´c, Phys. Rev. Lett. , 190504 (2020).[37] M. A. Rol et al., Phys. Rev. Lett. , 120502 (2019).[38] See Supplementary Material for additional details.[39] J. Peng, Z. Z. Ren, D. Braak, G. J. Guo, G. X. Ju, X.Zhang, and X. Y. Guo, J. Phys. A: Math. Theor. , 265303 (2014).[40] J. Peng et al., J. Phys. A: Math. Theor. , 285301(2015).[41] J. Peng et al., J. Phys. A: Math. Theor. , 174003(2017).[42] M. D. Hutchings et al., Phys. Rev. App. , 044003 (2017).[43] P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando, S.Gustavsson, and W. D. Oliver, Appl. Phys. Rev. , 021318(2019).[44] A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff,arXiv:2005.12667[45] R. Barends, J. Kelly, A. Megrant, A. Veitia, D. Sank, etal., Nature , 500 (2014).[46] M. Kjaergaard, M. E. Schwartz, J. Braum¨uller, P.Krantz, J. I.-J.Wang, S. Gustavsson, and W. D. Oliver,Annu. Rev. Condens. Matter Phys. , 369 (2020).[47] J. M. Martinis, and M. R. Geller, Phys. Rev. A ,022307 (2014).[48] Y. Chen, C. Neill, P. Roushan, N. Leung, et al., Phys.Rev. Lett. , 220502 (2014).[49] P. Ehrenfest, Ann. Phys. (Leipzig) , 327 (1916).[50] M. Born and V. Fock, Z. Phys. , 165 (1928).[51] Y. Yin et al., Phys. Rev. Lett. , 107001 (2013). Supplementary Material: One-photon Solutions to Multiqubit Multimode quantumRabi model
This supplementary material contains four parts: (1) Quasi-exact solution to the multiqubit multimode Rabi model;(2) Peculiarities of the special dark state | ψ i ; (3) Demonstration of the circuit design for the implementation ofthe two-qubit two-mode quantum Rabi model with variable couplings; (4) The Lindblad master equation we used fornumerical simulation. S1. QUASI-EXACT SOLUTION TO THE MULTIQUBIT MULTIMODE QUANTUM RABI MODEL D ± k and O ± k ( k = 0 , , , , . . . ) can be written as D ± k = ( h k M , ψ Nq, ± ( − k | ) T H pq ( | k M , ψ Nq, ± ( − k i ) , (S1) O ± k = ( h ( k + 1) M , ψ Nq, ∓ ( − k | ) T H pq ( | k M , ψ Nq, ± ( − k i ) , (S2)where ( | k M , ψ Nq, ± ( − j i ) is a 2 N − C kM + k − dimensional row vector, consisting of all possible qubit-photon productstate with k photons , N qubits and ± parity. Hence D ± k , O ± k are 2 N − C kM + k − × N − C kM + k − matrices. If eigenstate | ψ ± i = c ± M | M , ψ Nq, ± i + c ± M | M , ψ Nq, ± i + . . . + c ± LM | L M , ψ Nq, ± i exists, then according to Eq. (5), we obtain D ± − E ± O . . .O D ± − E ± O . . .. . . . . . . . . . . . . . . . . . O L − D ± L − − E ± O ± L − . . . O L − D ± L − E ± . . . O ± L c ± ,M c ± ,M . . .c ± L − ,M c ± L,M = 0 . (S3)Clearly, there are more equations than variables, but quasi-exact solution is possible if the parameters meet certaincondition. A necessary but not sufficient condition for Eq. (S3) is (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) D ± − E ± O . . .O D ± − E ± O . . .. . . . . . . . . . . . . . . . . . O L − D ± L − − E ± O L − . . . O L − D ± L − E ± (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 0 . (S4)Consequently, in general, E ± is dependent on the couplings. But if L = 1, Eq. (S3) reduces to D ± − E ± O ± D ± − E ± O (cid:18) c ± M c ± M (cid:19) = 0 , (S5)after elementary row transformations, which will be demonstrated later, such that the determinant in Eq. (S4) reducesto (cid:12)(cid:12)(cid:12)(cid:12) D ± − E ± O D ± − E ± (cid:12)(cid:12)(cid:12)(cid:12) = 0 . (S6)Here, E ± is independent of the couplings and its existence condition reduces to f ( ω i , ∆ j ) = 0 and f ( g ij ) = 0, whichwill be much easier to realize than f ( ω i , ∆ j , g ij ) = 0 generally. This character may have important applications inquantum information processing, such that we focus on the case L = 1.Let us start with the simplest case, the two-qubit two-mode quantum Rabi model, H q m = ω a † a + ω a † a + ( g σ x + g σ x )( a + a † ) + ( g σ x + g σ x )( a + a † ) + ∆ σ z + ∆ σ z . (S7)For even parity, in the basis formed by {| , , ↑ , ↑i , | , , ↓ , ↓i , | , , ↑ , ↓i , | , , ↓ , ↑i , | , , ↑ , ↓i , | , , ↓ , ↑i , | , , ↑ , ↑i , | , , ↓ , ↓i , | , , ↑ , ↑i , | , , ↓ , ↓i , | , , ↑ , ↑i , | , , ↓ , ↓i} , the coefficient matrix of Eq. (S3) reads − ∆ − ∆ − E + g g g g + ∆ − E + g g g g g g ω + ∆ − ∆ − E + g g ω − ∆ + ∆ − E + g g ω + ∆ − ∆ − E + g g ω − ∆ + ∆ − E + √ g √ g √ g √ g g g g g g g g g √ g √ g √ g √ g . Obviously, it takes the form of Eq. (S5) after elementary row transformations. Nontrivial solutions exist if thereare less nonzero rows than columns, which will be satisfied when ω = ω = E + = ∆ + ∆ , g = g = g , and g = g = g . The solution reads | ψ i = 1 N [(∆ − ∆ ) | , , ↑ , ↑i + ( g | , i + g | , i )( | ↓ , ↑i − | ↑ , ↓i )] . (S8)Extending our analysis to the M-mode case, it is easy to find the special dark state solution | ψ i (7) in the maintext. For odd parity, there are similar special dark states | ψ i − ,a = 1 N ′ [(∆ + ∆ ) | M , ↑ , ↓i + | W M i ( | ↓ , ↓i − | ↑ , ↑i )] , (S9) | ψ i − ,b = 1 N ′ [(∆ + ∆ ) | M , ↓ , ↑i + | W M i| ↓ , ↓i − | ↑ , ↑i ] , (S10)with the condition ∆ − ∆ = ω = E − and ∆ − ∆ = ω = E − respectively, and the condition for ω i , g ij is the sameas the even parity. Using a similar method, we can obtain the special dark state solution in Eq. (8) for the three-qubitM-mode quantum Rabi model. S2. PECULIARITIES OF THE SPECIAL DARK STATE | ψ i There are two peculiarities of | ψ i : 1. h ψ E = ω | ˙ H | ψ i = 0, no matter how fast the parameters changes. 2.There are C M +1 + 1 degenerate eigenstates | ψ E = ω i at Jaynes-Cummings coupling regime where the rotating waveapproximation is applied. Here we give a proof.First, we prove h ψ E = ω | ˙ H | ψ i = 0. In the adiabatic evolution of the two-qubit and two-mode quantum Rabi modelwith ∆ + ∆ = ω = ω = ω , g = g = g and g = g = g ,˙ H R = ˙∆ σ z − ˙∆ σ z + ˙ g ( a + a † ) σ x + ˙ g ( a + a † ) σ x + ˙ g ( a + a † ) σ x + ˙ g ( a + a † ) σ x . (S11)So it is easy to find˙ H R | ψ i = 1 N ( | ↓ , ↑i + | ↑ , ↓i ) h ( ˙ g (∆ − ∆ ) − g ) | , i + ( ˙ g (∆ − ∆ ) − g ) | , i i . (S12)Substituting E = ω into ( H R − E ) | ψ E = ω i = 0, we have(∆ − ∆ ) h , , ↓ , ↑ | ψ E = ω i = (∆ − ∆ ) h , , ↑ , ↓ | ψ E = ω i (∆ − ∆ ) h , , ↓ , ↑ | ψ E = ω i = (∆ − ∆ ) h , , ↑ , ↓ | ψ E = ω i . (S13)Therefore, h , , ↓ , ↑ | ψ E = ω i = −h , , ↑ , ↓ | ψ E = ω i , h , , ↓ , ↑ | ψ E = ω i = −h , , ↑ , ↓ | ψ E = ω i , (S14)such that h ψ E = ω | ˙ H R | ψ i = 0, no matter how fast the parameters change. In consequence, the adiabatic speed isnot restricted by the adiabatic theorem (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) h E m ( t ) | ˙ H | E n ( t ) i ( E m − E n ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≪ m = n t ∈ [0 , T ] (S15)at the degeneracy point. This result can be easily extended to the two-qubit and M-mode case.Next, we prove there are four degenerate eigenstates for the two-qubit two-mode Jaynes-Cummings (JC) model H J = ω a † a + ω a † a + g ( a σ +1 + a † σ )+ g ( a σ +2 + a † σ )+ g ( a σ +1 + a † σ )+ g ( a σ +2 + a † σ )+ ∆ σ z + ∆ σ z , (S16)when ∆ + ∆ = ω = ω , g = g = g and g = g = g . For this model, the excitation number operator C = a † a + a † a + ( σ z + σ z ) / C = 2 consisting of {| , , ↑ , ↑i , | , , ↑ , ↓i , | , , ↓ , ↑i , | , , ↑ , ↓i , | , , ↓ , ↑i , | , , ↓ , ↓i , | , , ↓ , ↓i , | , , ↓ , ↓i} , the Hamiltonian reads ω g g g g g ω + ∆ − ∆ √ g g g ω − ∆ + ∆ √ g g g ω + ∆ − ∆ g √ g g ω − ∆ + ∆ g √ g √ g √ g ω g g g g ω
00 0 0 √ g √ g ω , (S17)and the secular equation | H J − E | = 0 takes the form( E − ω ) f ( E, ∆ , ∆ , g , g , ω ) = 0 . (S18)We have then four degenerate eigenstates | ψ E = ω i with E = ω , one of which can be | ψ i in Eq. (7). It is easy tofind that Eqs. (S12) and (S14) still hold for H J , so we still obtain h ψ E = ω | ˙ H J | ψ i = 0. Therefore, the energy gaplimiting the adiabatic speed according to the adiabatic theorem will only be dependent on the other four energy levels,that is, | ∆ − ∆ | , which can be tuned to be ω when g ∼
0. This result can be easily extended to the two-qubit M-modequantum Rabi model, and that is one reason why the adiabatic generation of the multimode W state is exceptionallyfast. Since there are C M +1 states | M , ↓ , ↓i which has the same energy ω as | M , ↑ , ↑i when ∆ + ∆ = ω and g ij = g i ,it is easy to prove there are C M +1 + 1 degenerate eigenstates with E = ω . S3. DEMONSTRATION OF THE CIRCUIT DESIGN FOR THE IMPLEMENTATION OF THETWO-QUBIT TWO-MODE QUANTUM RABI MODEL WITH VARIABLE COUPLING
For the physical implementation of the two-qubit two-mode quantum Rabi model, we propose the circuit given byFig. 3 (b), which is described by the Lagrangian L = X j =1 (cid:20) C g j J j − V g j ) + C J j J j + E J j cos ϕ J j (cid:21) + X j =1 (cid:20) C r j r j − Φ r j L r j (cid:21) + X j =1 (cid:20) C s s j + 2 E J s cos ( ϕ ( j ) ext ) cos ϕ s j (cid:21) + C c J − ˙Φ s ) + C c J − ˙Φ s ) + C c r − ˙Φ s ) + C c r − ˙Φ s ) + C c J − ˙Φ s ) + C c J − ˙Φ s ) + C c r − ˙Φ s ) + C c r − ˙Φ s ) , (S19)where the Josephson phase is ϕ j = 2 π Φ j / Φ , with Φ = h/ e the superconducting flux quantum, and 2 e is theelectrical charge of a Cooper pair.By applying the Legendre transformation H (Φ j , Q j ) = P j Q j ˙Φ j − L and making some simplifications, we obtainthe Hamiltonian H = 1¯ C J ( Q J − e ¯ n g ) − E J cos ( ϕ J ) + γ J (sin ϕ J ) + 1¯ C J ( Q J − e ¯ n g ) − E J cos ( ϕ J ) + γ J (sin ϕ J ) + 12 ¯ C r Q r + Φ r L r + γ r Φ r + 12 ¯ C r Q r + Φ r L r + γ r Φ r − γ J r sin ϕ J Φ r − γ J r sin ϕ J Φ r − γ J r sin ϕ J Φ r − γ J r sin ϕ J Φ r , (S20)where Q j is the conjugate momenta (node charge) for each node,¯ C J j = C j + 2 C c , ¯ C r j = C r j + 2 C c , ¯ n g j = − C g j V g j e , (S21)and γ J = (cid:20) E J s cos ( ϕ (1) ext ) + 1 E J s cos ( ϕ (2) ext ) (cid:21) C c E J C + 2 C c ) , γ J = (cid:20) E J s cos ( ϕ (3) ext ) + 1 E J s cos ( ϕ (4) ext ) (cid:21) C c E J C + 2 C c ) ,γ r = (cid:20) E J s cos ( ϕ (1) ext ) + 1 E J s cos ( ϕ (3) ext ) (cid:21) C c Φ π ( C r + 2 C c ) L r ,γ r = (cid:20) E J s cos ( ϕ (2) ext ) + 1 E J s cos ( ϕ (4) ext ) (cid:21) Φ C c π ( C r + 2 C c ) L r ,γ J r = Φ C c E J π ( C r + 2 C c )( C + 2 C c ) L r E J s cos ( ϕ (1) ext ) , γ J r = Φ C c E J π ( C r + 2 C c )( C + 2 C c ) L r E J s cos ϕ (2) ext ,γ J r = Φ C c E J π ( C r + 2 C c )( C + 2 C c ) E J s cos ϕ (3) ext L r , γ J r = Φ C c E J π ( C r + 2 C c )( C + 2 C c ) E J s cos ϕ (4) ext L r . (S22)Then, we quantize this Hamiltonian by promoting the classical variables to quantum operators, Q j → ˆ Q j = 2 e ˆ n j and ϕ j → ˆ ϕ j with the commutation relation [ e i ˆ ϕ j , ˆ n j ] = e i ˆ ϕ j , and Q r j → ˆ Q r j = q ~ ω r j ¯ C r j / a † + a ), and Φ r j → ˆΦ r j = i q ~ ω r j ¯ L r j / a − a † ) with the commutation [ a, a † ] = 1, where ¯ L r j = /L r +2 γ rj , and ω r j = q / ¯ C r j ¯ L r + j .Thus, the quantum Hamiltonian readsˆ H = X j =1 ˆ H jq + X j =1 ˆ H jr − γ J r sin ϕ J Φ r − γ J r sin ϕ J Φ r − γ J r sin ϕ J Φ r − γ J r sin ϕ J Φ r . (S23)Here, ˆ H jq = 4 E C j (ˆ n j − ¯ n g j ) − E J j cos ( ˆ ϕ j ) + γ J j ( ϕ ext ) sin ( ˆ ϕ j ) (S24)0is the Hamiltonian of the j th qubit with charge energy E C j = e / C J i , and ˆ H jr = ~ ω r j a † j a j is the Hamiltonian of the j th resonator. Moreover, in the charge number basis, it means ˆ n j = P m m | m j ih m j | , where | m j i is the m th excitedstate of the j th subsystem, cos ( ˆ ϕ j ) and sin ( ˆ ϕ j ) readcos ( ˆ ϕ j ) = 12 X m | m j ih m j + 1 | + H . c . ! , sin ( ˆ ϕ j ) = − i X m | m j ih m j + 1 | − H . c . ! . (S25)In the following discussion, we consider ¯ n g = ¯ n g = 0 . ~ = 1. Note that the free Hamiltonian of the subsystem H jq in Eq. (S24) includes both the pure CPB free Hamiltonian and the nonlinear term proportional to sin ( ˆ ϕ J j ) term.However, as long as we keep E J j /E C j in the charge regime, the increase of γ J j will not destroy the anharmonicity ofour system. The level of anharmonicity still depends on the ratio E J j /E C j . Thus, in the charge regime, we can safelyperform the two-level approximation. Therefore, the operator sin ˆ ϕ J j in the subsystem basis now readssin ˆ ϕ J j = 12 σ yj , (S26)where σ αj is the Pauli matrix and I j is the identity operator. Accordingly, the nonlinear term γ J j sin ( ˆ ϕ J j ) , asdisplayed in Eq. (S24), can be approximated as γ J j sin ( ˆ ϕ J j ) ≈ γ J j I , (S27)which only provides a shift to the qubit frequency. Finally, we obtain the simplified Hamiltonian as followsˆ H = ω q σ z + ω q σ z + ω r a † a + ω r a † a − i ˜ γ J r σ y ( a − a † ) − i ˜ γ J r σ y ( a − a † ) − i ˜ γ J r σ y ( a − a † ) − i ˜ γ J r σ y ( a − a † ) , (S28)where ω q = E J and ω q = E J , and the effective coupling strength˜ γ J r = Φ C c E J π ( C r + 2 C c )( C + 2 C c ) L r E J s cos ( ϕ (1) ext ) r ~ ω r ¯ L r , ˜ γ J r = Φ C c E J π ( C r + 2 C c )( C + 2 C c ) L r E J s cos ϕ (2) ext r ~ ω r ¯ L r , ˜ γ J r = Φ C c E J π ( C r + 2 C c )( C + 2 C c ) E J s cos ϕ (3) ext L r r ~ ω r ¯ L r , ˜ γ J r = Φ C c E J π ( C r + 2 C c )( C + 2 C c ) E J s cos ϕ (4) ext L r r ~ ω r ¯ L r . (S29)Now, we consider the external flux ϕ ( j ) ext to be composed by a DC signal and a small AC signal as ϕ ( j ) ext = ϕ ( j ) ext ( t ) = ϕ ( j ) DC + ϕ ( j ) AC ( t ), where ϕ ( j ) AC ( t ) = A ( j )1 cos ( ν ( j )1 t + ˜ ϕ ( j )1 ) + A ( j )2 cos ( ν ( j )2 t + ˜ ϕ ( j )2 ) , (S30)with | A ( j )1 | , | A ( j )2 | ≪ | ϕ DC | . Then, we can approximate1 E J s cos ( ϕ ( j ) ext ) ≈ E J s " ϕ ( j ) DC )cos ( ϕ ( j ) DC ) ϕ ( j ) AC ( t ) , (S31)where ¯ E ( j ) J s = E J s cos ϕ ( j ) DC . By replacing Eq. (S31) in the Hamiltonian of Eq. (S23), we obtainˆ H = ω q σ z + ω q σ z + ω r a † a + ω r a † a − i h g (1)0 + g (1)1 ϕ (1) AC ( t ) i σ y ( a − a † ) − i h g (2)0 + g (2)1 ϕ (2) AC ( t ) i σ y ( a − a † ) − i h g (3)0 + g (3)1 ϕ (3) AC ( t ) i σ y ( a − a † ) − i h g (4)0 + g (4)1 ϕ (4) AC ( t ) i σ y ( a − a † ) , (S32)where the coupling strength g ( j )0 = Φ C c E J q ~ ω rj ¯ L rj π ( C r j + 2 C c )( C + 2 C c ) L r j ¯ E ( j ) J s , g ( j )1 = Φ C c E J q ~ ω rj ¯ L rj π ( C r j + 2 C c )( C + 2 C c ) L r j ¯ E ( j ) J s sin ( ϕ ( j ) DC )cos ( ϕ ( j ) DC ) j = { , } , (S33)1and g ( j )0 = Φ C c E J q ~ ω rj ¯ L rj π ( C r j + 2 C c )( C + 2 C c ) L r j ¯ E ( j ) J s , g ( j )1 = Φ C c E J q ~ ω rj ¯ L rj π ( C r j + 2 C c )( C + 2 C c ) L r j ¯ E ( j ) J s sin ( ϕ ( j ) DC )cos ( ϕ ( j ) DC ) j = { , } . (S34)To visualize the dynamics of our system, we go to the interaction picture and make the rotating-wave approximationto obtainˆ H I ≈ Ag (1)1 (cid:20) − σ x a ( e − i ˜ ϕ (1)1 − e i ˜ ϕ (1)2 ) − σ y a ( ie − i ˜ ϕ (1)1 + ie i ˜ ϕ (1)2 ) + σ x a † ( e − i ˜ ϕ (1)2 − e i ˜ ϕ (1)1 ) + σ y a † ( ie − i ˜ ϕ (1)2 + ie i ˜ ϕ (1)1 ) (cid:21) + Ag (2)1 (cid:20) − σ x a ( e − i ˜ ϕ (2)1 − e i ˜ ϕ (2)2 ) − σ y a ( ie − i ˜ ϕ (2)1 + ie i ˜ ϕ (2)2 ) + σ x a † ( e − i ˜ ϕ (2)2 − e i ˜ ϕ (2)1 ) + σ y a † ( ie − i ˜ ϕ (2)2 + ie i ˜ ϕ (2)1 ) (cid:21) + Ag (3)1 (cid:20) − σ x a ( e − i ˜ ϕ (3)1 − e i ˜ ϕ (3)2 ) − σ y a ( ie − i ˜ ϕ (3)1 + ie i ˜ ϕ (3)2 ) + σ x a † ( e − i ˜ ϕ (3)2 − e i ˜ ϕ (3)1 ) + σ y a † ( ie − i ˜ ϕ (3)2 + ie i ˜ ϕ (3)1 ) (cid:21) + Ag (4)1 (cid:20) − σ x a ( e − i ˜ ϕ (4)1 − e i ˜ ϕ (4)2 ) − σ y a ( ie − i ˜ ϕ (4)1 + ie i ˜ ϕ (4)2 ) + σ x a † ( e − i ˜ ϕ (4)2 − e i ˜ ϕ (4)1 ) + σ y a † ( ie − i ˜ ϕ (4)2 + ie i ˜ ϕ (4)1 ) (cid:21) . By considering ˜ ϕ ( j )1 = π , and ˜ ϕ ( j )2 = 2 π , we obtainˆ H I ≈ Ag (1)1 σ x ( a + a † ) + Ag (2)1 σ x ( a + a † ) + Ag (3)1 σ x ( a + a † ) + Ag (4)1 σ x ( a + a † ) , which corresponds to the two-qubit two-mode quantum Rabi model. We highlight that the coupling strengths Ag ( k )1 / A , such that they can be adiabatically changed in order to perform the W-stategeneration protocol proposed in the main text. S4. LINDBLAD MASTER EQUATION FOR NUMERICAL SIMULATION
We use the following Lindblad form master equation˙ ρ = − i [ H pq , ρ ] + M X i =1 κ a i ρa † i − a † i a i ρ − ρa † i a i ) + X j =1 γ j σ j ρσ † j − σ † j σ j ρ − ρσ † j σ j )+ X j =1 γ jφ ( σ jz ρσ jz − ρ ) (S35)to carry out the numerical simulation. Here, κ is the photon decay rate of the i th CWR, consisting of the intrinsicpart κ in and coupling part κ c with respect to the TL. γ j and γ jφ are the energy relaxation rate and the dephasingrate of the j th SQ, respectively. Although the ultrastrong coupling regime is reached, we use this Lindblad formmaster equation because the maximum coupling strength g we reach is just 0 .
29 and κ in = 10 − ω , γ j = 10 − ω , γ jφ = 10 − ω , which are extremely small and when κ c = 10 − ω is turned on, the qubit Bell state has been generated,being decoupled from the photon mode. We have testified it numerically by using a Markovian master equation [1, 2]˙ ρ ( t ) = − i [ H pq , ρ ( t )] + X j,k>j Γ jkm D ( | j ih k | ) ρ ( t ) , which corresponds to amplitude damping. Here, {| j i} j =0 , , .. are the eigenvectors of Hamiltonian H pq , with H | j i = ǫ j | j i , and D ( O ) ρ = 1 / OρO † − O † Oρ − ρO † O ). Also, index m = { , ..., M, M + 1 , ..., M + N } runs through allresonators and qubits, such that m ≤ M refers to photon modes and M < m ≤ M + N refers to qubits. The decayrates, Γ jkm , are taken as Γ jkm = ( κ m δ kj ω | C mkj | . C m = a m + a † m for m ≤ Mγ m δ kj ω | C mkj | . C m = σ mx for M < m ≤ M + N kj = ǫ k − ǫ j and C mkj = h k | C m | j i . κ m is the photon decay rate of the m -th CWR, while γ m is the energyrelaxation rate of the m -th qubit. The numerical results is almost the same as that obtained form the correspondingLindblad form master equation Eq. (S35). [1] F. Beaudoin, J. M. Gambetta, and A. Blais, Phys. Rev. A , 043832 (2011).[2] A. Ridolfo, M. Leib, S. Savasta, and M. J. Hartmann, Phys. Rev. Lett.109