Harvesting Multiqubit Entanglement from Ultrastrong Interactions in Circuit Quantum Electrodynamics
Federico Armata Giuseppe Calajo Tuomas Jaako M. S. Kim Peter Rabl
quant-ph cond-mat.supr-con


















^1QOLS and QuEST, Blackett Laboratory, Imperial College London, London SW7 2AZ, United Kingdom 
 
^2Vienna Center for Quantum Science and Technology, Atominstitut, TU Wien, 1040 Vienna, Austria

^3Korea Institute of Advanced Study, Dongdaemun-gu, Seoul, 02455, South Korea


We analyze a multiqubit circuit QED system in the regime where the qubit-photon coupling dominates over the system's bare energy scales. Under such conditions a manifold of low-energy states with a high degree of entanglement emerges. Here we describe a time-dependent protocol for extracting these quantum correlations and converting them into well-defined multipartite entangled states of noninteracting qubits. Based on a combination of various  ultrastrong-coupling effects the protocol can be operated in a fast and robust manner, while still being consistent with experimental constraints on switching times and typical energy scales encountered in superconducting circuits. Therefore, our scheme can serve as a probe for otherwise inaccessible correlations in strongly coupled circuit QED systems. It also shows how such correlations can potentially be exploited as a resource for entanglement-based applications. 





Harvesting Multiqubit Entanglement from Ultrastrong Interactions in Circuit Quantum Electrodynamics
    P. Rabl^2
    December 30, 2023
===================================================================================================



Cavity QED is the study of quantum light-matter interactions with real or artificial two-level atoms  coupled to a single radiation mode. In this context one is usually interested in strong interactions between excited atomic and electromagnetic states, while  the trivial ground state, i.e. the vacuum state with no atomic or photonic excitations, plays no essential role. This paradigm has recently been challenged by a number of experiments <cit.>, where interaction strengths comparable to the photon energy have been demonstrated.  In particular, in the field of circuit QED <cit.>, a single superconducting two-level system can already be coupled ultrastrongly <cit.> to a microwave resonator mode <cit.>. In this regime the physics changes drastically and even in the  ground state various nontrivial effects like spontaneous vacuum polarization <cit.>, light-matter decoupling <cit.> and different degrees of entanglement <cit.> can occur. However, compared to the vast literature on cavity QED systems in the weakly coupled regime, the opposite limit of extremely strong interactions is to a large extent still unexplored. As a consequence, ideas for how ultrastrong coupling (USC) effects can be controlled and exploited for practical applications are limited <cit.>.
 
 
In this Letter we consider a prototype circuit QED system consisting of multiple flux qubits coupled to a single mode of a microwave resonator. It has recently been shown that in the USC regime this circuit exhibits a manifold of nonsuperradiant ground and low-energy  states with a high degree of multiqubit entanglement <cit.>. This entanglement, however, is a priori not of any particular use, since any attempt to locally manipulate or measure the individual qubits would necessarily introduce a severe perturbation to the strongly coupled system. For this reason we describe  the implementation of an entanglement-harvesting  protocol <cit.>, which extracts quantum correlations from USC states and converts these correlations into equivalent multipartite entangled states of decoupled qubits. The protocol combines adiabatic and nonadiabatic parameter variations and exploits the counterintuitive decoupling of qubits and photons at very strong interactions <cit.> to make the entanglement extraction scheme intrinsically robust and consistent with experimentally available tuning capabilities. The extracted Dicke and singlet states belong to a family of robust multipartite entangled states <cit.> and form, for example,  a resource for Heisenberg-limited metrology applications <cit.>.  More generally, our analysis shows, how the interplay between different USC effects can contribute to the realization of non-trivial control tasks in a strongly interacting cavity QED system.

 



Model.—We consider a circuit QED system as shown in Fig. <ref>(a), where a single mode LC resonator with capacitance C and inductance L is coupled collectively to an even number of N=2,4,6,… flux qubits. This circuit is described by the Hamiltonian <cit.>

    ℋ=Q_r^2/2C+(Φ_r- Φ_0 ∑_i=1^N  φ_i)^2/2L+∑_i=1^N H_q^(i),

where Q_r and Φ_r are charge and generalized flux operators for the resonator obeying [Φ_r,Q_r]=iħ, and Φ_0=ħ/(2e) is the reduced flux quantum.  
For each qubit, H_q^(i) denotes the free Hamiltonian and φ_i is the difference of the superconducting phase across the qubit's subcircuit. As usual we assume that the qubit dynamics can be restricted to the two lowest tunneling states |↓⟩ and |↑⟩ of a symmetric double-well potential [cf. Fig. <ref>(b)].  Under this approximation and writing Φ_r=  √(ħ/(2C ω_r)) (a+a^†) and Q_r=i √(ħ C ω_r/2) (a^†-a), where ω_r=√(1/LC) is the resonator frequency and a and a^† are the annihilation and creation operators, we obtain 

    ℋ=   ħω_r a
    + ħ∑_i=1^N g_i/2(+a)σ^i_x
        +ħ∑_i=1^Nω_q^i/2σ_z^i + ħ∑_i,j=1^N g_ig_j/4ω_rσ^i_xσ_x^j.

Here σ_k^i are Pauli operators and ω^i_q are the qubit-level splittings. The second term in Eq. (<ref>) accounts for the collective qubit-resonator interaction with couplings g_i=Φ_0√(|φ_0^i|^2ω_r/(2ħ L)), where φ_0^i=2⟨↓_i|φ_i |↑_i⟩. The condition g_i>ω_r,ω^i_q can be reached with an appropriate flux-qubit design <cit.>, and the g_i(t) and ω_q^i(t) can be individually tuned by controlling the matrix element φ_0^i and the height of the tunnel barrier via local magnetic fluxes <cit.>. A specific four-junction qubit design <cit.>, which combines strong interactions with a high degree of tunability, is detailed in the Supplemental Material <cit.>. Finally, the last contribution in Eq. (<ref>)  represents an additional qubit-qubit interaction, which is usually neglected for cavity QED systems with weak or moderately strong couplings. However, this term is crucial in the USC regime and it is responsible for the nontrivial ground-state correlations that  are at the focus of the present Letter. 






USC spectrum.—We are primarily interested in a  symmetric configuration, i.e., g_i=g and ω_q^i=ω_q. In this case the Hamiltonian (<ref>) can be expressed in terms of collective angular momentum operators S_k=∑_i σ^i_k/2 and reduces to the extended Dicke Hamiltonian <cit.>

    ℋ=ħω_r a + ħ g(+a)S_x +ħω_qS_z + ħ g^2/ω_rS_x^2.
 
For g≪ω_r,ω_q we can make a rotating wave approximation and obtain the standard Tavis-Cummings model of cavity QED with a trivial ground state |G⟩=|n=0⟩⊗|↓⟩^⊗ N. If in addition |ω_q-ω_r|≫ g, all excited states are also essentially decoupled and the qubits can be individually prepared, manipulated, and measured by additional control fields. In the opposite limit, g≫ω_r,ω_q, the coupling terms ∼ S_x and ∼ S_x^2 dominate and the level structure  changes completely. This is illustrated in Fig. <ref>(c), which shows that for couplings g/ω_r≳ 3 the spectrum separates into manifolds of 2^N nearly degenerate states. The eigenstates in this regime are  displaced photon number states, |Ψ_s,m_x,n⟩≃ e^-g/ω_r(a^†-a)S_x|n⟩⊗|s,m_x⟩, with energies E_s,m_x,n≃ħω_r n + δ E^(n)_s,m_x <cit.>. Here  s is the total spin and m_x=-s,…,s the spin projection quantum number; i.e., S_x|s,m_x⟩=m_x |s,m_x⟩. Within the lowest manifold, the remaining level splittings are given by

    δ E^(0)_s,m_x=  ħΔ[m_x^2- s(s+1)],    Δ= ω_q^2 ω_r/2g^2,

and the resulting ordering of the states is shown in Fig. <ref>(d) for N=4 qubits. Thus, for even qubit numbers N, the ground state in the USC regime is of the form |G̃⟩≃ |n=0⟩⊗ |D_0⟩, where |D_0⟩ = |s=N/2,m_x=0⟩ denotes the fully symmetric Dicke state with vanishing projection along x. Importantly, this state exhibits a high degree of qubit-qubit entanglement, while it remains almost completely decoupled from the cavity field <cit.>. Our goal is now to identify a suitable protocol for converting this state into an equivalent state of the decoupled system, where it becomes available as an entanglement resource for further use.







Entanglement harvesting.—Figure <ref>(a) shows a general pulse sequence for implementing the  entanglement-harvesting protocol through variations of  ω_q(t) and g(t). For this protocol the system is initialized in the ground state |G⟩ of the weakly coupled system, where the qubits are far detuned from the cavity, ω_q=ω_ max≫ω_r, and the coupling is set to a minimal value, g=g_ min<ω_r. In the first two steps, T_1 and T_2, the system is adiabatically tuned into the USC regime with a maximal coupling g_ max>ω_r and a low value of the qubit frequency ω_ min≲ω_r. This process prepares the system in the USC ground state |G̃⟩. In the successive steps, T_3 and T_4, the qubits and the resonator mode are separated again, but now in the reverse order and using nonadiabatic parameter variations.  Ideally, during this part of the protocol the system simply remains in state |G̃⟩ and becomes the desired excited state of the weakly coupled system at the final time T_ f=∑_n=1^4 T_n. This general sequence achieves two main goals. First, the adiabatic preparation stage can be implemented very rapidly, since it must only be slow compared to the fast time scales set by ω_ max^-1 and g_ max^-1. At the same time the nonadiabatic decoupling processes only need to be fast compared to the slow time scales ω^-1_r, ω^-1_ min, and g_ min^-1. This second condition is most crucial for a time-dependent control of USC systems, since it makes the required switching times experimentally accessible and consistent with the two-level approximation assumed in our theoretical model. 



In Fig. <ref>(b) we plot the fidelity ℱ(t)={ρ(t) |D_0⟩⟨ D_0|}, where ρ(t) is the density operator of the full system, for a specific set of pulse parameters listed in the figure caption. We see that the entanglement extraction fidelity (EEF) ℱ_ E = max{ℱ(t)| t ≥ T_ f}, i.e., the maximal fidelity after the decoupling step, reaches near perfect values of ℱ_ E≃ 0.95-0.99 for different numbers of qubits, without any further fine-tuning of the control pulses.  Note that the fidelity oscillations at the end of the sequence are simply due to the fact that |D_0⟩ is not an eigenstate of the bare qubit Hamiltonian, H_q=ω_q S_z. However, this evolution does not affect the purity or the degree of entanglement of the final qubit state and can be undone by local qubit rotations. 









Experimental considerations.—For a possible experimental implementation of the protocol we consider qubits with a frequency of ω_ max/(2π)≈10 GHz coupled to a lumped-element  resonator of frequency ω_r/(2π)=500 MHz. The required maximal coupling strength of g_ max≃ 4.5ω_r ≈ 2π× 2.25  GHz is then consistent with experimentally demonstrated values <cit.>. For these parameters, the nonadiabatic switching times assumed in Fig. <ref>(b) correspond to T_3,4≃ 0.16 ns. These switching times are within reach of state-of-the-art waveform generators and a sinusoidal modulation of flux qubits on such time scales has already been demonstrated <cit.>. At the same time the duration of the whole protocol, T_ f=15/ω_r≈ 5 ns, is still much faster than typical flux qubit coherence times of  1-100 μs <cit.> or the lifetime of a photon, T_ ph=Q/ω_r, in a microwave resonator of quality factor Q=10^4-10^6.  Therefore, although many experimental techniques for implementing and operating circuit QED systems in the USC regime are still under development, these estimates clearly demonstrate the feasibility of realizing high-fidelity control operations in such devices.  


In practice additional limitations might arise from the lack of complete tunability of g(t) and ω_q(t). This is illustrated in Fig. <ref>(a), which shows the evolution of the lowest eigenenergies during different stages of the protocol for the case N=2 and a nonvanishing value of g_ min. In this case the appearance of several avoided crossings during the final ramp-up step prevents a fully nonadiabatic decoupling. In Fig. <ref>(b) we plot the resulting EEF for varying g_ min and T_4. This plot demonstrates the expected trade-off between the residual coupling and the minimal switching time, but also that the protocol is rather robust and fidelities of ℱ_ E∼ 0.9 are still possible for minimal  couplings of a few hundred MHz or switching times approaching ∼ 1 ns. 
Similar conclusions are obtained when a partial dependence between the pulses for g(t) and ω_q(t) or nonuniform couplings g_i(t) and frequencies ω_q^i(t) due to fabrication uncertainties are taken into account. Numerical simulations of the protocol under such realistic experimental conditions <cit.> demonstrate that no precise fine-tuning of the system parameters is required. 




Extracting entanglement from a thermal state.—The above-considered  protocol relies on a rather low resonator frequency ω_r in order to enhance both g/ω_r as well as  the nonadiabatic switching times. This implies that even at temperatures of T=20 mK the equilibrium populations of higher resonator states with n≥ 1 cannot be neglected.  In Fig. <ref>(c) we plot the EEF as a function of the  temperature T, assuming an initial resonator state ρ_ th=∑_n p_n|n⟩⟨n|, where p_n=n̅^n/(1+n̅)^n+1 is the thermal distribution for a mean excitation number n̅=1/(e^ħω_r/k_BT-1). We see that the EEF is significantly higher than one would naively expect from the initial population in the ground state |G⟩. The origin of this surprising effect can  be understood from the eigenvalue plot in Fig. <ref>(a). For example, the weak-coupling eigenstate |n=1⟩⊗|↓⟩^⊗ N is efficiently mapped on the corresponding USC state |n=1⟩⊗|s=N/2,m_x=0⟩, passing only through a weak, higher-order avoided crossing. Therefore, the intermediate—and as a result also the final—qubit state is one with the resonator being in state |1⟩. Although for higher photon numbers the avoided crossings become more relevant, the protocol still approximately implements the mapping |n⟩⊗|↓⟩^⊗ N→ |n⟩⊗|s=N/2,m_x=0⟩, independent of the resonator state |n⟩.  This feature makes it rather insensitive to thermal occupations and avoids additional active cooling methods for initializing the system in state |G⟩.




Entanglement protection.—Figure <ref>(d) shows that apart from the ground state |G̃⟩ there are many other highly entangled states within the lowest USC manifold. Of particular interest is the energetically highest state  |Ẽ⟩=|n=0⟩⊗|S⟩, where |S⟩ is a singlet state with total angular momentum s=0 and S_z|S⟩=S_x|S⟩=0. Therefore, once prepared, this state is an exact dark state of Hamiltonian (<ref>) and remains decoupled from the cavity field in all parameter regimes. Although this state is not connected to any of the bare qubit states in a simple adiabatic way, it can still be harvested by an adopted protocol, as described in Fig. <ref>(a) for the case N=4. For this protocol the system is initially prepared in the excited state |Ψ_0⟩=|0⟩⊗|↑↑↓↓⟩ and in a first step the qubit states are lowered below the resonator frequency in order to avoid further level crossings with higher-photon-number states. The increase of the coupling combined with a frequency offset to break the angular momentum conservation then evolves the system into a state with s=0 already for moderate couplings of g/ω_r≈ 1.8. Note that for N≥ 4 there are  multiple degenerate USC states with s=0 <cit.> [cf. Fig. <ref>(d)], out of which the protocol selects a specific superposition <cit.>.

Although the harvesting protocol for state |S⟩ loses some of the robustness of the ground-state protocol, it adds an important feature. By retaining a finite coupling g_ f=g(t=T_ f)∼ω_r at the end of the protocol, the extracted dark state |S⟩ is energetically separated from all other states with s≠0 and it is thereby protected against small frequency fluctuations. This effect is illustrated in Fig. <ref>(c), which shows the evolution of the extracted state |S⟩ in the presence of small random shifts of the individual qubit frequencies. For g_ f=0 this leads to dephasing of the qubits and a rapid transition out of the s=0 subspace. This dephasing can be substantially suppressed by keeping the coupling at a finite value. Thus, this example shows that USC effects can be used not only to generate complex multiqubit entangled states, but also to protect them. 

Conclusion.—We have presented a  protocol for extracting well-defined multiqubit-entangled states from the ground-state manifold of an ultrastrongly coupled circuit QED system. The detailed analysis of this protocol illustrates, how various–so far unexplored–USC effects can contribute to a robust generation and protection of complex multiqubit states. These principles can serve as a guideline for many other preparation, storage and control operations in upcoming USC circuit QED experiments with two or more qubits. 

Acknowledgments.—This work was supported by the Austrian Science Fund (FWF) through the SFB FoQuS, Grant No. F40, the DK CoQuS, Grant No. W 1210, and the START Grant No. Y 591-N16, and the People Programme (Marie Curie Actions) of the European Union's Seventh Framework Programme (FP7/2007-2013) under REA Grant No. 317232. F. A. wishes to express his gratitude to the Quantum Optics Theory Group of Atominstitut (TU Wien) for the warm hospitality received during his numerous visits. M. S. K. acknowledges support from the UK EPSRC Grant No. EP/K034480/1 and the Royal Society.



99


Todorov2010 Y. Todorov, A. M. Andrews, R. Colombelli, S. De Liberato, C. Ciuti, P. Klang, G.
Strasser, and C. Sirtori, “Ultrastrong Light-Matter Coupling Regime with Polariton Dots",
Phys. Rev. Lett. 105, 196402 (2010).

schwartz11 T. Schwartz, J. A. Hutchison, C. Genet, and T. W. Ebbesen, “Reversible Switching of Ultrastrong Light-Molecule Coupling", Phys. Rev. Lett. 106, 196405 (2011).

geiser12 M. Geiser, F. Castellano, G. Scalari, M. Beck, L. Nevou, and J. Faist, “Ultrastrong Coupling Regime and Plasmon Polaritons in Parabolic Semiconductor Quantum Wells", Phys. Rev. Lett. 108, 106402 (2012).


Scalari2015 G. Scalari, C. Maissen, S. Cibella, R. Leoni, C. Reichl, W. Wegscheider, M. Beck, and J. Faist, “THz ultrastrong light-matter coupling", Il Nuovo Saggiatore 31, 4 (2015). 

Zhang2016 Q. Zhang, M. Lou, X. Li, J. L. Reno, W. Pan, J. D. Watson, M. J. Manfra, and J. Kono,
“Collective non-perturbative coupling of 2D electrons with high-quality-factor terahertz
cavity photons", Nat. Phys. 12, 1005 (2016).




wallraff04 A. Wallraff, D. I. Schuster, A. Blais, L. Frunzio, R. S. Huang, J. Majer, S. Kumar, S. M. Girvin, and R. J. Schoelkopf, “Strong coupling of a single photon to a superconducting qubit using circuit quantum electrodynamics", Nature (London) 431, 162 (2004).

blais04 A. Blais, R.-S. Huang, A. Wallraff, S. M. Girvin and R. J. Schoelkopf, “Cavity quantum electrodynamics for superconducting electrical circuits: An architecture for quantum computation", Phys. Rev. A 69, 062320 (2004).



Ciuti2005 C.  Ciuti,  G.  Bastard,  and  I.  Carusotto, “Quantum vacuum properties of the intersubband cavity polariton field", Phys.  Rev.  B 72,  115303  (2005). 

Casanova2010 J. Casanova, G. Romero, I. Lizuain, J. J. García-Ripoll, and E. Solano, “Deep Strong Coupling Regime of the Jaynes-Cummings Model", Phys. Rev. Lett. 105, 263603 (2010).

USCRemark Note that in this Letter we do not further distinguish between the USC <cit.> and the deep strong-coupling regime <cit.>.   



niemczyk10 T. Niemczyk, et al., “Circuit quantum electrodynamics in the ultrastrong-coupling regime", Nat. Phys. 6, 772 (2010).

forndiaz10 P. Forn-Diaz, J. Lisenfeld, D. Marcos, J. J. Garcia-Ripoll, E. Solano, C. J. P. M. Harmans, and J. E. Mooij, “Observation of the Bloch-Siegert Shift in a Qubit-Oscillator System in the Ultrastrong Coupling Regime", Phys. Rev. Lett. 105, 237001 (2010).

baust16 A. Baust, E. Hoffmann, M. Haeberlein, M. J. Schwarz, P. Eder, J. Goetz, F. Wulschner, E. Xie, L. Zhong, F. Quijandría, D. Zueco, J.-J. García-Ripoll, L. García-Álvarez, G. Romero, E. Solano, K. G. Fedorov, E. P. Menzel, F. Deppe, A. Marx, and R. Gross, “Ultrastrong coupling in two-resonator circuit QED", Phys. Rev. B 93, 214501 (2016).

forndiaz17 P. Forn-Diaz, J. J. García-Ripoll, B. Peropadre, M. A. Yurtalan, J.-L. Orgiazzi, R. Belyansky, C. M. Wilson, and  A. Lupascu, “Ultrastrong coupling of a single artificial atom to an electromagnetic continuum", Nat. Phys. 13, 39 (2017).

yoshihara17 F. Yoshihara, T. Fuse, S. Ashhab, K. Kakuyanagi, S. Saito, and K. Semba, “Superconducting qubit-oscillator circuit beyond the ultrastrong-coupling regime", Nat. Phys. 13, 44 (2017).


chen17 Z. Chen, Y. Wang, T. Li, L. Tian, Y. Qiu, K. Inomata, F. Yoshihara, S. Han, F. Nori, J. S. Tsai, and J. Q. You, “Multi-photon sideband transitions in an ultrastrongly coupled circuit quantum electrodynamics system", Phys. Rev. A 96, 012325 (2017).

bosman17 S. J. Bosman, M. F. Gely, V. Singh, A. Bruno, D. Bothner, and G. A. Steele, “Multi-mode ultra-strong coupling in circuit quantum electrodynamics", arXiv:1704.06208 (2017).


nataff10a P. Nataf and C. Ciuti, “Vacuum Degeneracy of a Circuit QED System in the Ultrastrong Coupling Regime", Phys. Rev. Lett. 104, 023601 (2010).

nataff10b P. Nataf and C. Ciuti, “No-go theorem for superradiant quantum phase transitions in cavity QED and counter-example in circuit QED", Nat. Commun. 1, 72 (2010).

Bamba2016 M.  Bamba, K. Inomata, and Y. Nakamura, “Superradiant Phase Transition in a Superconducting Circuit in Thermal Equilibrium", Phys. Rev. Lett. 117, 173601 (2016).

deliberato14 S. De Liberato,“Light-Matter Decoupling in the Deep Strong Coupling Regime: The Breakdown of the Purcell Effect", Phys. Rev. Lett. 112, 016401 (2014).

jaako2016 T. Jaako, Z.-L. Xiang, J.J. Garcia-Ripoll, and P. Rabl, “Ultrastrong coupling phenomena beyond the Dicke model", Phys. Rev. A 94, 033850 (2016).

levine04 G. Levine and V. N. Muthukumar, “Entanglement of a qubit with a single oscillator mode", Phys. Rev. B 69, 113203 (2004).

hines04 A. P. Hines, C. M. Dawson, R. H. McKenzie, and G. J. Milburn, “Entanglement and bifurcations in Jahn-Teller models", Phys. Rev. A 70, 022303 (2004).

ashhab10 S. Ashhab and F. Nori, “Qubit-oscillator systems in the ultrastrong-coupling regime and their potential for preparing nonclassical states", Phys. Rev. A 81, 042311 (2010).

nataff11 P. Nataf and C. Ciuti, “Protected Quantum Computation with Multiple Resonators in Ultrastrong Coupling Circuit QED", Phys. Rev. Lett. 107, 190402 (2011).

romero12 G. Romero, D. Ballester, Y. M. Wang, V. Scarani, and E. Solano, “Ultrafast Quantum Gates in Circuit QED", Phys. Rev. Lett. 108, 120501 (2012).

kyaw15 T. H. Kyaw, S. Felicetti, G. Romero, E. Solano, and L.-C. Kwek, “Scalable quantum memory in the ultrastrong coupling regime", Sci. Rep. 5, 8621 (2015).

wang16 Y. Wang, J. Zhang, C. Wu, J. Q. You, and G. Romero, “Holonomic quantum computation in the ultrastrong-coupling regime of circuit QED", Phys. Rev. A 94, 012328 (2016).

wang17 Y. Wang, C. Guo, G.-Q. Zhang, G. Wang, and C. Wu, “Ultrafast quantum computation in ultrastrongly coupled circuit QED systems", Sci. Rep. 7, 44251 (2017).

stassi17 R. Stassi and F. Nori, “Quantum Memory in the Ultrastrong-Coupling Regime via Parity Symmetry Breaking", arXiv:1703.08951 (2017).







Reznik2005 B. Reznik, A. Retzker, and J. Silman, “Violating Bell's inequalities in vacuum", Phys. Rev. A 71, 042104 (2005).

Han2008 M. Han, S. J. Olson, and J. P. Dowling, “Generating entangled photons from the vacuum by accelerated measurements: Quantum-information theory and the Unruh-Davies effect", Phys. Rev. A 78, 022302 (2008).

auer2011 A. Auer and G. Burkard, “Entangled photons from the polariton vacuum in a switchable optical cavity", Phys. Rev. B 85, 235140 (2012).

olson2012 S. J. Olson and T. C. Ralph, “Extraction of timelike entanglement from the quantum vacuum", Phys. Rev. A 85, 012306 (2012).

sabin2012 C. Sabin, B. Peropadre, M. del Rey, and E. Martin-Martinez, “Extracting Past-Future Vacuum Correlations Using Circuit QED", Phys. Rev. Lett. 109, 033602 (2012).

Salton2015 G. Salton, R. B. Mann, and N. C. Menicucci, “Acceleration-assisted entanglement harvesting and rangefinding", New J. Phys. 17, 035001 (2015).

dai2015 L. Dai, W. Kuo and M. C. Chung, “Extracting entangled qubits from Majorana fermions in quantum dot chains through the measurement of parity", Sci. Rep. 5, 11188 (2015).


Guhne2008 O. Gühne, F. Bodoky, and M. Blaauboer, “Multiparticle entanglement under the influence of decoherence", Phys. Rev. A
78, 060301(R) (2008).

Bergmann2013 M. Bergmann and O. Gühne, “Entanglement criteria for Dicke states",  	J. Phys. A: Math. Theor. 46, 385304 (2013).


Toth2014 G. Toth and I. Apellaniz, “Quantum metrology from a quantum information science perspective", J. Phys. A: Math. Theor. 47, 424006 (2014).



vool16 U. Vool and M. Devoret, Introduction to Quantum Electromagnetic Circuits, Int. J. Circ. Theor. Appl. 45, 897 (2017).

SuppInfo See Supplemental Material for additional details about the circuit design and the performance of the protocol under nonideal conditions, which includes Ref. <cit.>.

beaudoin11s-main F. Beaudoin, J. M. Gambetta, and A. Blais, “Dissipation and ultrastrong coupling in circuit QED", Phys. Rev. A 84, 043832 (2011).

ridolfo12s-main A. Ridolfo, M. Leib, S. Savasta, and M. J. Hartmann. “Photon Blockade in the Ultrastrong Coupling Regime", Phys. Rev. Lett. 109, 193602 (2012).

Bourassa2009 J.  Bourassa,  J.  M.  Gambetta,  A.  A.  Abdumalikov,  O. Astafiev,  Y.  Nakamura,  and  A.  Blais,  “Ultrastrong  coupling regime of cavity QED with phase-biased flux qubits", Phys. Rev. A 80, 032109 (2009).

Peropadre2013 B.  Peropadre,  D.  Zueco,  D.  Porras,  and  J.  J.  Garcia-Ripoll,  “Nonequilibrium  and  Nonperturbative  Dynamics of Ultrastrong Coupling in Open Lines", Phys. Rev. Lett. 111, 243602 (2013).


Qiu2016 Y. Qiu, W. Xiong, X.-L. He, T.-F. Li, and J. Q. You,  “Four-junction superconducting circuit", Sci. Rep. 6, 28622 (2016).






Wilson2011 C. M.  Wilson,  G.  Johansson,  A.  Pourkabirian,  J. R. Johansson, T. Duty, F. Nori, and P. Delsing, “Observation of the dynamical Casimir effect in a superconducting circuit", Nature (London) 479, 376 (2011).

Yan2016 F. Yan, S. Gustavsson, A. Kamal, J. Birenbaum, A.P. Sears, D. Hover, D. Rosenberg, G. Samach, T. J. Gudmundsen, J. L. Yoder, T. P. Orlando, J. Clarke, A. J. Kerman, and W. D. Oliver, “The flux qubit revisited to enhance coherence and reproducibility", Nat. Commun. 7, 12964 (2016).

Dicke1954 R. H. Dicke, “Coherence in spontaneous radiation processes", Phys. Rev. 93, 99 (1954).

Arecchi1972 F. T. Arecchi, E. Courtens, R. Gilmore, and H. Thomas, “Atomic coherent states in quantum optics", Phys. Rev. A 6, 2211 (1972).

















    Supplemental material






§ USC EIGENSTATES

In the USC regime, the eigenstates of the extended Dicke Hamiltonian are labeled by the total spin s=0,1,…,N/2 and the spin projection quantum number m_x=-s,…,s, i.e., S⃗^2|s,m_x⟩=s(s+1)|s,m_x⟩ and S_x|s,m_x⟩=m_x |s,m_x⟩. For N>2 the states |s≠ N/2,m_x≠± N/2⟩ appear as multiplets <cit.>, due to the permutation symmetry of the Hamiltonian.  In the following we provide an overview of the relevant spin states for the case of N=2 and N=4 qubits. These states are most conveniently expressed in the rotated basis 

    |↓⟩_x=1/√(2)(|↓⟩-|↑⟩),    |↑⟩_x=1/√(2)(|↓⟩+|↑⟩).






 §.§ 2 qubits

For two qubits we have the usual three triplet states

    |s=1,m_x=1⟩=|↑↑⟩_x,
       |s=1,m_x=0⟩=1/√(2)(|↑↓⟩_x+|↓↑⟩_x),
       |s=1,m_x=-1⟩=|↓↓⟩_x,

and the singlet 

    |s=0,m_x=0⟩=1/√(2)(|↑↓⟩_x-|↓↑⟩_x).

When expressed in terms of the original qubit basis the two m_x=0 states of interest read

    |s=1,m_x=0⟩=1/√(2)(|↓↓⟩-|↑↑⟩),    |s=0,m_x=0⟩=1/√(2)(|↑↓⟩-|↓↑⟩).






 §.§ 4 qubits

For the case of N=4 qubits we obtain  a quintuplet for s=2, three triplets for s=1 and a two states for s=0. The maximally symmetric states with s=2 are the usual Dicke states in the x-basis, i.e.,  

    |s=2,m_x=2⟩=|↑↑↑↑⟩_x,
       |s=2,m_x=1⟩=1/2(|↓↑↑↑⟩_x+|↑↓↑↑⟩_x+|↑↑↓↑⟩_x+|↑↑↑↓⟩_x),
       |s=2,m_x=0⟩=1/√(6)(|↓↓↑↑⟩_x+|↑↑↓↓⟩_x+|↓↑↑↓⟩_x+|↑↓↓↑⟩_x+|↓↑↓↑⟩_x+|↑↓↑↓⟩_x),
       |s=2,m_x=-1⟩=1/2(|↓↑↓↓⟩_x+|↓↓↓↑⟩_x+|↓↓↑↓⟩_x+|↑↓↓↓⟩_x),
       |s=2,m_x=-2⟩=|↓↓↓↓⟩_x.

For the entanglement harvesting protocol, we are interested in the state |s=2,m_x=0⟩, which in the original qubit basis is given by  

    |s=2,m_x=0⟩=3/√(24)(|↑↑↑↑⟩+|↓↓↓↓⟩)-1/√(24)(|↑↑↓↓⟩+|↑↓↑↓⟩+|↑↓↓↑⟩+|↓↑↑↓⟩+|↓↑↓↑⟩+|↓↓↑↑⟩).

Each of the s=1 states is 3-fold degenerate and the corresponding states are 

    |s=1,m_x=1⟩   =1/2(|↑↑↑↓⟩_x+|↑↑↓↑⟩_x-|↑↓↑↑⟩_x-|↓↑↑↑⟩_x)
    1/2(|↑↑↑↓⟩_x-|↑↑↓↑⟩_x+|↑↓↑↑⟩_x-|↓↑↑↑⟩_x)
    1/2(|↑↑↑↓⟩_x-|↑↑↓↑⟩_x-|↑↓↑↑⟩_x+|↓↑↑↑⟩_x)
    
    |s=1,m_x=0⟩   =1/√(2)(|↑↑↓↓⟩_x-|↓↓↑↑⟩_x)
    1/√(2)(|↑↓↑↓⟩_x-|↓↑↓↑⟩_x)
    1/√(2)(|↑↓↓↑⟩_x-|↓↑↑↓⟩_x)
    
    |s=1,m_x=-1⟩   =1/2(|↑↓↓↓⟩_x+|↓↑↓↓⟩_x-|↓↓↑↓⟩_x-|↓↓↓↑⟩_x)
    1/2(|↑↓↓↓⟩_x-|↓↑↓↓⟩_x+|↓↓↑↓⟩_x-|↓↓↓↑⟩_x)
    1/2(|↑↓↓↓⟩_x-|↓↑↓↓⟩_x-|↓↓↑↓⟩_x+|↓↓↓↑⟩_x)


Finally, for N=4 we obtain two singlet states with s=0, which are given by

    |s=0,m_x=0⟩   =
    |S⟩=  1/√(3)(|↑↑↓↓⟩_x+|↓↓↑↑⟩_x)-1/√(12)(|↑↓↑↓⟩_x+|↑↓↓↑⟩_x+|↓↑↑↓⟩+|↓↑↓↑⟩_x),
    
    |S^'⟩= 1/2(|↑↓↑↓⟩_x-|↑↓↓↑⟩_x-|↓↑↑↓⟩_x+|↓↑↓↑⟩_x) .

Since the s=0 subspace is invariant under rotations, the states have the same form as in the original qubit basis,

    |s=0,m_x=0⟩   =
    |S⟩= 1/√(3)(|↑↑↓↓⟩+|↓↓↑↑⟩)-1/√(12)(|↑↓↑↓⟩+|↑↓↓↑⟩+|↓↑↑↓⟩+|↓↑↓↑⟩),
    
    |S^'⟩=1/2(|↑↓↑↓⟩-|↑↓↓↑⟩-|↓↑↑↓⟩+|↓↑↓↑⟩).

Note that in Eqs. (<ref>) and (<ref>) the specific choice of basis states has been used to match the state |S⟩ generated in the protocol described in Fig. 4 in the main text and in the following section.












§ PROTOCOL FOR THE GENERATION OF THE SINGLET QUBIT STATES WITH S=0



Compared to the state |D_0⟩ the singlet states |s=0,m_x=0⟩, defined by S⃗^2|s=0, m_x=0⟩=S_z |s=0, m_x=0⟩=S_x|s=0, m_x=0⟩=0, are not directly adiabatically connected to any of the bare states. Nevertheless these states can still be prepared by using an adapted protocol (Fig. 4 of the main text) that we are going to describe here in more detail. 

Similar to the ground-state protocol, we start from the decoupled regime where g≃ 0 and ω_q^i≫ω_r, but we initialize the system in the excited qubit state |Ψ_0⟩=|0⟩⊗|↑↑↓↓⟩ (|Ψ_0⟩=|0⟩⊗|↑↓⟩ for N=2 ), where half of the qubits are in the excited state and half in the ground state. Note that for N=4 and a fully symmetric system, the s=0 manifold is two-fold degenerate and spanned, e.g., by the basis states given in Eq. (<ref>). To prepare a well-defined state, we break the symmetry by creating an offset between the qubit frequencies, for example, by setting ω_q^1,2ω_q^3,4. Once the state |Ψ_0⟩ is prepared, all the qubit frequencies are lowered below the resonator frequency such that ω_q^i<ω_r/2. This is done in time step T_1 while keeping g≃ 0. As shown in Fig. <ref>, after this initial step all the relevant qubit states are below the first excited photon state. This configuration avoids undesired level crossings with higher photon number states during the next step of the protocol and only the n=0 manifold must be considered.  
  
During the second step ω_q^i≤ω_r/2, but we still keep a finite frequency difference between the qubits to separate the state |↑↑↓↓⟩ from other states with two qubits excited. This difference between the degenerate and non-degenerate qubit frequencies is visualized by  Fig. <ref>(a) and (b). As the coupling g is slowly increased while the difference in the qubit frequencies is tuned to zero, the state |↑↑↓↓⟩ is adiabatically transformed into the s=0 state |S⟩. During this process the state |S⟩ become almost completely degenerate with the other s=0 state |S^'⟩ [see Fig. <ref>(b)]. However, also the non-adiabatic coupling between these two states is almost negligible, such that the preparation process is still adiabatic on the timescale of the protocol. 
    
In the last step of the protocol, the qubit frequencies are ramped up to the initial values as shown in Fig. 4 of the main text. At this point maintaining a frequency offset is not  crucial anymore. Note that in this last protocol step there are not restrictions on the operational time T_3 because the system is now in the dark state |S⟩ and completely decoupled from the resonator mode. 







§ DISORDER

In the main text we have assumed that our qubits are perfectly identical with matching frequencies, and that each of them couples to the resonator with the same coupling constant. However, due to fabrication disorder and control imprecisions this assumption can be difficult to achieve in experiments with multiple qubits. To evaluate the influence of disorder on the entanglement-harvesting protocol, we show in Fig. <ref> the results obtained from numerical simulations of the protocol in the presence of frequency and coupling disorder. Figure <ref>(a) shows the average fidelity, assuming that in each run of the protocol the individual qubit frequencies evolve as ω_q^i(t)= ω_q(t)(1+ϵ_i), where ω_q(t) follows the ideal pulse given in Fig. 2 of the main text and the ϵ_i are randomly chosen from the interval [-0.1,0.1]. We see that the main part of the protocol is essentially unaffected by frequency disorder, since the system is initially in the ground state and in the USC regime the system is dominated by the interaction terms. Frequency disorder only becomes important in the final decoupled state, where it dephases the symmetric state |D_0⟩. Note, however, that for a fixed frequency distribution, this dephasing can be undone, since as shown in Fig. <ref>(a), it leads to almost no degradation of the purity or the degree of entanglement of the qubit state. In Fig. <ref>(c) and (d) the same plots are shown for the case of coupling disorder g_i(t)=g(t)(1+ϵ_i). Although, this type of disorder has a stronger influence on the evolution of the state, the plot shows that our protocol does not require strictly identical couplings and variation of around 10% still lead to EEF ℱ_ E≳ 0.9 and almost no degradation of the qubit-qubit entanglement. The main quantity affected is the entanglement entropy of the qubit subsystem, which does not approach the value of zero, thus showing that qubits and resonator are not perfectly decoupled. However, we note that this measure of entanglement is very sensitive in our case, since the qubit state we achieve at the end of the protocol coincides with our target state with fidelity above 90%.







§ IMPLEMENTATION OF THE PROTOCOL

In this section we describe a specific flux-qubit circuit, which can be operated in the USC regime and allows a high tunability of the qubit frequencies and couplings. We propose to achieve this goal by using four-junction flux qubits <cit.> with two of the junctions replaced by a SQUID-loop, effectively turning them into junctions with a flux-tunable Josephson energy. The flux qubit design is depicted in Fig. <ref>.

The tunable junctions are junctions 2 and 4 which have sizes α and β, respectively, with respect to junctions 1 and 3. We denote by φ_i = Φ_i/Φ_0 the jump of the superconducting phase across the junction i, where Φ_0 = ħ/2e is the reduced flux quantum and Φ_i the jump of the generalized flux across the junction. The conjugate charge to φ_i is denoted as n_i. The phases φ_i are not independent but constrained by the flux quantization condition for the three loops (the big loop and the two SQUID-loops α and β)

    ∑_i ∈{1,2,3,4}φ_i + f_ϵ   = 0,
    ∑_i ∈{2,6}φ_i + f_α   = 0,
    ∑_i ∈{4,5}φ_i + f_β   = 0,

where f_η = Φ_η/Φ_0 and η = α, β, ϵ is the magnetic frustration through the loop created by external magnetic fluxes Φ_η. Using the above equations we eliminate phase jumps φ_2, φ_5 and φ_6 from the problem. The standard quantization procedure for circuits then gives the Hamiltonian <cit.>

    H_q    = 4E_Cα + β + 2αβ[ (α + β + αβ) (n_1^2 + n_3^2) + (1 + 2α)n_4^2 - 2αβ n_1n_3 - 2α(n_1 + n_3)n_4 ]
           - E_J [ cos(φ_1) + αcos(f_α2)cos(φ_1 + φ_3 + φ̃_4 + f̃_ϵ) + cos(φ_3) + βcos(f_β2)cos(φ̃_4) ],

where φ̃_4 = φ_4 - f_β/2 and f̃_ϵ = f_ϵ + (f_β - f_α)/2. From the shape of the Hamiltonian we can see that if we tune the frustration parameters, f_α, f_β and f_ϵ, in unison such that f_ϵ = (2π + f_α - f_β)/2, the structure of the Hamiltonian stays the same except that the effective Josephson energy of the SQUID-loops vary sinusoidally with the frustration. This enables us to operate the flux qubits at the sweet spot, f̃_ϵ = π, while changing the potential landscape. Note that in practice the cross-talk between the magnetic fluxes may complicate the qubit control, but in principle it is always possible to measure this cross-talk and compensate it by appropriately chosen control pulses. 


The flux qubit couples to the resonator through the phase jump over the entire qubit (see the main text), which, with our notation, is given by Δφ = φ̃_4. The coupling constant g between the resonator and the qubits is proportional to the matrix element of Δφ between the ground and excited states of the qubits, Δφ_eg = ⟨ e | Δφ | g ⟩. Additionally, the coupling to the resonator renormalizes the qubit Hamiltonian by adding a term E_LΔφ^2/2, where E_L = Φ_0^2/L is the inductive energy related to the resonator inductance L, to the qubit Hamiltonian.

Now we are ready to demonstrate the tunability of the qubit frequency and qubit-resonator coupling. We diagonalize the qubit Hamiltonian H_q, plus the renormalization term coming from the coupling, numerically to find the eigenfrequencies and evaluate the transition matrix element Δφ_eg. We choose the following parameters for the simulation: α = 0.6, β = 6, E_L/h = 2.57 GHz, E_C/h = 4.99 GHz and E_J/h = 99.7 GHz. Our choice of E_L sets the resonator inductance to L = 63.7 nH. In addition we choose C = 1.59 pF which determines the resonator frequency and impedance to be ω_r = (LC)^-1/2 =  2π×500 MHz and Z_r = √(L/C) = 200 Ω, respectively. In the simulation we tune f_α from 0 to 0.70π and f_β from 0 to 0.96π (f_ϵ changes accordingly to keep the qubit at its sweet spot).



In Fig. <ref> a) we plot the transition frequency of the qubit, normalized to the resonator frequency against the external fluxes in the two SQUID-loops. The qubit frequency is highly tunable ranging from ∼ 50 ω_r all the way to ∼ 0.5 ω_r. The tunability of the normalized coupling constant, g/ω_r, is demonstrated in Fig. <ref> b). It is also highly tunable ranging from ∼ 0.17 up to ∼ 4.5. The only limitation we have here is that the qubit transition frequency and qubit resonator coupling cannot be tuned entirely separately. The path that we propose to take in the (f_α,  f_β) landscape is portrayed in Fig. <ref> with the red curves. We start from the point (0,  0) and then traverse the curve clockwise, as the arrows show, back to the origin. In Fig. <ref> a) we plot the proposed time-dependent pulses for g(t) and ω_q(t). We have chosen the shape of the pulse for g and from that computed the pulse shape of ω_q. The pulse sequence consists of two parts: ramping up g/ω_r from 0.25 to 4.5 adiabatically and then tuning it back down to its initial value non-adiabatically. The qubit frequency starts from 22.8ω_r, goes down to 0.7ω_r and then ramps back up to its initial value. Fig. <ref> b) shows the fidelity obtained with the pulse of a) in the case of four qubits coupled to the LC-resonator. Even with a limited individual tunability of g and ω_q we obtain fidelities of ℱ_ E > 0.9. For two qubits we can obtain a fidelity of ℱ_ E≈ 0.96 with the same pulse.





§ INFLUENCE OF HIGHER RESONATOR MODES

In our analysis in the main text we have assumed a single mode lumped-element resonator and have neglected the influence of all higher modes. For a typical lumped-element resonator of dimension d=500μm and a resonance frequency of 2π× 5 GHz, the next higher mode is estimated to be at ω_e≈ c/(6d)∼2π× 100 GHz. In Ref. <cit.> we have shown that for such a high ratio ω_e/ω_r > 20, the effect of a higher circuit mode has no relevant influence on the resulting USC physics. By simply scaling this circuit by a factor 10, we would obtain a fundamental resonance of ω_r/2π≈ 500 MHz and an excited mode at around ω_e/2π≈ 10 GHz. By using an optimized design to increase this frequency or simply changing the maximal qubit frequency to, e.g., 8 GHz would avoid a resonant coupling to this mode without affecting the protocol. 

During the first stage of the protocol the system is in the ground state. There might be some weak admixing with the excited mode ∼(g/ω_e)^2, which however, should not affect the adiabaticity condition. During the USC part of the protocol, the qubit frequency is tuned to a very low value. In this case the relation between g, ω_q and ω_r is very similar to the values assumed in the analysis in Ref. <cit.>. Therefore, according to this study the system properties should not be considerably affected. Finally, during the last stage of the protocol, where the qubit frequency is ramped-up again, the coupling has already been switched back to a very small value. The mixing with the excited mode will be small, ∼ g^2_ min/(ω_e-ω_q)^2. 

In summary, based on these estimates we do not expect a significant degrading of our protocol from interactions with higher-order modes of the lumped-element resonator.



§ NUMERICAL SIMULATIONS



 §.§ Coherent evolution

In this short paragraph we provide some details about the numerical simulations, which have been used to produce the plots in the main text.

For the plot of the eigenvalues in Fig. 1 in the main text we have diagonalized Hamiltonian (3) using a truncated set of 140 number states for the resonator mode. For the implementation of the entanglement harvesting protocols shown in Figs. 2, 3 and 4 in the main text we have numerically integrated the time-dependent Schrödinger equation using 100 resonator states. In all calculations we have verified that increasing this number of basis states does not affect our results. 

In Fig. 3(c) in the main text we plot the EEF for N = 4 and for a resonator mode initially in a thermal state ρ_ th=∑_n p(n)|n⟩⟨n| with p(n) the Gibbs distribution at temperature T. The fidelity has been obtained by solving the Schrödinger equation for each initial number state separately, and then averaging all the resulting fidelities according to the thermal probabilities, p_n.
For this plot we have included the first 10 resonator Fock-states and verified that averaging over a thermal distribution including more resonator states does not change the result.




 §.§ Master equation simulations

To make sure dissipation during the protocol does not spoil our scheme we perform simulations including cavity decay into the dynamics. We do so by modelling the coupling to the bath by a Markovian master equation. Note, however,  that in the ultrastrong coupling regime the dissipator is not given by the photon annihilation operator a and must be expressed in terms of the coupled eigenbasis states |k⟩  of the system Hamiltonian <cit.>. For the case of the extended Dicke model Hamiltonian and assuming an ohmic spectral density of the bath, the resulting master equation reads


    ρ̇ = -iħ[H, ρ] + ∑_k,  l > kΓ_kl(1 + n̅(Δ_lk, T)) 𝒟[|k⟩⟨l|]ρ + ∑_k,  l > kΓ_kln̅(Δ_lk, T)) 𝒟[|l⟩⟨k|]ρ,

where Γ_kl = κΔ_lkω_r |C^a_kl|^2, Δ_lk=E_l-E_k  is the energy difference between eigenstates l and k, C^a_kl = ⟨k| a + a^†|l⟩ and κ = ω_r/Q is the resonator decay rate in the weak coupling regime. n̅(Δ_lk, T) is the occupation of the environmental modes at the transition frequency Δ_lk at temperature T. The superoperator 𝒟 is defined in the standard way as 𝒟[O]∙ =(2O ∙ O^† - { O^†O, ∙})/2. Because we change ω_q and g in time the basis {|k⟩}, and thus also the jump operators and rates, changes in time. 
For the numerical simulation of Eq. (<ref>) we diagonalize the system Hamiltonian at each point in time to construct the correct jump operators and transition rates during the different stages of the protocol. After constructing the correct operators and transition rates we move back to the composite basis |n, s, m_z⟩ spanned by the Fock states |n⟩ and collective spin states |s, m_z⟩. This is done with a time-dependent unitary transformation U(t) = ∑_k, n, m_z|n, s, m_z⟩⟨k(t)|, where |k(t)⟩ is the instantaneous eigenstate of the extended Dicke Hamiltonian.



Using the above prescription we simulate our protocol in presence of photon decay for a two qubit system. In Fig. <ref> we plot the fidelity for an LC resonator of Q-factor Q=100 in contact with a bath at temperature T = 0, <ref>a), and T = ħω_r/k_B, <ref>b). We see that even for such a low Q-factor the cavity decay has no significant effect. At T = 0, where the system is essentially at all times in the ground state, the performance of the protocol is not affected at all and even at T = ħω_r/k_B the degradation is only on the level of ∼ 3.5%. In Fig. <ref>b) the dashed green line indicates the performance of the protocol at T = ħω_r/k_B without including cavity decay.


99
Arecchi1972s F. T. Arecchi, E. Courtens, R. Gilmore, and H. Thomas, “Atomic Coherent States in Quantum Optics", Phys. Rev. A 6, 2211 (1972).
Qiu2016s Y. Qiu, W. Xiong, X.-L. He, T.-F. Li, and J. Q. You,  “Four-junction superconducting circuit", Sci. Rep. 6, 28622 (2016).
vool16s  U. Vool and M. Devoret, “Introduction to Quantum Electromagnetic Circuits", arXiv:1610.03438 (2016).
Jaako2016s T. Jaako, Z.-L. Xiang, J.J. Garcia-Ripoll and P. Rabl, “Ultrastrong coupling phenomena beyond the Dicke model", Phys. Rev. A 94, 033850 (2016).
beaudoin11s F. Beaudoin, J. M. Gambetta, and A. Blais, “Dissipation and ultrastrong coupling in circuit QED", Phys. Rev. A 84, 043832 (2012).
ridolfo12s A. Ridolfo, M. Leib, S. Savasta, and M. J. Hartmann. “Photon Blockade in the Ultrastrong Coupling Regime", Phys. Rev. Lett. 109, 193602 (2012).


