Harvesting Multiqubit Entanglement from Ultrastrong Interactions in Circuit Quantum Electrodynamics

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.

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.
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 [1][2][3][4][5], where interaction strengths comparable to the photon energy have been demonstrated. In particular, in the field of circuit QED [6,7], a single superconducting two-level system can already be coupled ultrastrongly [8][9][10]] to a microwave resonator mode [11][12][13][14][15][16][17]. In this regime the physics changes drastically and even in the ground state various nontrivial effects like spontaneous vacuum polarization [18][19][20], light-matter decoupling [21,22] and different degrees of entanglement [22][23][24][25] 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 [26][27][28][29][30][31].
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 [22]. 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 [32][33][34][35][36][37][38], which extracts quantum correlations from USC states and converts these correlations into equivalent multipartite en-tangled 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 [22] 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 [39,40] and form, for example, a resource for Heisenberg-limited metrology applications [41]. 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. 1(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 [42,43] 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 (i) q 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. 1(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 arXiv:1707.08969v3 [quant-ph] 2 Nov 2017 and creation operators, we obtain Here σ i k are Pauli operators and ω i q are the qubit-level splittings. The second term in Eq. (2) accounts for the collective qubit-resonator interaction with couplings The condition g i > ω r , ω i q can be reached with an appropriate flux-qubit design [14,15,18,27,46,47], and the g i (t) and ω i q (t) can be individually tuned by controlling the matrix element ϕ i 0 and the height of the tunnel barrier via local magnetic fluxes [27,36]. A specific four-junction qubit design [47,48], which combines strong interactions with a high degree of tunability, is detailed in the Supplemental Material [43]. Finally, the last contribution in Eq. (2) 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 ω i q = ω q . In this case the Hamiltonian (2) can be expressed in terms of collective angular momentum operators S k = i σ i k /2 and reduces to the extended Dicke Hamiltonian [22] 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 2 x dominate and the level structure changes completely. This is illustrated in Fig. 1(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,mx,n e − g ωr (a † −a)Sx |n ⊗ |s, m x , with energies E s,mx,n ω r n + δE (n) s,mx [22]. 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 and the resulting ordering of the states is shown in Fig. 1(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 [22]. 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 2(a) shows a general pulse sequence for implementing the entanglementharvesting 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 = 4 n=1 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 ω −1 max and g −1 max . 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 −1 min . 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. 2 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) F E = max{F(t)|t ≥ T f }, i.e., the maximal fidelity after the decoupling step, reaches near perfect values of F 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 [14,15]. For these parameters, the nonadiabatic switching times assumed in Fig. 2 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 [49]. 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 [50] 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. 3(a), which shows the evolution of the lowest eigenenergies during different stages of the proto-col 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. 3(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 F 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 ω i q (t) due to fabrication uncertainties are taken into account. Numerical simulations of the protocol under such realistic experimental conditions [43] 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. 3(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 B T − 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. 3(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 1(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 (3) 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. 4(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-photonnumber 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 [51,52] [cf. Fig. 1(d)], out of which the protocol selects a specific superposition [43].
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. 4(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.

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 [S1], 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

qubits
For two qubits we have the usual three triplet states and the singlet When expressed in terms of the original qubit basis the two m x = 0 states of interest read

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., 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 Each of the s = 1 states is 3-fold degenerate and the corresponding states are Finally, for N = 4 we obtain two singlet states with s = 0, which are given by Since the s = 0 subspace is invariant under rotations, the states have the same form as in the original qubit basis, Note that in Eqs. (S8) and (S9) 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 ω i q ω 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. (S9). To prepare a well-defined state, we break the symmetry by creating an offset between the qubit frequencies, for example, by setting ω 1,2 q = ω 3,4 q . Once the state |Ψ 0 is prepared, all the qubit frequencies are lowered below the resonator frequency such that ω i q < ω r /2. This is done in time step T 1 while keeping g 0. As shown in Fig. S1, 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 ω i q ≤ ω 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. S1(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. S1(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. S2 the results obtained from numerical simulations of the protocol in the presence of frequency and coupling disorder. Figure S2(a) shows the average fidelity, assuming that in each run of the protocol the individual qubit frequencies evolve as ω i q (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. S2(a), it leads to almost no degradation of the purity or the degree of entanglement of the qubit state. In Fig. S2(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 F 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 [S2] with two of the junctions replaced by a SQUID-loop, effectively turning them into junctions with a fluxtunable Josephson energy. The flux qubit design is depicted in Fig. S3. 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} 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 [S2, S3] 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 = Φ 2 0 /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. S4 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. S4 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. S4 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. S5 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. S5 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 F E > 0.9. For two qubits we can obtain a fidelity of F 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. [S4] 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. [S4]. 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.

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 [S5, S6].
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 D is defined in the standard way as