Ultrastrong coupling circuit QED in the radio-frequency regime

TL;DR

This paper explores ultrastrong coupling in circuit QED with multiple qubits and a radio-frequency resonator, revealing novel symmetry-breaking transitions and particle-like quantum phenomena accessible through spectroscopic and Ramsey measurements.

Contribution

It demonstrates how ultrastrong coupling leads to new symmetry-breaking transitions and effective particle dynamics in circuit QED, expanding the understanding of quantum phenomena in these systems.

Findings

Observation of symmetry-breaking transitions in ground and excited states.

Probing wave packet dynamics via GHz spectroscopy and Ramsey measurements.

Identification of effective particle behavior in ultrastrong coupling regimes.

Abstract

We study a circuit QED setup where multiple superconducting qubits are ultrastrongly coupled to a single radio-frequency resonator. In this extreme parameter regime of cavity QED the dynamics of the electromagnetic mode is very slow compared to all other relevant timescales and can be described as an effective particle moving in an adiabatic energy landscape defined by the qubits. The focus of this work is placed on settings with two or multiple qubits, where different types of symmetry-breaking transitions in the ground- and excited-state potentials can occur. Specifically, we show how the change in the level structure and the wave packet dynamics associated with these transition points can be probed via conventional excitation spectra and Ramsey measurements performed at GHz frequencies. More generally, this analysis demonstrates that state-of-the-art circuit QED systems can be used to access a whole range of particle-like quantum mechanical phenomena beyond the usual paradigm of coupled qubits and oscillators.

Tables

Table 1. Table 1: The set of circuit parameters for the qubits and the L ​ C 𝐿 𝐶 LC resonator, which are used for the numerical simulations shown in Fig.  7 (b).

Qubit	Resonator
$C_J=2.21\mathrm{fF}$	$C=79.58\mathrm{pF}$
$E_J/h=336.8\mathrm{GHz}$	$L=127.3\mathrm{nH}$
$\alpha =0.74$	$C_s=1.06\mathrm{fF}$
	$L_s=1.27\mathrm{nH}$

Equations

$$\begin{array}{rl}{H}_{\mathrm{cQED}}={\omega }_r{a}^{†}a+& \sum _{i=1}^N\frac{{\omega }_q^i}{2}{\sigma }_z^i+\sum _{i=1}^N\frac{g_i}{2}(a^{†}+a){\sigma }_x^i\\ +& \sum _{i,j=1}^N\left(\frac{g_i{g}_j}{4{\omega }_r}+J_{ij}\right){\sigma }_x^i{\sigma }_x^j.\end{array}$$

$$H_{\mathrm{mag}}=\frac{({\varphi }_r-{\varphi }_q{)}^2}{2L}.$$

$$\begin{array}{r}{H}_{\mathrm{EDM}}={\omega }_r{a}^{†}a+{\omega }_q{S}_z+g(a^{†}+a)S_x+(1+\epsilon )\frac{g^2}{{\omega }_r}{S}_x^2,\end{array}$$

$$X=\sqrt{\frac{{\omega }_r}{2{\omega }_q}}(a^{†}+a),P=i\sqrt{\frac{{\omega }_q}{2{\omega }_r}}(a^{†}-a),$$

$${\tilde{H}}_{\mathrm{EDM}}=\frac{P^2}{2\mu }+\frac{X^2}{2}+{\tilde{H}}_q(X),$$

$${\tilde{H}}_q(X)=S_z+\sqrt{2}\lambda X{S}_x+(1+\epsilon ){\lambda }^2{S}_x^2.$$

$${|\psi ⟩}_{n,k}=\int dX{\varphi }_{n,k}(X)|X⟩|{\chi }_n(X)⟩,$$

$${\tilde{H}}_q(X)|{\chi }_n(X)⟩={\tilde{E}}_n(X)|{\chi }_n(X)⟩$$

$$\left[-\frac{1}{2\mu }\frac{{\partial }^2}{\partial X^2}+{\tilde{V}}_n(X)\right]{\varphi }_{n,k}(X)={\widetilde{ϵ}}_{n,k}{\varphi }_{n,k}(X).$$

$$C_{n,m}=\frac{1}{\mu }\frac{\partial {\varphi }_{n,k}}{\partial X}\frac{⟨{\chi }_m(X)|\sqrt{2}\lambda S_x|{\chi }_n(X)⟩}{{\tilde{E}}_n(X)-{\tilde{E}}_m(X)},$$

$${\tilde{V}}_{s,m_z}^{(\mathrm{DM})}(X)=\frac{X^2}{2}+m_z\sqrt{1+2{\lambda }^2{X}^2},$$

$$\Delta \widetilde{ϵ}{|}_{\lambda >{\lambda }_c}\approx \frac{8}{\sqrt{\pi }}Ae{}^{-S_0},$$

$$\Delta \widetilde{ϵ}{|}_{\lambda ={\lambda }_c}=c\sqrt[3]{\frac{{\lambda }_c^2}{{\mu }^2}},$$

$$\Delta \tilde{E}(X=0)=(1+\epsilon ){\lambda }^2$$

$${\tilde{V}}_{s,m_x}^{(0)}(X)=\frac{1}{2}{\left(X+\sqrt{2}\lambda m_x\right)}^2+\epsilon {\lambda }^2{m}_x^2,$$

$$H_{\mathrm{eff}}$$

$$+\left[(1+\epsilon )\frac{g^2}{2{\omega }_r}-\frac{g^2}{4{\omega }_q}\right]({\vec{S}}^2-S_z^2),$$

$$S(\omega )=\frac{\Gamma }{2}\mathrm{Re}{\int }_0^{\infty }d\tau ⟨{\sigma }_x^1(\tau ){\sigma }_x^1(0)⟩e{}^{i\omega \tau },$$

$$S(\omega )=\sum _f\frac{1}{4}\frac{{\Gamma }^2|⟨f|{\sigma }_x^1|\mathrm{GS}⟩{|}^2}{(\omega -{\omega }_f{)}^2+{\Gamma }^2/4},$$

$$g>\sqrt{\Gamma {\omega }_r}$$

$$S(\omega )\simeq C\sum _k\frac{{\Gamma }^2|⟨{\varphi }_{1,k}|{\varphi }_{0,0}⟩{|}^2}{(\omega -{\omega }_{k0}{)}^2+{\Gamma }^2/4},$$

$$H_{\Omega }(t)=\Omega \cos ({\omega }_{d}t+\theta ){\sigma }_x^1.$$

$$H=H_{LC}+H_{\mathrm{int}}+H_q.$$

$$\begin{array}{rl}{H}_{LC}=& \frac{Q_a^2}{2C}+\frac{Q_b^2}{2C_b}+\frac{({\varphi }_a-{\varphi }_b{)}^2}{2L_s}+\frac{{\varphi }_b^2}{2L}\end{array}$$

$$\begin{array}{r}{H}_{\mathrm{int}}=\frac{Q_b{Q}_q}{C_q}+\frac{{\varphi }_q{\varphi }_b}{L}\end{array}$$

$$H_{\mathrm{int}}=\sum _{\eta =\pm }i\frac{g_{Q,\eta }}{2}(c_{\eta }^{†}-c_{\eta })n_q+\frac{g_{\varphi ,\eta }}{2}(c_{\eta }^{†}+c_{\eta }){\phi }_q.$$

$$\frac{g_{\varphi ,-}}{{\omega }_{-}}=\sqrt{\frac{R_Q}{\pi Z}},\frac{g_{Q,-}}{{\omega }_{-}}=2\sqrt{\frac{\pi Z}{R_Q}}\frac{C_b}{C_q}$$

$$\frac{g_{\varphi ,+}}{{\omega }_{+}}=\sqrt{\frac{R_Q}{\pi Z_b}}\frac{L_s}{L},\frac{g_{Q,+}}{{\omega }_{+}}=2\sqrt{\frac{\pi Z_b}{R_Q}}\frac{C_b}{C_q}$$

$$\tilde{H}=\frac{P^2}{2\mu }+\frac{X^2}{2}+{\tilde{H}}_q(X,P),$$

$${\tilde{H}}_q(X,P)={\tilde{\omega }}_{+}{c}_{+}^{†}{c}_{+}+{\tilde{H}}_{\mathrm{int}}(X,P)+{\tilde{H}}_q,$$

$${\tilde{H}}_0$$

Full text

I Introduction

Circuit quantum electrodynamics (QED) is a rapidly developing field where fundamental processes of quantum light-matter interactions are studied with superconducting qubits (‘artificial atoms’) coupled to microwave resonators and transmission lines Blais2004 ; Wallraff2004 ; Gu2017 . Due to the extraordinary combination of strong coupling and very low losses, many quantum optical phenomena, such as vacuum Rabi splittings Wallraff2004 , photon blockade Lang2011 ; Hoffman2011 or super- and subradiant decay mlynek2014 , have already been demonstrated in these systems with very high precision. Moreover, by using high-impedance resonators or by employing galvanic instead of electric coupling schemes, circuit QED systems can overcome fundamental bounds on the coupling strength in conventional cavity QED systems Devoret2007 ; debernardis2018 . It is then possible to access so-called ultrastrong (USC) or deep-strong coupling regimes Ciuti2005 ; Casanova2010 ; forndiaz18 ; kockum19 , where the qubit-photon coupling is comparable to the photon energy and light-matter interactions become non-perturbative. These conditions have recently been demonstrated in experiments with superconducting qubits coupled to microwave resonators and transmission lines forndiaz10 ; baust16 ; forndiaz17 ; yoshihara17 ; Chen2017 ; Bosman2017 . When extended to multiple qubits, this regime could enable new applications such as protected quantum memories nataff11 , ultrafast gate operations romero12 or entanglement harvesting schemes sabin2012 ; Armata2017 .

In regular cavity- and circuit QED systems with weak or moderate coupling, light-matter interactions are only effective close to resonance, where energy-conserving transitions between photons and atomic excitations can take place. This constraint does no longer apply in the USC regime, where photons and qubits with very dissimilar frequencies can still strongly influence each other. One specific limit of interest in this context is the low-mode-frequency or adiabatic regime Graham1984 ; Crisp1992 ; Liberti2006 ; liberti2006b ; Ashhab2010 , where the bare oscillation frequency of the resonator mode, $\omega_{r}$, is much smaller than the qubit transition frequency, $\omega_{q}$. In this regime, the qubit state adjusts instantaneously to the slowly varying field amplitude and provides in turn an effective adiabatic potential for the photon mode. This separation of timescales is similar to the Born-Oppenheimer (BO) approximation in the description of nucleus-electron systems in molecular physics born1927 . In the weak-coupling limit, this physics finds applications, for example, for the readout of superconducting qubits or quantum dots through the off-resonant coupling to a low-frequency mode johansson2006 ; nigg2009 ; lahay2009 . In the USC regime, this adiabatic picture is frequently employed in quantum optics and solid state physics to discuss symmetry-breaking effects in Rabi-, Dicke- and Jahn-Teller models, where the ground-state potential surface changes from a shape with a single minimum to a double-well or mexican-hat potential larson2017 ; bersuker2006 . In Refs. Bakemeier2012 ; Ashhab2013 ; hwang2015 ; puebla2016 ; hwang2018 it has been shown in more detail that even for a single qubit this change in the adiabatic potential reproduces the properties of quantum-, excited-state- or dissipative phase transitions when $\omega_{r}\rightarrow 0$. However, reaching this regime with circuit QED systems requires electromagnetic modes in the radio-frequency regime where $\omega_{r}$ is only a few or a few tens of MHz. At these frequencies the mode is thermally occupied even at mK temperatures and electronic measurement techniques, which operate efficiently only above $\sim 1\,\mathrm{GHz}$, are not readily available. There is, nevertheless, experimental progress towards the quantum control of such modes. For example, active ground state cooling and the readout of individual photon number states of a $170\,\mathrm{MHz}$ electromagnetic resonator has recently been demonstrated gely2019 .

In this paper we study the properties of circuit QED systems in the low-mode-frequency regime by considering a generic setup of two or multiple flux qubits coupled to a radio-frequency resonator. We derive an effective description of this system in terms of adiabatic BO potentials for the electromagnetic mode and discuss the characteristics of the resulting potential energy surfaces for different qubit numbers and coupling parameters. Compared to related previous works Graham1984 ; Crisp1992 ; Liberti2006 ; liberti2006b ; Ashhab2010 ; larson2017 ; Bakemeier2012 ; Ashhab2013 ; hwang2015 ; puebla2016 ; hwang2018 , our main interest here is in the excited potential curves. These potentials can exhibit multiple first- and second-order symmetry-breaking transitions, even under conditions where the ground state potential still has a single minimum. This leads to a qualitatively new situation where properties of the symmetry-breaking transition occurring at MHz frequencies can be probed via regular qubit-readout techniques operated in the GHz domain. As two specific examples, we describe how the change in the level structure and the dynamics of the wave packet splitting near the transition point can be detected via excitation spectra and Ramsey coherence measurements. This analysis demonstrates that quantum dynamics and phase-transition physics in the radio-frequency regime can be controlled and detected using state-of-the-art superconducting-circuit technology.

The remainder of the paper is structured as follows. After introducing in Sec. II the basic model and its approximate treatment in the low-frequency limit, we provide in Sec. III a general overview of the characteristic features that can appear in the resulting adiabatic potential curves. In Sec. IV and Sec. V we then discuss two schemes for detecting the symmetry-breaking transition in the excited potential. Finally, in Sec. VI we perform a more rigorous justification for the single-mode approximation in a realistic setting and conclude our findings in Sec. VII.

II Circuit QED in the radio-frequency regime

We consider a circuit QED setup as shown in Fig. 1, where $N$ flux qubits are coupled inductively to a lumped-element resonator with frequency $\omega_{r}$ and bosonic annihilation (creation) operator $a$ $(a^{\dagger})$. By taking only the two energetically lowest states $|g\rangle$ and $|e\rangle$ of each flux qubit into account, the quantized dynamics of this circuit is described by the general cavity QED Hamiltonian jaako2016 ; debernardis2018 ($\hbar=1$)

[TABLE]

Here $\sigma_{x,y,z}$ are the usual Pauli operators and $\omega^{i}_{q}$ and $g_{i}$ denote the transition frequency and the coupling strength of each qubit. Apart from the collective qubit-photon coupling in the first line of Eq. (1), this Hamiltonian also contains two direct qubit-qubit interaction terms. The so-called depolarization term $\sim g_{i}g_{j}$ represents the gauge-dependent part of the qubit-photon interaction debernardis2018 ; debernardis2018b and therefore should not be interpreted as a physical coupling between the qubits. In the current setup its origin can be understood as follows. When expressed in terms of the generalized flux variables indicated in Fig. 1(a), the magnetic energy stored in the inductor reads

[TABLE]

Here $\phi_{r}\sim(a+a^{\dagger})$ is the oscillator flux variable and $\phi_{q}=\sum_{i=1}^{N}\Delta\phi_{i}$, where the $\Delta\phi_{i}\sim\sigma_{x}^{i}$ represent the flux jumps across each qubit. Therefore, after expanding this energy in terms of the canonical variables, we obtain the photon-qubit interaction together with an apparent all-to-all coupling between the qubits jaako2016 ; Armata2017 .

Equation (1) also includes additional direct qubit-qubit couplings with strengths $J_{ij}$. Such couplings arise, for example, from the mutual inductance between nearby flux qubits, or, more generally, from a common coupling of the flux qubits to auxiliary superconducting quantum interference devices (SQUIDs), see, for example, Ref. Plourde2004 . In the latter case, the range, the sign and the strength of the elements $J_{ij}$ can be fully engineered and controlled by external bias currents. Although the presence of such qubit-qubit couplings is not essential for the main effects described in this work, they provide an additional tuning nob for the resulting potential surfaces discussed below.

II.1 Extended Dicke model

For concreteness and to simplify the discussion, we will focus in the remainder of our analysis on the case of identical qubits, $\omega_{q}^{i}=\omega_{q}$ and $g_{i}=g$, and all-to-all inter-qubit interactions, $J_{ij}=J$. In this case, Hamiltonian (1) reduced to the extended Dicke model (EDM) jaako2016 ; Todorov2014

[TABLE]

where $S_{\alpha}=1/2\sum_{i=1}^{N}\sigma_{\alpha}^{i}$ are collective spin operators. Here we have adopted the convention $J=\varepsilon g^{2}/(4\omega_{r})$ debernardis2018 , such that the dimensionless parameter $\varepsilon$ characterizes the relative strength between qubit-qubit and qubit-photon interactions. We emphasize that none of the qualitative findings of this work rely on the assumptions of identical qubits or purely collective couplings. We assume, however, that even in the presence of imperfections the model preserves its parity symmetry, i.e., it remains invariant under the transformation $\sigma_{x}^{i}\rightarrow-\sigma_{x}^{i}$ and $a\rightarrow-a$. For a more detailed derivation of the EDM for two basic circuit configurations, which correspond to $\varepsilon\geq 0$ and $\varepsilon=-1$, the reader is referred to Refs. jaako2016 ; Armata2017 and Bamba2016 , respectively. Further details on the validity of the single-mode approximation assumed in our model are postponed to Sec. VI.

II.2 Born-Oppenheimer approximation

In regular cavity and circuit QED setups the coupling $g$ between (artificial) atoms and photons is usually small and only resonant processes, where $\omega_{r}\approx\omega_{q}$, are relevant. In this regime, the physics of light-matter interactions is most conveniently described in terms of the bare photon number eigenstates, $|n\rangle$, where $a^{\dagger}a|n\rangle=n|n\rangle$. Close to resonance, these photons can hybridize with matter excitations, which results, for example, in the appearance of a vacuum Rabi-splitting $\sim g\sqrt{N}$ in the excitation spectrum. In this work we are interested in a very different parameter regime, $\omega_{r}\ll\omega_{q},g$ and $g^{2}/\omega_{r}\sim\omega_{q}$. Under such conditions the oscillation of the electromagnetic mode is slow compared to the qubit dynamics, while at the same time the coupling to the qubits has a substantial effect on the cavity and vice versa. Therefore, the representation of the electromagnetic field in terms of the bare photon number states is no longer useful in this regime. Instead, it is more convenient to describe the electromagnetic field as an effective massive particle moving in a set of adiabatic potential surfaces.

To model the static and dynamical properties of the circuit QED system in the limit $\omega_{r}\rightarrow 0$, we start by introducing the rescaled quadrature variables

[TABLE]

which correspond to the usual position and momentum operators of a harmonic oscillator. After normalizing all energies with respect to the qubit frequency, i.e., $\tilde{H}_{\rm EDM}=H_{\rm EDM}/\omega_{q}$, we obtain

[TABLE]

where $\mu=\omega^{2}_{q}/\omega^{2}_{r}$ is the effective mass and

[TABLE]

Here we have defined $\lambda=\sqrt{g^{2}/(\omega_{r}\omega_{q})}$ as the relevant dimensionless coupling parameter.

The decomposition used in Eq. (5) shows that for $\omega_{r}\ll\omega_{q}$ the ‘kinetic’ energy term $\sim P^{2}$ is small compared to the potential and qubit energies. Thus, in direct analogy to the treatment of nucleus-electron systems in molecular physics, we can apply a BO approximation to separate the fast dynamics of the qubits and the slow motion of the resonator. Under this approximation the eigenstates of the combined system are given by Liberti2006 ; liberti2006b

[TABLE]

where $\ket{\chi_{n}(X)}$ is the adiabatic qubit eigenstate. It obeys the Schrödinger equation

[TABLE]

for a fixed value of $X$. The dependence of $\tilde{E}_{n}(X)$ on the position coordinate provides an additional effective potential for the resonator wave function $\phi_{n,k}(X)$, which is a solution of the eigenvalue equation

[TABLE]

Here, $\tilde{V}_{n}(X)=X^{2}/2+\tilde{E}_{n}(X)$ is the total BO potential associated with the $n$-th qubit eigenstate. These adiabatic potentials are symmetric around $X=0$ due to the invariance of Eq. (6) under flipping the signs of $X$ and $S_{x}$ simultaneously.

The BO approximation neglects couplings between the adiabatic energy eigenstates, which are induced by the momentum operator. The main off-diagonal correction terms, which couple Eq. (9) for different $n$, are of the form

[TABLE]

where we used the relation $\langle\chi_{m}(X)|\frac{\partial}{\partial_{X}}|\chi_{n}(X)\rangle=\langle\chi_{m}(X)|\frac{\partial\tilde{H}_{q}(X)}{\partial X}|\chi_{n}(X)\rangle/[\tilde{E}_{n}(X)-\tilde{E}_{m}(X)]$. For moderate couplings and electromagnetic eigenstates near the potential minimum, where $\frac{\partial\phi_{n,k}}{\partial X}\sim\sqrt{\mu}$, we obtain a scaling $C_{n,m}\sim\lambda/\sqrt{\mu}$. In all the numerical examples discussed in this work we consider the parameter regime $\mu=10^{4}-10^{6}$, where the adiabatic condition $|C_{n,m}|\ll 1$ is well-satisfied for the relevant potential curves. However, as discussed in Sec. III.3 below, for large couplings, $\lambda>1$, some of the excited potential surfaces are only separated by higher-order avoided crossings. In this case, $|\tilde{E}_{n}(X)-\tilde{E}_{m}(X)|\ll 1$ and corrections beyond the BO approximation may become relevant.

III Adiabatic potential surfaces

The energy landscape formed by all the $V_{n}$ is fully determined by the qubit Hamiltonian $H_{q}(X)$ and depends on the interaction parameters $\lambda$ and $\varepsilon$ as well as the number of qubits, $N$. In this section we summarize the characteristic features of these potentials in different parameter regimes.

III.1 Ground-state symmetry breaking

In Fig. 2(a) and (b) we first plot $V_{n=0,1}(X)$ for the simplest case of a single qubit, where $H_{\rm EDM}$ reduces to the quantum Rabi model. For $\lambda\ll 1$ the original harmonic potentials $\tilde{V}_{0,1}(X)\simeq X^{2}/2\pm 1/2$ are only weakly perturbed and have a single minimum at $X_{\rm min}=0$. For increasing coupling strength the ground state potential becomes shallower until it transitions into a double-well potential with two degenerate minima at $X_{\rm min}\approx\pm\sqrt{\lambda^{4}-\lambda_{c}^{4}}/(\sqrt{2}\lambda\lambda_{c}^{2})$ above the critical coupling $\lambda_{c}=1$. As indicated in Fig. 2(c) for $N=2$, the same qualitative behavior is also found for the $N$-qubit Dicke model (DM), which corresponds to the case $\varepsilon=-1$ in the current notation. In this specific situation the potentials have the simple analytic form Liberti2006 ; larson2017

[TABLE]

and can be labelled by the total spin quantum number $s$ and the spin projection $m_{z}=-s,s+1,\dots,s-1,s$ associated with the qubit states $|\chi_{s,m_{z}}(X=0)\rangle$. This result shows that also for larger $N$ the transition occurs first in the ground-state potential, $m_{z}=-N/2$, but at a reduced coupling parameter $\lambda_{c}=1/\sqrt{N}$.

This change from a single-well to a double-well structure of the adiabatic potential is familiar from studies of Rabi-, Dicke- and Jahn-Teller-type models, where even at larger $\omega_{r}$ this picture explains the observed symmetry breaking in the ground state. Here, symmetry breaking means that for $\lambda>\lambda_{c}$ the tunnel-splitting between the two lowest resonator states, $\Delta\epsilon$, is exponentially suppressed such that in a realistic setting any weak perturbation will randomly localize the system in one of the degenerate minima. Specifically, by approximating the ground-state potential $V_{N/2,-N/2}^{({\rm DM})}(X)$ by an equivalent quartic double-well potential with the same minima and the same barrier height, we obtain (see e.g. Ref. Garg2000 )

[TABLE]

with an exponent $S_{0}=\frac{2}{3}\sqrt{\mu(\lambda^{2}-\lambda_{c}^{2})^{3}(\lambda^{2}+\lambda_{c}^{2})}/(\lambda^{2}\lambda_{c}^{4})$ and prefactor $A=\sqrt[4]{(\lambda^{2}-\lambda_{c}^{2})^{5}/(\mu(\lambda^{2}+\lambda_{c}^{2}))}/(\lambda\lambda_{c}^{2})$. This means that for $\lambda\gg\lambda_{c}$ the tunneling is suppressed by $\Delta\tilde{\epsilon}\sim\operatorname{e}\!^{-2\sqrt{\mu}\lambda^{2}N^{2}/3}$. When passing from the single-well to the double-well configuration the potential becomes purely quartic at $\lambda=\lambda_{c}$, in which case the minimal energy splitting is vranicar2000

[TABLE]

with a numerical prefactor of $c\approx 1.1$. This means that the density of states at the transition point, $\tilde{\nu}_{\lambda=\lambda_{c}}=1/\Delta\tilde{\epsilon}$, scales as $\sim\sqrt[3]{N\mu^{2}}$ and diverges in the classical limit $\mu\rightarrow\infty$. Note, however, that when compared to the density of states of the unperturbed harmonic potential, $\tilde{\nu}_{\lambda=0}=\sqrt{\mu}$, the peak at the transition point is not very pronounced. This illustrates the necessity to use either many qubits or very large ratios $\mu=\omega_{q}^{2}/\omega_{r}^{2}$ to detect sharp experimental signatures associated with this quasi-divergence. It has been estimated that the regime $\mu>10^{4}$ can be reached using effective implementations of the Rabi model in trapped ion systems Puebla2017 . In circuit QED, reaching this parameter regime requires resonator frequencies below $100$ MHz.

III.2 Non-interacting qubits

While the DM has been the subject of many theoretical studies, it only represents a very special class of cavity and circuit QED setups with strong ferroelectric interactions debernardis2018 . In Fig. 2(d-f) we now show the potential curves for another relevant example of two non-interacting (or only weakly interacting) qubits, where $\varepsilon\simeq 0$. In this case we see that as the coupling parameter $\lambda$ increases, the absolute minimum of the ground state potential remains at $X_{\rm min}=0$, while for $\lambda\gtrsim 1$ two additional local minima appear around $X\approx\pm 2$. More important for the current work, already at an intermediate coupling strength of $\lambda=1/\sqrt{2}$ a symmetry-breaking transition occurs in the first excited potential, see App. A. As explained in more detail in Secs. IV and V below, this has important practical implications. The symmetry-breaking effect now occurs at an absolute frequency scale set by the qubit frequency $\sim\omega_{q}$, at which efficient electronic readout techniques are readily available.

Compared to the DM, another qualitative difference is that for $\varepsilon=0$ (as well as for any $\varepsilon\neq-1$) the degeneracy between potential curves of different spin quantum numbers is lifted by the $S_{x}^{2}$ term. For example, for $N=2$ we obtain a splitting of (see App. A),

[TABLE]

between the triplet and singlet potentials at $X=0$. Therefore, even though the oscillator coordinate is zero, there is still an energy penalty for the triplet state compared to the singlet. This somewhat counterintuitive result can be understood by looking at the original magnetic interaction Hamiltonian given in Eq. (2). At $X=0$, i.e. $\phi_{r}=0$, there will be an energy cost for states for which $\phi_{q}\neq 0$, such as the triplet state. For the singlet this contribution vanishes.

III.3 General structure

For a larger number of qubits the energy potential landscapes become considerably more involved and depending on $\lambda$ and $\varepsilon$ the individual potentials curves can exhibit various local and global minima. This is illustrated in Fig. 3(a) and (b) for the example of $N=3$ and $\varepsilon=0.02$. Apart from the formation a double-well structure in the ground- and second excited potential curve, in this case we also obtain a triple-well potential with three degenerate minima at a value of $\lambda^{2}\simeq 1.86$. This situation corresponds to a first-order phase transition point, where the potential at $X=0$ remains stable, but the two outer wells become lower in energy after the transition point [cf. Fig. 3(b)].

For a fully symmetric system the total number of distinct potential curves scales as $(N/2+1)^{2}$ (for an even number of qubits) and as $2^{N}$, if different parameters for each qubit are taken into account. To obtain a better intuition about the basic potential landscapes that can arise, it is instructive to consider the limit of very large coupling, $\lambda\gg 1$, following the analysis presented in Ref. jaako2016 . In this limit, the terms $\sim S_{x}$ and $\sim S_{x}^{2}$ dominate over the bare qubit splitting $\sim S_{z}$. By neglecting the contribution of the $S_{z}$ term completely, the eigenstates $|\chi^{0}_{s,m_{x}}(X)\rangle$ are also eigenstates of $S_{x}$ and can be labelled by the total spin $s$ and its projection along the $x$-axis, $m_{x}$. The corresponding zeroth-order potential curves are

[TABLE]

which are just a set of displaced parabolas indicated by the dashed lines in Fig. 3(c). As detailed in Ref. jaako2016 and App. A, the presence of the $S_{z}$ term in the full model leads to a coupling between the states with $m_{x}$ and $m_{x}\pm 1$. Focusing on the limit of very small values for $\varepsilon$, this coupling lowers, first of all, each well by $\Delta\tilde{V}_{s,m_{x}}=(m_{x}^{2}-s(s+1))/(2\lambda^{2})$, due to second-order couplings to energetically higher potentials. Further, at the crossing points $X\simeq-\lambda(2m_{x}+1)/\sqrt{2}$, the $S_{z}$ terms lifts the degeneracy and leads to an avoided crossing with a splitting of $\tilde{\Delta}_{m_{x},m_{x}+1}=\sqrt{s(s+1)-m_{x}(m_{x}+1)}$. As indicated in Fig. 3(c), the actual adiabatic potentials $V_{n}(X)$ are the resulting connected potential curves. Note that in this picture, which holds for $\lambda\gg 1$, the excited potentials are only separated through higher-order avoided crossings. For general $\lambda$, the potential energy surface can then be understood as a smooth interpolation between the stack of parabolas $\tilde{V}_{s,m_{z}}(X)=m_{z}+X^{2}/2$ for the uncoupled system into the degenerate-well structure depicted in Fig. 3(c).

III.4 Weak-coupling limit

Finally, let us briefly remark on the limit $\lambda\ll 1$, where the system is far away from the instability points. In this case we can expand the adiabatic potential to second order in $\lambda$ to obtain a state-dependent shift of the oscillation frequency, $\omega_{r}^{m_{z}}=\omega_{r}+m_{z}\delta\omega_{r}$ with $\delta\omega_{r}=g^{2}/(2\omega_{q})$. As long as $\delta\omega_{r}<\omega_{r}$, we can make a rotating wave approximation and describe this shift in terms of an effective Hamiltonian

[TABLE]

which corresponds to the usual Stark-shift Hamiltonian for a far detuned cavity QED system. The third term in the first line of Eq. (16) represents a shift of the qubit frequency proportional to the cavity photon number. This type of interaction is frequently encountered in measurement schemes for quantum dots or superconducting qubits johansson2006 ; nigg2009 ; lahay2009 ; Gu2017 and has been exploited in Ref. gely2019 to resolve individual photon number states of a $170\,\mathrm{MHz}$ resonator footnote . At slightly larger couplings or even lower resonator frequencies, where $\delta\omega_{r}\gtrsim\omega_{r}$, the photon number is no longer a conserved quantity and this picture breaks down.

IV Excitation spectrum

In regular circuit QED systems both $\omega_{r}$ and $\omega_{q}$ are in the GHz regime and both the qubit and the resonator mode can be probed using efficient electronic readout techniques. In the RF domain this is no longer possible using conventional methods and therefore all the properties of the electromagnetic mode must be inferred through measurements performed on the qubits. As a first example, we consider in this section the single-qubit excitation spectrum $S(\omega)$, which is proportional to the excited state population of the first qubit when being driven by a weak external field of frequency $\omega\sim\omega_{q}$. The excitation spectrum is then given by

[TABLE]

where the expectation value is taken with respect to the equilibrium state of the systems in the absence of the driving field. In Eq. (17) we have introduced a characteristic decay rate $\Gamma$, with which the correlation function decays. Since in the parameter regime of interest the qubits are not strongly perturbed by the coupling to the resonator mode, $\Gamma$ can simply be taken as the bare decay rate of the excited qubit state. In the numerical examples discussed in this section we assume moderate values of $\Gamma\sim\omega_{r}$, which is still large compared to decoherence rates achieved with state-of-the-art flux qubits Yan2016 .

IV.1 Mode splitting in the low-frequency regime

In a first step we assume that the combined system is prepared close to the absolute ground state $\ket{\rm GS}$ (for example by actively cooling the resonator mode gely2019 ). In this case the excitation spectrum is given by

[TABLE]

where the sum runs over all final states $\bra{f}$, which are separated by a frequency $\omega_{f}$ from the ground state.

In Fig. 4 (a) and (b) we first compare the excitation spectrum of a conventional, i.e., resonant cavity QED system with the case of a circuit QED system in the low-frequency regime. For both plots $N=2$ and we use the more common convention to plot the spectrum as a function of the coupling strength $g$ instead of $\lambda$. In this case the resonant system exhibits the familiar Rabi splitting $\sim\sqrt{2}g$, which arises from the hybridization of the photon with the triplet state $|T\rangle$. This splitting is visible for $\sqrt{2}g>\Gamma$, which defines the strong coupling regime for a resonant cavity QED system. The additional resonance in the middle is the singlet state $|S\rangle$. The singlet state is decoupled from the cavity mode, but its energy relative to the ground state still decreases for larger coupling strengths.

In the radio-frequency regime, where $\omega_{r}/\omega_{q}=0.01$, the mode frequency is no longer resolved on the scale of this plot. However, the cavity mode has still a drastic influence on the spectrum through the induced splitting between the triplet and the singlet state [see Eq. (14)]. To be observable, this splitting must exceed the qubit decay rate $\Gamma$ and therefore we identify

[TABLE]

as the minimal strong-coupling condition for the low-frequency regime.

IV.2 Probing the symmetry-breaking transition

The dashed lines in Fig. 4(b) indicate the excitation frequencies derived from the energies $V_{n}(X=0)$. The perfect agreement demonstrates that the dominant peaks of $S(\omega)$ at low temperatures probe the adiabatic potentials curves at $X=0$. This can be understood in terms of the BO wave function in Eq. (7) and the fact that for small and moderate couplings the qubit states $|\chi_{n}(X)\rangle$ are only weakly dependent on $X$. The excitation spectrum can then be approximated as

[TABLE]

with a constant prefactor $C=|\bra{\chi_{1}(0)}\sigma^{1}_{x}\ket{\chi_{0}(0)}|^{2}/4$ and $\omega_{k0}=\epsilon_{1,k}-\epsilon_{0,0}$. We see that the spectrum is mainly determined by excited states with a big overlap with the ground-state wave function $\phi_{0,0}(X)$, which only extends over a scale $\sim 1/\sqrt{\mu}$ around $X=0$. As a result, the spectral lines shown in Fig. 4(b) exhibit no particular feature at the transition point $\lambda^{2}=0.5$, which corresponds to the value of $g/\omega_{q}\simeq 0.07$ in this plot.

To identify spectral signatures of the symmetry-breaking transition in the excited state, Figs. 5(a) and (b) now show a zoom of the triplet line for $\lambda\geq\lambda_{c}$ together with the lineshapes of $S(\omega)$ for different values of $\lambda$ below and above the transition point. For $\lambda<\lambda_{c}$ we observe a small decrease of the height of the resonance peak, which can be mainly attributed to the decrease in $C$, i.e., the change in the qubit transition matrix element. However, at $\lambda=\lambda_{c}$ there is an additional sudden drop in the height of the resonance, which for $\lambda>\lambda_{c}$ also gets substantially broader. This suppression can again be understood from the overlap between $\phi_{0,0}(X)$ and the few lowest wave functions in the excited state potential, which become considerably broader at the transition point. For sufficiently small $\Gamma$, the appearance of additional sidebands below the main spectral line indicate the excitation of motional states located in one of the two displaced potential wells.

IV.3 Spectrum at finite temperature

As a second approach to observe the structural change in the excited state potential more directly, we consider in Fig. 5(c) the excitation spectrum at finite temperature. In this case, Eq. (20) must be averaged over a thermal distribution of initial states $|\phi_{0,k}\rangle$. This means that also resonator states further away from the center contribute and $S(\omega)$ probes the excited-state potential over a much larger range. For the considered temperature of $k_{B}T/\omega_{q}=0.1$, the qubits are initially still in the ground state with high probability, while a large number of resonator states $\sim k_{B}T/\omega_{r}=10$ are occupied. We see that under such conditions, the triplet line splits into two distinct branches after the transition point. As illustrated in Fig. 5(d), these two lines correspond to the energy separations between the potential curves evaluated at $X=0$ and at the minima $X=X_{\rm min}$ of the excited state potential, respectively. Since this measurement does not require any pre-cooling of the resonator, it provides a simple way to detect first signatures of the structural change of the potential curve.

V Wave packet dynamics

In this section we discuss another technique to detect changes in the potential structure through the corresponding change in the wave packet dynamics. For this purpose we consider again a situation where the ground state potential remains stable while the first excited potential curve exhibits a symmetry-breaking transition. The basic idea is to perform a Ramsey-type measurement as depicted in Fig. 6(a). Here the system is initially prepared in a superposition between the ground and the first excited qubit state such that the electromagnetic mode moves along two different potential curves simultaneously. Measurements of the qubit coherence can then be used to determine the overlap between the two wave packets as a function of time.

In Fig. 6 we illustrate this measurement protocol for the same setting as in the previous section, $N=2$ and $\varepsilon=0$. In this case, the symmetry-breaking transition occurs first in the excited potential at $\lambda=1/\sqrt{2}$. For this parameter, the ground-state potential is still approximately harmonic and active cooling methods can be applied to initialize the system in the absolute ground state $|\Psi(t=0)\rangle=|{\rm GS}\rangle\simeq|\phi_{0,0}\rangle\otimes|\chi_{0}(0)\rangle$. In a first step a microwave field of frequency $\omega_{d}$ is used to implement a $\pi/2$-rotation, which prepares the qubit state in an equal superposition between $|\chi_{0}(0)\rangle$ and the triplet state $|\chi_{T}(0)\rangle$. During this time the systems evolves according to the Hamiltonian $H(t)=H_{\rm EDM}+H_{\Omega}(t)$, where

[TABLE]

To prepare the superposition with a fidelity of about $0.95$, we set $\theta=0$, tune the frequency into resonance, $\omega_{d}=V_{T}(0)-V_{0}(0)$ and optimize the pulse time, $\tau_{\frac{\pi}{2}}\sim\pi/(2\Omega)$, for each set of parameters. Next, the system evolves freely for a waiting time $\tau_{w}$ during which the wave packet propagates along two different potential curves. In a last step, a second $\pi/2$-rotation, now with an optimized phase $\theta$, is applied, such that for $\lambda=0$ the qubits would be rotated back into the ground state. In the interacting system, variations of the return probability $P_{0}(\tau_{w})=\mathrm{Tr}\{\rho(t_{f})\ket{\chi_{0}(0)}\bra{\chi_{0}(0)}\}$, where $\rho(t_{f})$ is the density operator of the system at the final time $t_{f}=2\tau_{\frac{\pi}{2}}+\tau_{w}$, can be used to probe the wave packet dynamics during the free evolution time. This probability can be measured, for example, via regular dispersive readout schemes for qubits Blais2004 .

In Fig. 6(b) we plot the return probability $P_{0}$ as a function of $\tau_{w}$ and for varying coupling strengths. Note that in the numerical simulations the phase of the second $\pi/2$-pulse, $\theta$, has been adjusted for each parameter set to compensate trivial phase rotations due to a fixed energy offset between the ground and the excited state. The remaining variations then depend only on the wave packet dynamics and clearly distinguish the regime below and above the transition point. In the former case both potential curves are approximately quadratic and both wave packets remain localized around the origin. Above the transition the part of the wave packet promoted to the excited potential curve is expelled from the central region and undergoes oscillations. The resulting periodic decrease and revival of the wave packet overlap is clearly seen in the Ramsey fringes.

We verify that the observed modulation frequency in the Ramsey signal, which increases with larger $\lambda$, is consistent with the dynamics of a wave packet that is initialized at the center of the corresponding double-well potential. This confirms that the described Ramsey-protocol probes the wave packet dynamics in the excited state potential, more precisely, its overlap with the Gaussian ground state at $X=0$. To further illustrate this point, Fig. 6(c) shows four snapshots of the actual resonator wave function during the protocol, after projecting the system on the excited triplet state $|\chi_{T}(0)\rangle$. Initially, at $\tau_{w}=0$, the wave packet is a Gaussian centered around $X=0$. After a waiting time $\tau_{w}\sim\omega_{r}^{-1}$ most of the wave packet has propagated away from the central region, reducing the overlap with the ground state wave function. At even longer times, the reflected wave packets return, but due to the nonlinearity of the potential there is no perfect revival. This explains also the overall decay of the Ramsey fringes over a few oscillation periods. Note that in order to resolve these revivals, the qubit decoherence time $T_{2}$ must be longer than $\tau_{w}\sim\omega_{r}^{-1}$. Even for a resonator frequency as low as $\omega_{r}/(2\pi)\approx 10\,\mathrm{MHz}$, this condition can be fulfilled with realistic coherence times of $T_{2}=1-100\,\mu$s Yan2016 .

VI Implementation

In our model introduced in Sec. II we have considered the coupling of multiple flux qubits to a single resonator mode. For resonant systems, $\omega_{r}\approx\omega_{q}$, and moderate couplings such a situation can be realized by using a lumped-element $LC$ resonator. In this case the frequency of the fundamental mode, $\omega_{r}$, can be well separated from the next higher excitation mode with frequency $\omega_{\rm ex}$, such that even under USC conditions $\omega_{\rm ex}\gg\omega_{r},\omega_{q},g$ yoshihara17 . However, to observe the physics described in this work, we are interested in rather extreme ratios $\omega_{q}/\omega_{r}\gtrsim 100$, where the validity of the single-mode approximation must be evaluated in more detail.

VI.1 Two-mode circuit

Figure 7(a) shows a more realistic model for a two-qubit circuit QED system, where the self-capacitance $C_{s}$ of the inductor and the self-inductance $L_{s}$ of the capacitor are taken into account caloz2004 . In this case the resonator supports a second high-frequency mode with $\omega_{\rm ex}\sim 1/\sqrt{L_{s}C_{s}}$. In addition, we model each qubit by a superconducting loop with three Josephson junctions, as implemented in many experiments orlando1999 ; macha2014 ; schwarz2013 ; deppe2008 .

By choosing the generalized flux variables $\phi_{a}$ and $\phi_{b}$ indicated in Fig. 7(a) and the phase jumps $\Delta\phi_{1,2}$ across each qubit as dynamical degrees of freedom, the Hamiltonian for this circuit can be written as (see App. B)

[TABLE]

Here

[TABLE]

is the Hamiltonian for the $LC$ resonator and

[TABLE]

is the qubit-resonator coupling. In Eq. (23) and (24) the operators $Q_{a,b}$ and $Q_{q}$ are the conjugate charges for the flux variables $\phi_{a,b}$ and the collective qubit flux $\phi_{q}=\Delta\phi_{1}+\Delta\phi_{2}$, respectively, and $C_{q}=C_{J}(\alpha+1/2)$ and $C_{b}=C_{s}C_{q}/(C_{q}+2C_{s})$ are the relevant capacitances. Finally, $H_{q}$ is the Hamiltonian for the qubit degrees of freedom, which includes a correction term $\phi_{q}^{2}/(2L)$ from the inductive coupling. Its precise from is given in App. B, but it is not essential for following discussion.

The quadratic Hamiltonian for the electromagnetic modes can be diagonalized and written in terms of the normal mode operators $c_{\pm}$ as $H_{LC}=\sum_{\eta=\pm}\omega_{\eta}c^{\dagger}_{\eta}c_{\eta}$. In the limit of interest, $L_{s}\ll L$ and $C_{s}\ll C$, we find that $\omega_{-}\approx 1/\sqrt{LC}=\omega_{r}$ and $\omega_{+}\approx 1/\sqrt{L_{s}C_{b}}\gg\omega_{r}$. In terms of these eigenmodes, the coupling to the qubits is given by

[TABLE]

Here we have defined the dimensionless qubit variables $\varphi_{q}=\phi_{q}/\Phi_{0}$ and $n_{q}=Q_{q}/(2e)$, where $\Phi_{0}$ is the reduced flux quantum. In general, the inductive and capacitive couplings $g_{\phi,\pm}$ and $g_{Q,\pm}$ have a complicated dependence on all the system parameters (see App. B), but in the limit $\omega_{-}\ll\omega_{+}$ we obtain the approximate scaling for the low-frequency mode

[TABLE]

and

[TABLE]

for the high-frequency mode. Here $R_{Q}=h/(2e)^{2}\approx 6450\,\Omega$ is the resistance quantum. We see that for an $LC$ resonator with an impedance $Z=\sqrt{L/C}\approx 50\,\Omega$, a strong, primarily inductive coupling with $g_{\phi,-}/\omega_{-}>1$ can be achieved. Assuming $Z_{b}=\sqrt{L_{s}/C_{b}}\sim Z$, the coupling to the high frequency mode is still substantial in absolute numbers, but the relevant ratio $g_{\phi,+}/\omega_{+}$ is suppressed by a factor $L_{s}/L$.

After introducing rescaled quadratures for the low-frequency mode, as in Eq. (4), we arrive at the rescaled Hamiltonian $\tilde{H}=H/\omega_{q}$, where

[TABLE]

$\omega_{q}$ is the characteristic qubit frequency and $\mu=\omega_{q}^{2}/\omega_{-}^{2}$, as above. The Hamiltonian for the high-frequency dynamics,

[TABLE]

now also includes the high-frequency mode and depends in general also on the momentum coordinate $P$ due to the capacitive coupling. Note that here we do not make a two-level approximation such that ${H}_{q}(X,P)$ still contains the exact dynamics of all high-frequency degrees of freedom, which depend parametrically on the slowly varying coordinates $X$ and $P$.

VI.2 Discussion

In Table 1 we present a set of parameters for a realistic setup of two qubits, which result in $\omega_{q}/(2\pi)\approx 8\,\mathrm{GHz}$, $\omega_{-}/(2\pi)=50\,\mathrm{MHz}$ and an effective mass of $\mu\approx 2.5\times 10^{4}$. The coupling parameters are $g_{\phi,-}/\omega_{-}\approx 7.15$ and $g_{Q,-}/\omega_{-}\approx 0.06$. In a pure single-mode model these values translate into a coupling parameter of $\lambda^{2}\approx 0.7$, which is well above the symmetry-breaking transition in the first excited state. Note that by varying the qubit parameters by external magnetic fluxes, this coupling could also be gradually tuned across the transition forndiaz17 ; Armata2017 . For the chosen values for $L_{s}$ and $C_{s}$ we obtain $\omega_{+}/(2\pi)\approx 160$ GHz, well above the qubit energy, and $g_{Q,+}/\omega_{+}\approx 0.37$ and $g_{\phi,+}/\omega_{+}\approx 0.01$. The value of $L_{s}$ is extrapolated from the parameters reported in Ref. mckenziesell19 and could even be much lower. The assumed value of $C_{s}$ is within a factor of $\sim$2-3 of simulated and measured values for high-impedance coil inductors fink2016 ; barzanjeh17 ; barzanjeh19 . For these parameters we compare in Fig. 7 (b) the BO potentials obtained from the full two-mode Hamiltonian, Eq. (29), to the potential of the simplified one-mode EDM used in the previous sections. The potentials are qualitatively similar with only small shifts near $X=0$. Note that in our simplified model, the main limitation for $C_{s}$ arises from the renormalization of the qubit parameters, which partially is an artefact of modelling the distributed self- and stray capacitances of the real device in terms of a single capacitor. Leaving this aspect aside, much higher values of $C_{s}$ and much lower values of $\omega_{+}$ can be tolerated without substantially affecting the properties of the low-frequency mode.

As mentioned above, in the full model the adiabatic qubit energies obtained from of Eq. (29) also depends on the momentum quadrature $P$. However, since this effect arises from the much smaller charge coupling and is proportional to $S_{y}$, it only induces a small second order term $\sim(g_{Q,-}/\omega_{q})^{2}P^{2}$. This correction only results in a tiny renormalization of the effective mass, which is negligible for the current set of parameters. This analysis shows that the operation of circuit QED systems in the extreme regime of ultrastrong coupling and very low frequency is experimentally feasible. Although the current model is still an oversimplification of a real device, we expect that as long as $\omega_{\rm ex}>\omega_{q}$, the presence of higher excitation modes renormalizes the qubit parameters, but will not substantially change the qualitative features of the adiabatic potential curves and the measurement schemes discussed in this work.

VII Conclusions

In summary, we have analyzed the ultrastrong coupling of multiple superconducting flux qubits to a radio-frequency electromagnetic mode. In this regime, the dynamics of the resonator mode can be modelled as an effective particle, which moves along the adiabatic potential curves generated through the coupling to the qubits. We have shown that already for a simple two-qubit setting, the first excited potential exhibits a transition from a single- to a double-well configuration. In the limit $\omega_{r}\rightarrow 0$ such a transition can be used as a minimal instance to study properties of quantum phase transitions, as discussed in several previous works. Importantly, characteristic signatures of this transition can be probed with regular spectroscopic measurements performed only on the qubits.

These basic findings of this work clearly demonstrate that circuit QED systems in the USC regime are not limited to tests of conventional cavity QED physics. By designing more complex potentials and measurement schemes, this adiabatic approach can be used to access a large variety of particle-like physics and wave packet phenomena with a quasi-continuum of available states. This is very different from conventional superconducting quantum circuits, where one usually only has access to a few discrete modes. While in this work we have explicitly evaluated these capabilities for a low- frequency resonator mode, an alternative strategy to access this regime could be to work with $LC$ resonators of only slightly below $1\,\mathrm{GHz}$ and push the qubit frequencies to several tens of GHz. Also in this case, values of $\mu\sim 10^{3}-10^{4}$ and $\lambda\gg 1$ are feasible, but different types of measurement strategies will be required.

Acknowledgements.

We thank Georg Arnold for valuable feedback on circuit parameters. This work was supported by the Austrian Science Fund (FWF) through Grant No. P 31701-N27 and DK CoQuS, Grant No. W 1210, and by an ESQ Discovery Grant of the Austrian Academy of Sciences (ÖAW). J. J. García-Ripoll acknowledges support from AEI Project PGC2018-094792-B-I00, CSIC Research Platform PTI-001, and CAM/FEDER Project No. S2018/TCS-4342 (QUITEMAD-CM).

Appendix A Adiabatic potentials

In this appendix we provide approximate analytic results for the adiabatic potentials for certain limiting cases of interest.

A.1 Double well transition of two qubit EDM

For two qubits we can obtain the instability point of the excited state analytically through fourth-order perturbation theory. Our starting point is the qubit Hamiltonian $\tilde{H}_{q}$ given in Eq. (6). Including the bare potential $X^{2}/2$ for convenience it can be written as $\tilde{H}_{q}=\tilde{H}_{0}+\tilde{H}_{1}$, where

[TABLE]

In the case of two qubits we can diagonalize $\tilde{H}_{0}$ analytically because the triplet state $|m_{z}=0\rangle\equiv|T\rangle$, is decoupled from the other two states $|m_{z}=\pm 1\rangle\equiv|\pm 1\rangle$. The eigenenergies and eigenstates of $\tilde{H}_{0}$ are

[TABLE]

and

[TABLE]

Here we have introduced the abbreviation $\bar{\varepsilon}=1+\varepsilon$ and the normalization factor $N_{\pm}=\sqrt{1+\left(\sqrt{1+4\bar{\varepsilon}^{-2}\lambda^{-4}}\pm 2\bar{\varepsilon}^{-1}\lambda^{-2}\right)^{2}}$. The singlet state $\ket{S}$ is decoupled from the resonator and thus has the energy $\tilde{E}_{S}=X^{2}/2$, independent of $\lambda$. This leads to the splitting of the singlet and triplet states given in Eq. (14).

To obtain the instability point for the excited state $\ket{\chi_{T}}$ we calculate the perturbative corrections resulting from $\tilde{H}_{1}$ up to fourth order. The odd order perturbative corrections vanish since the semiclassical Hamiltonian is symmetric under the transformation $X\to-X,\,S_{x}\to-S_{x}$. Thus, the energy of $|\psi_{T}\rangle$ to fourth order is given by

[TABLE]

where

[TABLE]

is the second order correction. Note that other than the last term in Eq. (35), all other fourth-order correction terms vanish. From the zeroth-order results form above we find

[TABLE]

Therefore, in total we obtain

[TABLE]

such that

[TABLE]

We see that $\tilde{E}_{T}$ transforms from a single well to a double well when the quadratic term in $X$ changes sign, i.e., at

[TABLE]

The positions of the new minima above $\lambda^{2}>1/(2\sqrt{1+\varepsilon})$ are not predicted accurately by the above perturbative model since at the minimum the condition $\lambda X\ll 1$ is not respected anymore.

A.2 Structure of the ground- and first excited-state potentials for $\lambda\gg 1$

In Sec. III.3 we showed that in the limit $\lambda\gg 1$ the effective potentials are to a first approximation given by

[TABLE]

We calculate the effect of the free qubit Hamiltonian $S_{z}$ on these potentials at the minima $X=-\sqrt{2}\lambda m_{x}$ and at the points where two neighbouring minima meet $X=-(1+\varepsilon)(2m_{x}+1)\lambda/\sqrt{2}$. The second order perturbative correction to the potentials is given by

[TABLE]

which holds for $X$ away from the degeneracy points. At the minima $X=-\sqrt{2}\lambda m_{x}$ we obtain

[TABLE]

where $s_{\pm}(m_{x})=\sqrt{s(s+1)-m_{x}(m_{x}\pm 1)}$, and to obtain the simplified expression we have assumed weak interactions $|\varepsilon(2m_{x}\pm 1)|\ll 1$. Thus, for weakly interacting systems the $m_{x}=0$ state emerges as the ground state. Note that this semiclassical result agrees with the predictions from a strong-coupling perturbation theory of the full model for $\varepsilon=0$ jaako2016 .

At the degeneracy points we can diagonalize the two by two matrix describing the two spin states $m_{x}$ and $m_{x}+1$. The reduced Hamiltonian is given by

[TABLE]

where $\xi=\lambda^{2}(1+\varepsilon)\left[1+(2m_{x}+1)^{2}\varepsilon\right]$ and $\mathrm{I}_{2\times 2}$ is the $2\times 2$ identity matrix. The eigenvalues are simply

[TABLE]

Thus, the splitting between the originally degenerate states is $\tilde{\Delta}_{m_{x},m_{x}+1}=s_{+}(m_{x})$.

Appendix B Two-mode circuit QED Hamiltonian

In this appendix we present additional details about the derivation of the full Hamiltonian for the circuit depicted in Fig. 7(a). We follow the standard quantization approach vool2016 and define the generalized flux variables

[TABLE]

where $V_{a,b}(t)$ is the voltage at the nodes indicated in Fig. 7(a). The qubits are described by the phase variables $\varphi_{1,2}=\Delta\phi_{1,2}/\Phi_{0}$, which represent the phase difference across the upper junction with Josephson energy $\alpha E_{J}$. Depending on the flux qubit design, there can be additional internal dynamical degrees of freedom. For the considered three-junction design, this additional dynamical variable is the generalized flux $\phi_{-,k}$ between the two junctions in the lower arm.

The classical dynamics of this circuit is then described by the Lagrangian

[TABLE]

where the effective capacitances are $C_{q}=C_{J}(\alpha+1/2)$ and $C_{-}=C_{J}/2$ and $\phi_{q}=\sum_{k}\Delta\phi_{k}$, and $\Phi_{e}$ is an external magnetic flux threading the loop formed by the three junctions of the qubits. By introducing the canonical charges $Q_{\eta}=\partial\mathcal{L}/\partial\dot{\phi}_{\eta}$ for all flux variables and performing a Legendre transformation, we obtain the Hamiltonian given in Eq. (22) with

[TABLE]

By diagonalizing the harmonic Hamiltonian $H_{LC}$ given in Eq. (23) we obtain the mode frequencies

[TABLE]

with the bare frequencies $\omega_{a}=1/\sqrt{L_{s}C}$, $\omega_{b}=1/\sqrt{L_{b}C_{b}}$ and the coupling $g_{ab}=\sqrt{Z_{b}/Z_{a}}\omega_{a}$. Using these normal modes, we obtain the coupling Hamiltonian $H_{\rm int}$ given in Eq. (25), where the couplings are

[TABLE]

In these expressions the mixing angle $\xi$ is given by

[TABLE]

For $C_{s},L_{s}\ll C,L$, respectively, we can identify $\omega_{-}$ as a low-frequency mode with $\omega_{-}\simeq\omega_{r}$. Using the rescaled quadratures of Eq. (4) we can write the total Hamiltonian as in Eq. (28).