Hamiltonian of a flux qubit-LC oscillator circuit in the deep-strong-coupling regime

Following the standard quantization procedure, nodes are assigned to the circuit as shown in Fig. 1a. Before deriving the circuit Hamiltonian, the circuit in Fig. 1a is transformed to the one shown in Fig. 1b by applying the so-called Y-$\Delta$ transformation, by which a $\Delta$-shaped network of electrical elements is converted to an equivalent Y-shaped network or vice versa, to the inductor network. Thus, node 3 surrounded by the inductors is eliminated. The inductances of the new set of inductors are given as

and

The Lagrangian of the circuit can now be obtained relatively easily Vool and Devoret (2017):

where $C_{J}$ and $E_{J}=I_{c}\Phi_{0}/(2\pi)$ are the capacitance and the Josephson energy of the Josephson junction, $I_{c}$ is the critical current of the Josephson junction, $\Phi_{0}=h/(2e)$ is the superconducting flux quantum, and $\hat{\Phi}_{k}$ and $\dot{\hat{\Phi}}_{k}$ $(k=1,2)$ are the node flux and its time derivative for node $k$. The node charges are defined as the conjugate momenta of the node fluxes as $\hat{q}_{k}=\partial\hat{\mathcal{L}}_{circ}/\partial\dot{\hat{\Phi}}_{k}$. After the Legendre transformation, the Hamiltonian is obtained as

where

and

As can be seen from Eq. (circuit Hamiltonian), the Hamiltonian $\hat{\mathcal{H}}_{circ}$ can be separated into three parts: the first part $\hat{\mathcal{H}}_{1}$ consisting of the charge and flux operators of node 1, the second part $\hat{\mathcal{H}}_{2}$ consisting of node 2 operators, and the third part $\hat{\mathcal{H}}_{12}$ containing the product of the two flux operators.

Let us consider an alternative treatment of the circuit shown in Fig. 1a, where the circuit is assumed to be naively divided into two well-defined components. The capacitor and the inductors in the outer loop of the circuit in Fig. 1a form an LC oscillator [Fig. 1c]. The Josephson junction and the inductors in the inner loop form a single-Josephson-junction flux qubit [Fig. 1d]. The LC oscillator and the single-Josephson-junction flux qubit share the inductor $L_{c}$ at the common part of the two loops. It is now instructive to investigate the relation of the following pairs of Hamiltonians: $\hat{\mathcal{H}}_{1}$ and the Hamiltonian of the LC oscillator shown in Fig. 1c, $\hat{\mathcal{H}}_{2}$ and the Hamiltonian of the flux qubit shown in Fig. 1d, and $\hat{\mathcal{H}}_{12}$ and the Hamiltonian of the inductive coupling between the LC oscillator and the flux qubit $M\hat{I}_{LC}\hat{I}_{q}=-L_{c}\hat{\Phi}_{1}\hat{\Phi}_{2}/[(L_{c}+L_{1})(L_{c}+L_{2})]$, where we have used the relations of the oscillator current $\hat{I}_{LC}=\hat{\Phi}_{1}/(L_{c}+L_{1})$, the qubit current $\hat{I}_{q}=\hat{\Phi}_{2}/(L_{c}+L_{2})$, and the mutual inductance $M=L_{c}$. Actually, only the inductances are different in the Hamiltonians of the two pictures: The inductances of the LC oscillator, the flux qubit, and the coupling Hamiltonians derived from the separate treatment are respectively $L_{c}+L_{1}$, $L_{c}+L_{2}$, and $-(L_{c}+L_{1})(L_{c}+L_{2})/L_{c}$, while those in $\hat{\mathcal{H}}_{circ}$ are $L_{LC}$, $L_{FQ}$, and $-L_{12}$. Figure 2 shows the inductances in the separate treatment and in $\hat{\mathcal{H}}_{circ}$ as functions of $L_{c}$ on condition that the inductance sums are kept constant at $L_{c}+L_{1}=800$ pH and $L_{c}+L_{2}=2050$ pH. As $L_{c}$ approaches 0, or more specifically when $L_{c}\ll L_{1},L_{2}$, we obtain $L_{c}+L_{1}$ $\sim$ $L_{LC}$, $L_{c}+L_{2}$ $\sim$ $L_{FQ}$, and $(L_{c}+L_{1})(L_{c}+L_{2})/L_{c}$ $\sim$ $L_{12}$. In this way, the Hamiltonian derived from the separate treatment has the same form as $\hat{\mathcal{H}}_{circ}$, and the inductances in the separate treatment approach those of $\hat{\mathcal{H}}_{circ}$.

We compare the Hamiltonian $\hat{\mathcal{H}}_{circ}$ with the generalized quantum Rabi Hamiltonian:

The first part $\hat{\mathcal{H}}_{LC}$ represents the energy of the LC oscillator, where $\hat{a}^{\dagger}$ and $\hat{a}$ are the creation and annihilation operators, respectively. The second part $\hat{\mathcal{H}}_{FQ}$ represents the energy of the flux qubit written in the energy eigenbasis at $\varepsilon=0$. The operators $\hat{\sigma}_{x,z}$ are the standard Pauli operators. The parameters $\hbar\Delta_{q}$ and $\hbar\varepsilon$ are the tunnel splitting and the energy bias between the two states with persistent currents flowing in opposite directions around the qubit loop. The third part $\hat{\mathcal{H}}_{coup}$ represents the inductive coupling energy.

The relation between $\hat{\mathcal{H}}_{1}$ and $\hat{\mathcal{H}}_{LC}$ is straightforward. The resonance frequency and the operators in $\hat{\mathcal{H}}_{LC}$ can be analytically described by the variables and operators in $\hat{\mathcal{H}}_{1}$ as $\omega=1/\sqrt{L_{LC}C}$ and $\hat{a}+\hat{a}^{\dagger}\to\hat{\Phi}_{1}/(L_{LC}I_{zpf})$, where $I_{zpf}=\sqrt{\hbar\omega/(2L_{LC})}$ is the zero-point-fluctuation current.

To see the relation between $\hat{\mathcal{H}}_{2}$ and $\hat{\mathcal{H}}_{FQ}$, we numerically calculated the eigenenergies of $\hat{\mathcal{H}}_{2}$ as functions of $\Phi_{x}$. In the calculation, we used the following parameters: $L_{c}+L_{2}=2050$ pH, $L_{J}=990$ pH ($E_{J}/h=165.1$ GHz), and $C_{J}=4.84$ fF ($E_{C}/h=4.0$ GHz). As shown in Fig. 3a the lowest two energy levels are well separated from the higher levels, which are more than 40 GHz higher in frequency. The lowest two energy levels of $\hat{\mathcal{H}}_{2}$ are well approximated by $\hat{\mathcal{H}}_{FQ}$ [Fig. 3b], which gives almost identical results obtained by the local basis reduction method Consani and Warburton (2020), with $\hat{\sigma}_{x}$ and $\hat{\sigma}_{z}$ swapped. The eigen frequencies of the ground state $\omega_{0}$ and the first excited state $\omega_{1}$ are respectively fitted by $\omega_{os}-\sqrt{\varepsilon^{2}+\Delta_{q}^{2}}/2$ and $\omega_{os}+\sqrt{\varepsilon^{2}+\Delta_{q}^{2}}/2$. Besides the offset $\omega_{os}$, the fitting parameters are $\Delta_{q}$ and the maximum persistent current $I_{p}$, which is determined as the proportionality constant between the energy bias and the flux bias, $\varepsilon=2I_{p}(\Phi_{x}-0.5\Phi_{0})$. We also numerically calculated the expectation values of the flux operator, $\left\langle g\left|\hat{\Phi}_{2}\right|g\right\rangle$ and $\left\langle e\left|\hat{\Phi}_{2}\right|e\right\rangle$, which are well approximated by $\left\langle g\left|\hat{\sigma}_{x}\right|g\right\rangle$ and $\left\langle e\left|\hat{\sigma}_{x}\right|e\right\rangle$, respectively [Fig. 3c]. The expectation values of the flux operator $\left\langle g\left|\hat{\Phi}_{2}\right|g\right\rangle$ and $\left\langle e\left|\hat{\Phi}_{2}\right|e\right\rangle$ are respectively fitted by $-\Phi_{2max}\varepsilon/\sqrt{\varepsilon^{2}+\Delta_{q}^{2}}$ and $\Phi_{2max}\varepsilon/\sqrt{\varepsilon^{2}+\Delta_{q}^{2}}$. Here, $\left|g\right\rangle$ and $\left|e\right\rangle$ are respectively the ground and excited states of the qubit Hamiltonian $\hat{\mathcal{H}}_{2}$. The only fitting parameter is determined by the ratio $\Phi_{2max}=-\left\langle i\left|\hat{\Phi}_{2}\right|i\right\rangle/\left\langle i\left|\hat{\sigma}_{x}\right|i\right\rangle$ ($i=g,e$). The Pauli operator $\hat{\sigma}_{x}$ is therefore identified as being proportional to the flux operator $\hat{\sigma}_{x}\to-\hat{\Phi}_{2}/\Phi_{2max}$. The relation between $\hat{\mathcal{H}}_{12}$ and $\hat{\mathcal{H}}_{coup}$ can now be obtained by using the relations for the oscillator and qubit operators identified above. This way we find that the Hamiltonian $\hat{\mathcal{H}}_{12}$ can be expressed as $-(L_{LC}/L_{12})I_{zpf}\Phi_{2max}\hat{\sigma}_{x}(\hat{a}+\hat{a}^{\dagger})$, which is exactly the same form as $\hat{\mathcal{H}}_{coup}$, with the coupling strength $\hbar g=-(L_{LC}/L_{12})I_{zpf}\Phi_{2max}$. Note that $\mathcal{H}_{coup}$ directly derived from the Lagrangian $\hat{\mathcal{L}}_{circ}$ is in the flux gauge, as the qubit-oscillator coupling term is of the form $\hat{\Phi}_{1}\hat{\Phi}_{2}$, which is optimal for our system with a single oscillator mode Stokes and Nazir (2019); Roth et al. (2019); Di Stefano et al. (2019). Excluding the qubit’s energy levels higher than the first excited states, $\hat{\mathcal{H}}_{circ}$ takes the form of $\hat{\mathcal{H}}_{R}$. In other words, once the circuit parameters, i.e. $\Phi_{x}$, $L_{c}$, $L_{1}$, $L_{2}$, $C$, $E_{J}$, and $C_{J}$, are given, the corresponding parameters in $\hat{\mathcal{H}}_{R}$ ($\omega$, $\varepsilon$, $\Delta_{q}$, and $g$) are set. The relation between $\hat{\mathcal{H}}_{circ}$ and $\hat{\mathcal{H}}_{R}$ is summarized in Table 1.

As previously mentioned, $\hat{\mathcal{H}}_{R}$ considers only the lowest two energy levels of the flux qubit, while $\hat{\mathcal{H}}_{circ}$ includes all energy levels. To investigate the effect of the qubit’s higher energy levels, we perform numerical calculations and compare the energy levels calculated by $\hat{\mathcal{H}}_{R}$ and $\hat{\mathcal{H}}_{circ}$. The details of the numerical diagonalization of $\hat{\mathcal{H}}_{circ}$ are given in methods. Since the contributions from the qubit’s higher energy levels are expected to become larger as the coupling strength increases, we consider parameters that cover a wide range of coupling strengths from the weak-coupling to the deep-strong-coupling regime. In the calculations, we fix $L_{c}+L_{1}=800$ pH, $L_{c}+L_{2}=2050$ pH, $C=0.87$ pF, $L_{J}=990$ pH ($E_{J}/h=165.1$ GHz), and $C_{J}=4.84$ fF ($E_{C}/h=4.0$ GHz), and sweep the flux bias $\Phi_{x}$ around $\Phi_{0}/2$ at various values of $L_{c}$. These parameters are used in all the calculations for Figs. 4-7. Some of our calculations were performed using the QuTiP simulation package Johansson et al. (2013).

Transition frequencies of the qubit-oscillator circuit $\omega_{ij}$ corresponding to the transition $\left|i\right\rangle\to\left|j\right\rangle$ numerically calculated from $\hat{\mathcal{H}}_{circ}$ around the resonance frequency of the oscillator $\omega$ are plotted in Fig. 4 for a $L_{c}=20$ pH and b $L_{c}=350$ pH, where the indices $i$ and $j$ label the energy eigenstates according to their order in the energy-level ladder, with the index 0 denoting the ground state. The same circuit parameters as in Fig. 3 are used. In the case $L_{c}=20$ pH, two characteristic features are observed: avoided level crossings between the qubit and oscillator transition signals approximately at $\Phi_{x}/\Phi_{0}=0.497$ and 0.503, and the dispersive shift that for example creates the separation between the frequencies of the transitions $\left|0\right\rangle\to\left|2\right\rangle$ and $\left|1\right\rangle\to\left|3\right\rangle$, leading to the peak/dip feature at the symmetry point, i.e. $\Phi_{x}/\Phi_{0}=0.5$. Note that transitions $\left|0\right\rangle\to\left|2\right\rangle$ and $\left|1\right\rangle\to\left|3\right\rangle$ around $\Phi_{x}/\Phi_{0}=0.5$ respectively correspond to transitions $\left|g0\right\rangle\to\left|g1\right\rangle$ and $\left|e0\right\rangle\to\left|e1\right\rangle$, where “$g$” and “$e$” denote, respectively, the ground and excited states of the qubit, and “0” and “1” the number of photons in the oscillator’s Fock state. In the case $L_{c}=350$ pH, the characteristic spectrum indicates that the qubit-oscillator circuit is in the deep-strong-coupling regime Yoshihara et al. (2017b). Transition frequencies of $\hat{\mathcal{H}}_{R}$ are also plotted in Figs. 4a and b. It should be mentioned that the parameters of $\hat{\mathcal{H}}_{R}$ are obtained in two different ways. Here, the parameters of $\hat{\mathcal{H}}_{R}$ are obtained from the relations described in Table 1. The overall shapes of the spectra of $\hat{\mathcal{H}}_{R}$ and $\hat{\mathcal{H}}_{circ}$ look similar. On the other hand, the shift of the entire spectrum becomes as large as more than 200 MHz for $L_{c}=350$ pH. To quantify the difference of the spectra between $\hat{\mathcal{H}}_{R}$ and $\hat{\mathcal{H}}_{circ}$, transition frequencies up to the third excited state numerically calculated from $\hat{\mathcal{H}}_{circ}$ are fitted by $\hat{\mathcal{H}}_{R}$. In the fitting, $\hat{\mathcal{H}}_{R}$ with an initial approximate set of parameters is numerically diagonalized and then the parameters, $\omega$, $\Delta_{q}$, $g$, and $I_{p}(=\varepsilon/[2(\Phi_{x}-0.5\Phi_{0})])$, are varied to obtain the best fit. Figures 4c and d show that the same spectra of $\hat{\mathcal{H}}_{circ}$ in Figs. 4a and b are well fitted by $\hat{\mathcal{H}}_{R}$, but with a different parameter set.

Figure 5 shows parameters of $\hat{\mathcal{H}}_{R}$ obtained from the relations described in Table 1, and by fitting the spectra of $\hat{\mathcal{H}}_{circ}$ up to the third excited state to $\hat{\mathcal{H}}_{R}$ as a function of $L_{c}$. The parameters of $\hat{\mathcal{H}}_{R}$ obtained in these two different ways become quite different from each other at large values of $L_{c}$ due to contributions of the qubit’s higher energy levels. Interestingly, the average of the squares of the residuals in the least-squares method for obtaining the parameters of $\hat{\mathcal{H}}_{R}$ by fitting, $\overline{[\delta\omega_{0i}/(2\pi)]^{2}}$ ($i$ = 1, 2, and 3), remains rather small, at most $25~{}$MHz², which is consistent with the good fitting shown in Figs. 4c and d. In fact, the energy level structure of $\hat{\mathcal{H}}_{circ}$ up to the seventh excited state can still be fitted well by the quantum Rabi Hamiltonian as shown in Supplementary Information SI .

The circuit Hamiltonian $\hat{\mathcal{H}}_{circ}$ is in the flux gauge, as the qubit-oscillator coupling term is of the form $\hat{\Phi}_{1}\hat{\Phi}_{2}$. The Hamiltonian can be transformed into the charge gauge:

where

The details of the calculations are given in Supplementary Information SI . As can be seen in Eq. (Hamiltonian in the charge gauge), the Hamiltonian $\hat{\mathcal{H}}_{circ}^{\prime}$ can be separated into three parts: the first part $\hat{\mathcal{H}}_{1}^{\prime}$ consisting of the charge and flux operators of node 1, the second part $\hat{\mathcal{H}}_{2}^{\prime}$ consisting of the charge and flux operators of node 2, and the third part $\hat{\mathcal{H}}_{12}^{\prime}$ containing the product of the two charge operators. It is worth mentioning that the inductance in $\hat{\mathcal{H}}_{2}^{\prime}$ is equal to $L_{c}+L_{2}$, which is the inductance of the flux qubit in the separate treatment:

The overall form of $\hat{\mathcal{H}}^{\prime}_{circ}$ [Eq. (Hamiltonian in the charge gauge)] is almost the same as that of $\hat{\mathcal{H}}_{circ}$ [Eq. (circuit Hamiltonian)], with the exception that the coupling term is of the form $\hat{q}_{1}\hat{q}_{2}$ instead of $\hat{\Phi}_{1}\hat{\Phi}_{2}$. When we examined the mapping between $\hat{\mathcal{H}}_{2}$ and $\hat{\mathcal{H}}_{FQ}$, we explained that the operator $\hat{\Phi}_{2}$ can be identified as $\hat{\sigma}_{x}$. Similarly, if we calculate the matrix elements of the operator $\hat{q}_{2}$ for the two lowest qubit states, we find that $\left\langle g\left|\hat{q}_{2}\right|g\right\rangle=\left\langle e\left|\hat{q}_{2}\right|e\right\rangle=0$ and $\left\langle g\left|\hat{q}_{2}\right|e\right\rangle=-\left\langle e\left|\hat{q}_{2}\right|g\right\rangle\neq 0$ [Fig. 6a]. We can therefore identify the operator $\hat{q}_{2}$ as the qubit operator $\hat{\sigma}_{y}$.

It then seems natural to look for a mapping of $\hat{\mathcal{H}}^{\prime}_{circ}$ to the variant of the quantum Rabi Hamiltonian:

Compared with $\hat{\mathcal{H}}_{R}$, only $\hat{\mathcal{H}}_{coup}$ is replaced by $\hat{\mathcal{H}}_{coup}^{\prime}$. It should be noted, however, that there is a physical difference between the two Hamiltonians and that there is no simple mapping between them. The operator $\hat{\sigma}_{x}$ in $\hat{\mathcal{H}}_{coup}$ is the same as one of the Pauli operators in $\hat{\mathcal{H}}_{FQ}$. In contrast, the operator $\hat{\sigma}_{y}$ in $\hat{\mathcal{H}}_{coup}^{\prime}$ is different from the two operators ($\hat{\sigma}_{x}$ and $\hat{\sigma}_{z}$) in $\hat{\mathcal{H}}_{FQ}$. As a result, $\hat{\mathcal{H}}_{R}$ and $\hat{\mathcal{H}}^{\prime}_{R}$ are physically different and for example produce different spectra.

In Fig. 4 we plot a few of the transition frequencies in the circuit’s spectrum along with the corresponding transition frequencies obtained from $\hat{\mathcal{H}}_{R}^{\prime}$. When the parameters of $\hat{\mathcal{H}}_{R}^{\prime}$ are read off $\hat{\mathcal{H}}_{circ}^{\prime}$, similarly to what is shown in Table I, Figs. 4e and 4f show that the fitting is poor in most parts of the spectrum, even for the relatively weak coupling case $L_{c}=20$ pH. Here, $g^{\prime}=q_{1zpf}q_{2max}L_{LC}/(C_{J}L_{12})$, $q_{1zpf}=\sqrt{\hbar\omega^{\prime}C^{\prime}/2}$, $\omega^{\prime}=1/\sqrt{L_{LC}}\times\sqrt{(1/C)+(L_{LC}^{2}/C_{J}L_{12}^{2})}$, $1/C^{\prime}=(1/C)+[L_{LC}^{2}/(C_{J}L^{2}_{12})]$, and $q_{2max}$ is numerically calculated as shown in Fig. 6a. Contrary to the open grey and red circles obtained from $\hat{\mathcal{H}}_{circ}^{\prime}$, the solid grey and red lines obtained from $\hat{\mathcal{H}}_{R}^{\prime}$ are larger and do not show the peaks and dips around the symmetry point. Even with a numerical optimization of the fitting parameters [Figs. 4g and 4h], only parts of the spectrum can be fitted well. In particular, the peaks and dips that occur in the spectrum at $\Phi_{x}/\Phi_{0}=0.5$ are not reproduced by $\hat{\mathcal{H}}_{R}^{\prime}$.

Here it is useful to consider the two characteristic features in the spectrum, i.e. the avoided crossings and dispersive shifts, as being the result of energy shifts in the spectrum of an uncoupled circuit when the coupling term is added. The gap of an avoided level crossing, or in other words the Rabi splitting, is proportional to the matrix element of the coupling term between the relevant energy eigenstates of the uncoupled system. The dispersive shift of one energy level caused by another energy level is proportional to the square of the matrix element between the two energy eigenstates according to perturbation theory. The details of the dispersive shifts up to second order in perturbation theory are described in Supplementary information SI .

The difference between the spectra shown in the solid lines in Figs. 4a and 4e is attributed to the difference between the matrix elements of the qubit’s flux and charge operators, $\hat{\Phi}_{2}$ and $\hat{q}_{2}$. The flux bias dependences of the numerically calculated matrix elements of the qubit’s charge and flux operators are shown in Fig. 6. Regarding the matrix elements of the qubit’s flux operator $\left|\left\langle j\left|\hat{\Phi}_{2}\right|i\right\rangle\right|$ ($i=g,e$), those involving the higher qubit levels $j=f,h,k,l$ are smaller than those of $j=g,e$. Here, $\left|f\right\rangle$, $\left|h\right\rangle$, $\left|k\right\rangle$, and $\left|l\right\rangle$ respectively represent the second, third, fourth, and fifth excited states of $\hat{\mathcal{H}}_{2}$. Regarding the matrix elements of the qubit’s charge operator $\left|\left\langle j\left|\hat{q}_{2}\right|i\right\rangle\right|$ ($i=g,e$), on the other hand, some of those involving the higher qubit levels $j=f,h,k,l$ are significantly larger than those of $j=g,e$ in most of the flux bias range. This difference stems from the fact that the qubit Hamiltonian involves a double-well potential of the flux variables, as explained in Supplementary information SI .

The matrix elements $\left\langle g\left|\hat{\Phi}_{2}\right|e\right\rangle$ and $\left\langle e\left|\hat{\Phi}_{2}\right|g\right\rangle$ have a peak at $\Phi/\Phi_{0}=0.5$, which directly leads to the peak/dip feature at the symmetry point as shown in the solid lines in Fig. 4a. On the other hand, the matrix elements $\left\langle g\left|\hat{q}_{2}\right|e\right\rangle$ and $\left\langle e\left|\hat{q}_{2}\right|g\right\rangle$ are almost constant in the flux bias range of Fig. 6, which is consistent with the absence of a peak/dip feature at the symmetry point in the solid lines in Fig. 4e. Instead, the matrix elements $\left\langle h\left|\hat{q}_{2}\right|g\right\rangle$ and $\left\langle f\left|\hat{q}_{2}\right|e\right\rangle$, which have somewhat similar flux-bias dependence to those of $\left\langle g\left|\hat{\Phi}_{2}\right|e\right\rangle$ and $\left\langle e\left|\hat{\Phi}_{2}\right|g\right\rangle$, have peaks at $\Phi/\Phi_{0}=0.5$. As a result, for the flux gauge, the contribution from higher levels is just a small correction, while it cannot be neglected in the charge gauge. For this reason, the flux gauge turns out to be more convenient for purposes of mapping $\hat{\mathcal{H}}_{circ}$ to the quantum Rabi Hamiltonian.

It is worth mentioning that we are dealing with an inductively coupled circuit, and one might think that the nature of the coupling, i.e. inductive or capacitive, will determine which gauge is more suitable, especially because the two gauges differ mainly by the form of the coupling term. We show that the deciding factor is the qubit rather than the coupling. As a result, if we have a flux qubit capacitively coupled to the oscillator, the flux gauge will be more suitable for the purpose of approximating the circuit Hamiltonian by the quantum Rabi Hamiltonian. It is also worth mentioning that our results should apply to fluxonium-resonator circuits, since the fluxonium Hamiltonian close to the degeneracy point also involves a double-well potential of the flux variables. It should be noted, however, that the fluxonium’s transition frequencies to the higher energy levels, $\omega_{02}$ and $\omega_{03}$, are smaller than those of flux qubits, and, hence, the contribution from higher levels would be larger in fluxonium-resonator circuits.

We note here that the differences between the energy spectra obtained from $\hat{\mathcal{H}}_{circ}$, $\hat{\mathcal{H}}_{R}$, and $\hat{\mathcal{H}}_{R}^{\prime}$ in the case $\varepsilon=0$, $\Delta_{q}=\omega$, and a quadratic-plus-quartic potential energy function were extensively studied in Ref. De Bernardis et al. (2018). Contrary to the case of $\Delta_{q}<\omega$, the spectra obtained from $\hat{\mathcal{H}}_{circ}$ and $\hat{\mathcal{H}}_{R}$ are almost identical for the lowest energy levels in the whole range of $g/\omega$ while that of $\hat{\mathcal{H}}_{R}^{\prime}$ deviates from the other two when $g/\omega\gtrsim 0.1$. The Rabi splitting between the first two excited energy levels, which is symmetric and is proportional to $g$, is clearly visible for $g/\omega\ll 0.1$, while the dispersive shift due to the higher energy levels is proportional to $g^{2}$, and can be observed only when $g/\omega\gtrsim 0.1$. We also note that for $L_{c}=20$ pH, $g^{\prime}/(2\pi)=0.043$ GHz, which is much smaller than $g/(2\pi)=0.424$ GHz. On the other hand, the gap of the avoided level crossings plotted as lines in Fig. 4e is not much smaller than that in Fig. 4a considering $g^{\prime}\sim 0.1g$. The gap of the avoided level crossings for $\hat{\mathcal{H}}_{R}^{\prime}$ is $2g^{\prime}$ while that for $\hat{\mathcal{H}}_{R}$ is $2g\Delta_{q}/\omega$. The factor $\Delta_{q}/\omega=0.206$ partly explains why the avoided level crossings for $\hat{\mathcal{H}}_{R}$ and $\hat{\mathcal{H}}_{R}^{\prime}$ are not very different from each other.

One of the most paradoxical features of $\hat{\mathcal{H}}_{R}$ in the deep-strong-coupling regime is the non-negligible number of photons in the oscillator in the ground state Ashhab and Nori (2010). In terms of the creation and annihilation operators, the photon number operator is $\hat{a}^{\dagger}\hat{a}$. Considering the mapping between $\hat{\mathcal{H}}_{circ}$ and $\hat{\mathcal{H}}_{R}$, the photon number operator in $\hat{\mathcal{H}}_{circ}$ is $\hat{\mathcal{H}}_{1}/(1/\sqrt{L_{LC}C})-0.5$. As shown in Fig. 7a, a nonzero number of photons in the oscillator is obtained. Another paradoxical feature of $\hat{\mathcal{H}}_{R}$ in the deep-strong-coupling regime is a non-negligible expectation value of the field operator $\left\langle\hat{a}+\hat{a}^{\dagger}\right\rangle$ when the qubit flux bias is $\Phi_{x}\neq 0.5\Phi_{0}$. The corresponding operator in $\hat{\mathcal{H}}_{circ}$ is $\hat{\Phi}_{1}/(L_{LC}I_{zpf})$. Fig. 7b shows the expectation values of the fluxes, $\left\langle\hat{\Phi}_{1}\right\rangle$ and $\left\langle\hat{\Phi}_{2}\right\rangle$, as functions of the flux bias $\Phi_{x}$ for the ground and the first excited states in the case $L_{c}=350$ pH. At $\Phi_{x}/\Phi_{0}\neq 0.5$, a nonzero expectation value $\left\langle\hat{\Phi}_{1}\right\rangle$ is demonstrated. One may suspect that a nonzero $\left\langle\hat{\Phi}_{1}\right\rangle$ is unphysical because it might result in a nonzero DC current through the inductor of an LC oscillator in an energy eigenstate. The relation between the flux and current operators is explicitly given in the circuit model [Fig. 1a]:

where $\hat{I}_{k}$ ($k=1,2$) is defined as the operator of the current flowing from the ground node to node $k$. In the case $L_{c}=0$, Eq. (19) reduces to $\hat{I}_{1}=\hat{\Phi}_{1}/L_{1}$ and $\hat{I}_{2}=\hat{\Phi}_{2}/L_{2}$ which means that the operator $\hat{\Phi}_{1}$ can indeed be understood as a current operator. As shown in Fig. 7c, at $\Phi_{x}/\Phi_{0}=0.498$, the expectation value $2\pi\left\langle\hat{\Phi}_{1}\right\rangle/\Phi_{0}$ in the ground state is proportional to $L_{c}$ and is zero only when $L_{c}=0$. Fig. 7d shows the expectation values of the currents $\left\langle\hat{I}_{1}\right\rangle$ and $\left\langle\hat{I}_{2}\right\rangle$ as functions of $\Phi_{x}$ in the case $L_{c}=350$ pH. The expectation value $\left\langle\hat{I}_{2}\right\rangle$ is nonzero at $\Phi_{x}/\Phi_{0}\neq 0.5$, while $\left\langle\hat{I}_{1}\right\rangle$ is exactly zero at all values of $\Phi_{x}$. In this way, although $\left\langle\hat{\Phi}_{1}\right\rangle$ is nonzero in the case $L_{c}\neq 0$, an unphysical DC current through the inductor of an LC oscillator is not predicted.

Next, we consider the case of $\hat{\mathcal{H}}_{circ}^{\prime}$. Since unitary transformations do not change the eigenenergies of Hamiltonians, $\hat{\mathcal{H}}_{circ}$ and $\hat{\mathcal{H}}_{circ}^{\prime}$ have exactly the same eigenenergies. On the other hand, the photon-number and flux operators do not commute with the gauge transformation:

and

The corresponding photon number operator in $\hat{\mathcal{H}}_{circ}^{\prime}$ is $\hat{\mathcal{H}}_{1}^{\prime}/(\hbar\omega^{\prime})-0.5$, where $\omega^{\prime}=1/\sqrt{L_{LC}}\times\sqrt{(1/C)+(L_{LC}^{2}/C_{J}L_{12}^{2})}$ is the resonance frequency of $\hat{\mathcal{H}}_{1}^{\prime}$. The corresponding field operator in $\hat{\mathcal{H}}_{circ}^{\prime}$ is $\hat{\Phi}_{1}/(L_{LC}I_{zpf})$. As shown in Fig. 7e, the expectation number of photons in the charge gauge is much smaller than that in the flux gauge, and the non-negligible expectation number of photons in the oscillator in the ground state arises only in the flux-gauge. We note here that the resonance frequency of $\hat{\mathcal{H}}_{1}^{\prime}$ is larger than $1/\sqrt{L_{LC}C}\simeq 6.03~{}$GHz: $\omega^{\prime}/(2\pi)=15.66$ GHz in the case $L_{c}=350$ pH. We can see that the expectation value of the operator $\hat{\Phi}_{1}$ in the case $L_{c}=350$ pH is zero at all values of $\Phi_{x}$ [Fig. 7f], which are also different from the case of the flux-gauge Hamiltonian. The expectation value of the current operator $\hat{I}_{1}$ is zero at all values of $\Phi_{x}$, which is true in the flux gauge as well.