Superconducting qubits beyond the dispersive regime

TL;DR

This paper introduces a new formalism for analyzing superconducting circuits beyond the dispersive regime, enabling accurate modeling of qubit interactions as circuits scale up and detuning decreases.

Contribution

The authors develop a formalism that diagonalizes superconducting circuit Hamiltonians beyond the dispersive regime, allowing for accurate analysis of qubit interactions at small or negligible detuning.

Findings

Formalism reproduces perturbative results in the dispersive regime.

Formalism extrapolates to regimes with small or zero detuning.

Closed-form formulas for circuit characteristics like dressed frequencies and Kerr couplings.

Abstract

Superconducting circuits consisting of a few low-anharmonic transmons coupled to readout and bus resonators can perform basic quantum computations. Since the number of qubits in such circuits is limited to not more than a few tens, the qubits can be designed to operate within the dispersive regime, where frequency detuning are much stronger than coupling strengths. However, scaling up the number of qubits will bring the circuit out of this regime and invalidates current theories. We develop a formalism that allows to consistently diagonalize superconducting circuit hamiltonian beyond dispersive regime. This will allow to study qubit-qubit interaction unperturbatively, therefore our formalism remains valid and accurate at small or even negligible frequency detuning; thus our formalism serves as a theoretical ground for designing qubit characteristics for scaling up the number of qubits in superconducting circuits. We study the most important circuits with single- and two-qubit gates, i.e. a single transmon coupled to a resonator and two transmons sharing a bus resonator. Surprisingly our formalism allows to determine the circuit characteristics, such as dressed frequencies and Kerr couplings, in closed-form formulas that not only reproduce perturbative results but also extrapolate beyond the dispersive regime and can ultimately reproduce (and even modify) the Jaynes-Cumming results at resonant frequencies.

Equations

$$H=H_{\text{har}}+H_{\text{anhar}},H_{\text{anhar}}=-\frac{E_C}{3Z_t^2{\hslash }^2}{\varphi }_t^4$$

$$H_{\text{har}}=\sum _{i=r,t}\frac{q_i^2}{2C_i}+\frac{{\varphi }_i^2}{2L_i}+H_{int}$$

$$H_{\textup{har}}=\frac{1}{2}\mathbf{Q}^{\textup{T}}\mathbf{MQ}+\frac{1}{2}\mathbf{X}^{\textup{T}}\mathbf{X};\ \mathbf{M}=\left[\begin{array}[]{ccccc}\omega_{t}^{2}&g\sqrt{4\omega_{t}\omega_{r}}\\ g\sqrt{4\omega_{t}\omega_{r}}&\omega_{r}^{2}\end{array}\right],$$

$$H$$

$${\tilde{\omega }}_t=K_{-}^2-\frac{{\chi }_t}{2}-{\chi }_{rt},{\tilde{\omega }}_r=K_{+}^2-\frac{{\chi }_r}{2}-{\chi }_{rt}$$

$$\displaystyle{\chi}_{t}=\delta{\left(1+s\right)^{2}}K_{-}^{4}/{4\omega_{t}^{2}},\ \ {\chi}_{r}=\delta{\left(1-s\right)^{2}}{K_{+}^{4}}/{4\omega_{t}^{2}}\ \ \ \ \$$

$${\hat{a}}_i^{†}-{\hat{a}}_i$$

$$\mathbf{M}=\left[\begin{array}[]{ccc}\omega_{1}^{2}&0&g_{1}\sqrt{4\omega_{1}\omega_{3}}\\ 0&\omega_{2}^{2}&g_{2}\sqrt{4\omega_{2}\omega_{3}}\\ g_{1}\sqrt{4\omega_{1}\omega_{3}}&g_{2}\sqrt{4\omega_{2}\omega_{3}}&\omega_{3}^{2}\end{array}\right]$$

$${\bar{\omega }}_k^2=2\sqrt{\frac{f}{3}}\cos \frac{{\text{cos}}^{-1}\left(-\frac{h}{2}{\left(\frac{3}{f}\right)}^{\frac{3}{2}}\right)-2\pi (k-1)}{3}-\frac{b}{3}$$

$$H=\sum _{i=1,2,3}{\tilde{\omega }}_i{\hat{\alpha }}_i^{†}{\hat{\alpha }}_i$$

$$\displaystyle-\sum_{i=1,2,3}\Bigg{\{}\frac{\chi_{i}}{2}\left(\hat{\alpha}_{i}^{\dagger}\hat{\alpha}_{i}\right)^{2}+2\sum_{k>i}\chi_{ik}\left(\hat{\alpha}_{i}^{\dagger}\hat{\alpha}_{i}\right)\left(\hat{\alpha}_{k}^{\dagger}\hat{\alpha}_{k}\right)\Bigg{.}$$

$$+\sum _{k>i}\left({\mathcal{J}}_{ik}+\sum _{l\ne i,k}{\mathcal{L}}_{ikl}{\hat{\alpha }}_l^{†}{\hat{\alpha }}_l\right)\left(\widehat{{\alpha }_i}{\hat{\alpha }}_k^{†}+{\hat{\alpha }}_i^{†}{\hat{\alpha }}_k\right)$$

$$\displaystyle\Bigg{.}\ \ \ +\sum_{k\neq i}\mathcal{K}_{ik}\left[\left(\hat{\alpha}_{i}^{\dagger}\hat{\alpha}_{i}\right)\hat{\alpha}_{i}\hat{\alpha}_{k}^{\dagger}+\hat{\alpha}_{i}^{\dagger}\hat{\alpha}_{k}\left(\hat{\alpha}_{i}^{\dagger}\hat{\alpha}_{i}\right)\right]\Bigg{\}}$$

$${\tilde{\omega }}_i={\bar{\omega }}_i-\frac{{\chi }_i}{2}-\sum _{j(\ne i)}{\chi }_{ij}.$$

$$2\sqrt{\frac{f}{3}}\cos \frac{{\text{cos}}^{-1}\left(-\frac{h}{2}{\left(\frac{3}{f}\right)}^{\frac{3}{2}}\right)-2\pi (k-1)}{3}-\frac{b}{3}\ge 0$$

$${\omega }_1^4+{\omega }_2^4+{\omega }_3^4\ge {\omega }_1^2{\omega }_2^2+{\omega }_1^2{\omega }_3^2+{\omega }_2^2{\omega }_3^2-4g_1^2{\omega }_1{\omega }_2-4g_2^2{\omega }_2{\omega }_3$$

$$({\omega }_1^2-{\omega }_2^2{)}^2+({\omega }_2^2-{\omega }_3^2{)}^2+({\omega }_3^2-{\omega }_1^2{)}^2\ge -4g_1^2{\omega }_1{\omega }_2-4g_2^2{\omega }_2{\omega }_3$$

$$S=\left[\begin{array}[]{ccc}\frac{V_{1}\gamma_{12}}{N_{1}}&\frac{V_{1}\gamma_{22}}{N_{2}}&\frac{V_{1}\gamma_{32}}{N_{3}}\\ \frac{V_{2}\gamma_{11}}{N_{1}}&\frac{V_{2}\gamma_{21}}{N_{2}}&\frac{V_{2}\gamma_{31}}{N_{3}}\\ \frac{\gamma_{11}\gamma_{12}}{N_{1}}&\frac{\gamma_{21}\gamma_{22}}{N_{2}}&\frac{\gamma_{31}\gamma_{32}}{N_{3}}\end{array}\right]$$

$${\mathcal{J}}_{ik}=\sum _{j=1,2}{\delta }_j\left[\frac{1}{3}{U}_{ji}^3{U}_{jk}+U_{jk}^3{U}_{ji}+\frac{2}{3}\frac{(U_{ji}{U}_{jk}{U}_{j3}{)}^2}{U_{j1}{U}_{j2}}\right],$$

$${\mathcal{K}}_{ik}=\sum _{j=1,2}{\delta }_j{U}_{ji}^3{U}_{jk},{\mathcal{S}}_{ikl}=\frac{4}{3}\sum _{j=1,2}{\delta }_j\frac{(U_{ji}{U}_{jk}{U}_{j3}{)}^2}{U_{ji}{U}_{jk}},$$

$$H=H_0+ϵH_{int},H_0\equiv \sum _{i=1,2,3}{\omega }_i{\hat{\alpha }}_i^{†}{\hat{\alpha }}_i,$$

$$H_{int}\equiv \sum _{k=1,2}{g}_k\left({\hat{\alpha }}_3{\hat{\alpha }}_k^{†}+{\hat{\alpha }}_3^{†}{\hat{\alpha }}_k\right)$$

$$S_1=-\sum _{k=1,2}{g}_k\left(\widehat{{\alpha }_3}{\widehat{{\alpha }^{\prime }}}_k^{†}-{\hat{\alpha }}_3^{†}{\widehat{{\alpha }^{\prime }}}_k\right).$$

$$H_{BD}=H_0-\sum _{i,j=1,2;(i\ne j)}\frac{g_i{g}_j}{2}\left(\widehat{{\alpha }_i}{\widehat{{\alpha }^{\prime }}}_j^{†}+{\hat{\alpha }}_i^{†}{\widehat{{\alpha }^{\prime }}}_j\right)$$

$$\mathbf{M}=\left[\begin{array}[]{ccc}\omega^{2}&0&V\\ 0&\omega^{2}&V\\ V&V&\omega_{r}^{2}\end{array}\right]$$

$$\omega ,\sqrt{\frac{{\omega }^2+{\omega }_r^2\pm \sqrt{({\omega }^2-{\omega }_r^2{)}^2+8V^2}}{2}}$$

$$\omega ,\omega \sqrt{1\pm \frac{2\sqrt{2}g}{\omega }}$$

$$\omega ,\omega \pm \sqrt{2}g$$

$${\hat{a}}_n=\sum _m{A}_{nm}{\hat{\alpha }}_m+B_{nm}{\hat{\alpha }}_m^{†}$$

$$A_{nm}=\left(\sqrt{\frac{{\omega }_n}{8{\tilde{\omega }}_m}}+\sqrt{\frac{{\tilde{\omega }}_m}{8{\omega }_n}}\right)S_{nm},$$

$$B_{nm}=\left(\sqrt{\frac{{\omega }_n}{8{\tilde{\omega }}_m}}-\sqrt{\frac{{\tilde{\omega }}_m}{8{\omega }_n}}\right)S_{nm}$$

$${\left(a_n-a_n^{†}\right)}^4=$$

$$6\sum _{m=1}^3{\left(A_{nm}-B_{nm}\right)}^4\left[{\left({\hat{\alpha }}_m^{†}{\hat{\alpha }}_m\right)}^2+{\hat{\alpha }}_m^{†}{\hat{\alpha }}_m\right]$$

$$+6\sum _{m<k}{\left(A_{nm}-B_{nm}\right)}^2{\left(A_{nk}-B_{nk}\right)}^2\left[{\hat{\alpha }}_m^2{\hat{\alpha }}_k^{†2}+{\hat{\alpha }}_m^{†2}{\hat{\alpha }}_k^2+4{\hat{\alpha }}_m^{†}{\hat{\alpha }}_m{\hat{\alpha }}_k^{†}{\hat{\alpha }}_k+2{\hat{\alpha }}_m^{†}{\hat{\alpha }}_m+2{\hat{\alpha }}_k^{†}{\hat{\alpha }}_k\right]$$

$$+4\sum _{m\ne k}{\left(A_{nm}-B_{nm}\right)}^3\left(A_{nk}-B_{nk}\right)\left({\hat{\alpha }}_m^2{\hat{\alpha }}_m^{†}{\hat{\alpha }}_k^{†}+{\hat{\alpha }}_m^{†2}\widehat{{\alpha }_m}{\hat{\alpha }}_k+2{\hat{\alpha }}_m^{†}{\hat{\alpha }}_m\widehat{{\alpha }_m}{\hat{\alpha }}_k^{†}+2{\hat{\alpha }}_m^{†}{\hat{\alpha }}_m{\hat{\alpha }}_m^{†}{\hat{\alpha }}_k+\widehat{{\alpha }_m}{\hat{\alpha }}_k^{†}+{\hat{\alpha }}_m^{†}{\hat{\alpha }}_k\right)$$

$$+8\sum _{m\ne k\ne l}{\left(A_{nm}-B_{nm}\right)}^2\left(A_{nk}-B_{nk}\right)\left(A_{nl}-B_{nl}\right)\left({\hat{\alpha }}_m^2{\hat{\alpha }}_l^{†}{\hat{\alpha }}_k^{†}+{\hat{\alpha }}_m^{†2}\widehat{{\alpha }_l}{\hat{\alpha }}_k+2{\hat{\alpha }}_m^{†}{\hat{\alpha }}_m\widehat{{\alpha }_l}{\hat{\alpha }}_k^{†}+2{\hat{\alpha }}_m^{†}{\hat{\alpha }}_m{\hat{\alpha }}_l^{†}{\hat{\alpha }}_k+\widehat{{\alpha }_l}{\hat{\alpha }}_k^{†}+{\hat{\alpha }}_l^{†}{\hat{\alpha }}_k\right)$$

Full text

I Introduction

Quantum computation is rapidly progressing toward practical technology Lucero ; Kandala17 ; Cai ; Monz . So far quantum bits have been well developed on superconducting circuits Wilhelm . When cooled to milikelvin temperatures, coherent tunneling of the Cooper pairs through the Josephson junction (JJ) exhibits slightly nonlinear harmonic oscillations with addressable energy levels Gambetta16 . Such quantum states have long coherent times and can operate in nanosecond scales. Moreover, they are compatible with the well-established microwave control technology and can scale up in large numbers. All these features makes superconducting qubits one of the prominent platforms for constructing a multiqubit quantum processor Chow ; Niemczyk .

The state of art superconducting circuits contains a few tens of qubits with operational gate error rate about 0.1$\%$ for single qubit gates Barends and 1$\%$ for two-qubit gates Sheldon yet marginally below the threshold for error detection in the surface code Gambetta17 . Scaling up the qubit number to be more than a few tens will dramatically increases errors to unacceptable values and the key milestone of next few years is to reduce the errors Preskill ; Martinis18 . Achieving this not only requires further enhancements in the circuit quality Martinins16 but also needs progressive advancement in theory theory . So far the Jaynes-Cummings model, originally introduced in quantum optics JC , and its generalization have been routinely applied on superconducting circuits. These models have been so far well-studied for parameters admissible by perturbation theory, namely within ‘dispersive regime’ Wallraff11 ; Nataf . Moreover, special resonant frequency solutions are also known Haroche ; Blais04 . However, scaling up the number of qubits within the narrow domain of parameters will introduce new issues, such as circuit frequency crowding that must be avoided Plourde17 . Recently it has been discussed that significant advantages can be made in engineering circuits outside of the dispersive regime Bultkin . However, those studies have been performed numerically in the absence of established theory.

Here, motivated by the ‘black box quantization’ method, we develop a formalism for evaluating the qubit characteristics in circuits consisting of transmons and multimode resonators/cavities at arbitrary frequencies and coupling strengths. Black box quantization has been recently introduced in Ref. Nigg12 for circuits consisting of low anharmonic transmons coupled to resonators. The low anharmonicity allows dividing the circuit hamiltonian into harmonic and anharmonic sectors. In the absence of anharmonicity the transmons and the resonators can be treated on equal footing, thus the Foster decomposition Foster can replace the harmonic circuit with a set of lumped imaginary impedances seen by the anharmonic sector. Identifying the characteristic impedances is of the central importance in this method for which Ref. Nigg12 proposes iterative feedbacks between experiment and theory. This formalism has initiated so far several progressive improvements for extracting circuit parameters from electromagnetic simulation tureci17 ; Firat17 ; Richer17 . In Fig.(1) the harmonic sector is made of $N$ qubits (in blue boxes) coupled to cavity modes (in the gray area). The curly (red) crosses denote anharmonic sector.

After introducing a unitary transformation matrix in the space of total number of qubits and resonators, we find a normal-mode basis for the harmonic sector of a multiqubit circuit. Using the transformation we determine all dressed frequencies and Kerr nonlinear terms in the leading order of anharmonicity. The simplicity and accuracy of this method allows us to present he results of a transmon and two transmons in closed form formulas. For complex circuits this method provides insightful hamiltonian diagonalization inside and outside of the dispersive regime, which will be progressively useful in scaling up the number of qubits. We explain the formalism first in single transmon, then we generalize it before we solve another example of two transmons.

II Single transmon coupled to resonator

A transmon coupled to a resonator – The canonical variables are charges and phases Clerk , i.e. $({q}_{i},{\phi}_{i})$ with $i$ being $t,r$ for transmon and resonator. The transmon is coupled to the center conductor of resonator by the capacitance $C_{g}$. The dipole interaction $H_{int}=\beta V_{r}q_{t}$ couples the transmon charge and the resonator voltage $V_{r}=q_{r}/C_{r}$, with $C_{r/t}$ being the resonator/transmon capacitance and $\beta\equiv C_{g}/C_{t}$. Keeping $\beta\ll 1$ guarantees the increase in qubit coherence time Manucharyan . The circuit harmonic and anharmonic sectors sum to define the circuit classical hamiltonian:

[TABLE]

The characteristic impedances and the harmonic frequencies in the circuit are $Z_{i}=\sqrt{L_{i}/C_{i}}$ and $\omega_{i}=1/\sqrt{L_{i}C_{i}}$, respectively, with $E_{C}$ being total capacitive energy of transmon (including JJ and shunt capacitances as well as capacitive coupling between transmon and voltage sources), and $\hbar$ the reduced Planck constant. We define canonical variables $({Q}_{i},{X}_{i})\equiv({q}_{i}\sqrt{L_{i}},{\phi}_{i}/\sqrt{L_{i}})$ such that the harmonic part of Eq. (1) can be transformed to

[TABLE]

with $g\equiv\beta\omega_{r}\sqrt{Z_{r}/4Z_{t}}$ and $\mathbf{Q}\equiv(Q_{t},Q_{r})$ and $\mathbf{X}\equiv(X_{t},X_{r})$.

This Hamiltonian can be diagonalized by unitarily transforming $\mathbf{Q}$ and $\mathbf{X}$ into new canonical variable $\mathcal{Q}$ and $\mathcal{X}$, i.e. ${Q}_{i}=\sum_{j}S_{ij}\mathcal{Q}_{j}$ and ${X}_{i}=\sum_{j}T_{ij}\mathcal{X}_{j}$. Given that the variables in the new and the old frames must satisfy the Poisson brackets of canonical coordinates, i.e. $\{\mathcal{Q}_{i},\mathcal{X}_{j}\}=\{Q_{i},X_{j}\}=\delta_{ij}$, one can find that $T_{ij}=S_{ij}$ (see Appendix A). The only term in Eq. (2) that needs diagonalization is $\mathbf{Q}^{\textup{T}}\mathbf{MQ}$, which in the new basis must look like $\mathcal{Q}^{\textup{T}}\Omega\mathcal{Q}$, with $\Omega$ being a diagonal matrix, $\Omega_{tt}=\bar{\omega}_{t}^{2}$, $\Omega_{rr}=\bar{\omega}_{r}^{2}$, and zero otherwise. The unitary transformation $\mathbf{S}$ is therefore the matrix of columns of normalized eigenvectors of $\mathbf{M}$.

In the new basis the following dressed frequencies can be found in the linear sector: $\bar{\omega}_{t}\equiv K_{-}^{2}$ and $\bar{\omega}_{r}\equiv K_{+}^{2}$, with $K_{\pm}\equiv 2^{-\frac{1}{4}}\left(\omega_{t}^{2}+\omega_{r}^{2}\pm\Delta\Sigma s^{-1}\right)^{\frac{1}{4}}$ and $s\equiv({1+{16\left({g}/{\Delta}\right)^{2}\omega_{r}\omega_{t}}/{\Sigma^{2}}})^{-1/2}$, $\Sigma\equiv\omega_{r}+\omega_{t}$ and $\Delta\equiv\omega_{r}-\omega_{t}$. The unitary transformation matrix $\mathbf{S}$ is made of columns of the following normalized eigenvectors $(\pm\sqrt{({1\mp s})/2},\sqrt{({1\pm s})/2})^{T}$ associate to the eigenvalues $K_{\pm}^{2}$. (In Appendix H similar results have been found using the Bogoliubov transformations BV ) In this basis the anharmonic term proportional to $X_{t}^{4}$ should be transformed using the phase transformation $X_{t}=-\sqrt{(1+s)/2}\mathcal{X}_{t}+\sqrt{(1-s)/2}\mathcal{X}_{r}$, and this can make many types of terms possible, e.g. $C_{m}\mathcal{X}_{r}^{m}\mathcal{X}_{t}^{4-m}$ with coupling strengths $C_{m}(s)$ and $m=0,1,2,3,4$. In the original eigenbasis $|n_{j}\rangle$, with $j=t,r$, the ladder operators $\hat{a}_{j}=\sum_{n_{j}}\sqrt{n_{j}+1}\ket{n_{j}}\bra{n_{j}+1}$ can help to rewrite the charge operator $\hat{Q}_{j}=\sqrt{\hbar/2\omega_{j}}(\hat{a}^{\dagger}_{j}+\hat{a}_{j})$ and the phase operator $\hat{X}_{j}=i\sqrt{\hbar\omega_{j}/2}(\hat{a}^{\dagger}_{j}-\hat{a}_{j})$ Clerk . Similarly in the normal mode basis the ladder operators $\hat{\alpha}_{k}$ determine the new charge and phase operators: $\hat{\mathcal{Q}}_{j}$ and $\hat{\mathcal{X}}_{j}$. These two basis can be transformed into one another using the following Bogoliubov-Velatin transformation: $\hat{a}^{\dagger}_{t}-\hat{a}_{t}=U_{tt}(\hat{\alpha}^{\dagger}_{t}-\hat{\alpha}_{t})+U_{tr}(\hat{\alpha}^{\dagger}_{r}-\hat{\alpha}_{r})$, with $U_{tt}=-[{(1+s){\bar{\omega}_{t}}/{2}{\omega_{t}}}]^{\frac{1}{2}}$ and $U_{tr}=[{{(1-s)}{\bar{\omega}_{r}}/2{{\omega}_{t}}}]^{\frac{1}{2}}$.

The anharmonic quantum hamiltonian from Eq. (1) can be written as $H_{\textup{anhar.}}=-(\delta/12)(\hat{a}^{\dagger}_{t}-\hat{a}_{t})^{4}$, with $\delta\equiv E_{C}$ being the anharmonicity coefficient. In the new basis, this hamiltonian is transformed to $-\frac{\delta}{12}[U_{tt}(\hat{\alpha}^{\dagger}_{t}-\hat{\alpha}_{t})+U_{tr}(\hat{\alpha}^{\dagger}_{r}-\hat{\alpha}_{r})]^{4}$, defining self-Kerr coefficient Boissonneault of the transmon ${\chi}_{t}=\delta U_{tt}^{4}$ and that of the resonator ${\chi}_{r}=\delta U_{tr}^{4}$. Note that the anharmonic hamiltonian is not diagonal in the normal mode basis, however we can simplify it by ignoring irrelevant terms to first order and applying secular approximation. This reformulates total hamiltonian to

[TABLE]

The transmon state in normal mode basis makes a shift proportional to ${\chi}_{rt}$, namely the cross-Kerr coefficient, in the resonator frequency. It is simple to show that ${\chi}_{rt}=\sqrt{{\chi}_{r}{\chi}_{t}}$ and therefore it linearly scales with the anharmonicity $\delta$. Defining ‘dressed frequency’ $\tilde{\omega}_{i}$ to be the coefficient of $\hat{\alpha}_{i}^{\dagger}\hat{\alpha}_{i}$, after summing over all relevant terms we find the following closed form formula for the dressed frequencies:

[TABLE]

which indicate $E_{n_{t}n_{r}}=\sum_{i}\tilde{\omega}_{i}n_{i}-\chi_{i}^{2}n_{i}^{2}/2-2\chi_{rt}n_{t}n_{r}$. The validity of these formulas are much wider than dispersive regime, in fact they are valid for arbitrary coupling strength and frequency detuning.

Let us compare our results with other models. In general two sets of analytical results are known for the circuit: i) within the dispersive regime, which defines the validity of perturbation theory, or ii) at resonant frequencies in the Jaynes-Cummings model. In dispersive regime the detuning frequency is much stronger than the coupling strength, i.e. $g/\Delta\ll 1$, and within this regime Eqs. (4) and (5) are expanded in polynomials of $g/\Delta$. This will result in the following dressed frequencies: $\tilde{\omega}_{t}\approx\omega_{t}-2\omega_{r}{g}/{\Delta}\Sigma$, $\tilde{\omega}_{r}\approx\omega_{r}+2{g^{2}}\omega_{t}/{\Delta}\Sigma$ and the self-Kerr nonlinearities ${\chi}_{t}=\delta[1-4{g}^{2}\omega_{r}(\omega_{r}^{2}+\omega_{t}^{2})/\omega_{t}\Sigma^{2}{\Delta}^{2}]$ and ${\chi}_{r}=16\delta\left({g}/{\Delta}\right)^{4}\omega_{r}^{4}/\Sigma^{4}$. These expression are in agreement with the non-Rotating Wave Approximation (non-RWA) results recently reported in Gely using perturbation theory. In the RWA regime $g\ll\Delta\ll\Sigma$, therefore we can simplify these expressions by using Gely : $\omega_{r}/\Sigma\approx\omega_{t}/\Sigma\approx 1/2$. The RWA dressed frequencies are $\tilde{\omega}_{t}^{\textup{RWA}}\approx\omega_{t}-{\delta}/{2}-{g^{2}}/{\Delta}-\delta g^{2}/\Delta^{2}$, $\tilde{\omega}_{r}^{\textup{RWA}}\approx\omega_{r}+{g^{2}}/{\Delta}-\delta g^{2}/\Delta^{2}$, ${\chi}_{t}^{\textup{RWA}}\approx\delta[1-2({g}/{\Delta})^{2}]$, ${\chi}_{r}^{\textup{RWA}}\approx\delta\left({g}/{\Delta}\right)^{4}$. These results are in agreement with the original perturbative Lamb and AC-Stark shifts reported by J. Koch, et.al. in Ref. Koch07 and experimentally observed Fragner .

Fig(2a) shows the transmon and the resonator dressed frequencies at different detuning frequencies $\Delta$ and a fixed coupling strength $g$. For the choice of circuit parameters, there is a negligible mismatch between RWA and non-RWA perturbative results, therefore, in the logarithmic scales the lines labelled by Perturbation can be plotted using both formulations. In large detuning $\Delta$ the exact dressed frequencies of Eq. (4), in solid lines, are in good agreement with the perturbative (dotted) results. However, as $g/\Delta$ increases almost above $\sim$ 1/3 the difference between exact and perturbative results starts to appear. Another regime of interest is the special solution of the resonant point, where bare frequency of the transmon and the resonator are the same. Perturbation theory at this point diverges, however, the Jaynes-Cumming model predicts that due to atom-photon coupling a $2g$ frequency gap between the two dressed frequencies is produced Haroche ; Blais04 . Our exact formalism in Eq. (4) not only confirms this result but also provides a modification in it due to the presence of finite anharmonicity in transmons, which makes the dressed frequency detuning $2g(1-\delta/4\omega_{r})$. Fig. (2b inset) shows this gap rescaled by $2g$ at different coupling strengths $g$, which is unity for infinite anharmonicity (labelled J.C.); and is nonlinear for finite anharmonicity. Fig. (2b) shows cross-Kerr coefficient –defined below Eq. (3)– in solid line and compares it with the perturbative results in dashed lines. At the resonant point as expected the perturbation theory diverges, however, in contrast, the exact solution reveals the finite value of $\delta/4+o(g^{2})$ for any choice of bare frequencies.

III General method: N-atoms coupled to a resonator

$N$ transmons coupled to $M$ resonator. – Black box quantization in the original form has been proposed as an experimental method to get theoretical feedback on fitting parameters. Here we study a purely theoretical approach to generalize it to $N$ modes and $M$ transmons. This will help not only to scale up quantum circuits, but also it helps to study nonperturbative solutions.

Scaling up entanglement is one of the purely quantum phenomenon that is most crucial for quantum computing. Such phenomenon can take place in large scale quantum circuits with $N$ transmons interacting with $M$ resonators. A total of $N+M$ pairs of canonical variables can be defined: the charge vector $\mathbf{Q}=\left(Q_{1},\cdots,Q_{N+M}\right)^{T}$ and the phase vector $\mathbf{X}=\left(X_{1},\cdots,X_{N+M}\right)^{T}$. The circuit hamiltonian can be divided into a harmonic sector and a weakly anharmonic sector. The harmonic hamiltonian is $H_{\textup{har.}}=\frac{1}{2}\sum_{i=1}^{N+M}\omega_{i}^{2}Q_{i}^{2}+\frac{1}{2}X_{i}^{2}+\sum_{i=1}^{N}\sum_{j=N+1}^{N+M}g_{ij}\sqrt{4\omega_{i}\omega_{j}}Q_{i}Q_{j}$. Using a generalization of $\mathbf{M}$ matrix in Eq. (2) this hamiltonian is simplified to $H_{\textup{har.}}=\frac{1}{2}\mathbf{Q}^{\textup{T}}\mathbf{MQ}+\frac{1}{2}\mathbf{X}^{\textup{T}}\mathbf{X}$ with the matrix $\mathbf{M}$ being nonzero only at $M_{RR}=\omega_{r}^{2}$, $M_{TT}=\omega_{a}^{2}$, $M_{TR}=M_{RT}=g_{t}\sqrt{4\omega_{t}\omega_{r}}$; Sub-indices T labels transmons $\{1,2,\cdots,N\}$ and R the resonators $\{N+1,\cdots,M\}$. Consider that the following unitary transformations charges ${Q}_{i}=\sum_{j}S_{ij}{\mathcal{Q}}_{j}$ and phases ${X}_{i}=\sum_{j}T_{ij}{\mathcal{X}}_{j}$ take them to a normal mode basis. As discussed in Appendix A, these unitary transformations are identical, i.e. $T_{ij}=S_{ij}$. They transform the harmonic hamiltonian to $\frac{1}{2}\sum_{i}\bar{\omega}_{i}\mathcal{Q}_{i}^{2}+\mathcal{X}_{i}^{2}$. Detailed analysis show that $S$ is the matrix of normalized eigenvectors of $\mathbf{M}$-matrix. This evaluates dressed frequencies in the absence of anharmonicity, which in this paper we determine them exactly in closed form formula for circuits with one and two transmons and one resonator, however for larger circuits the M-matrix can be evaluated numerically and this determines all exact dressed frequencies.

Once $S$ is found charges and phases can be promoted to operators and rewritten in terms of ladder operator $\hat{a}$ in the original basis and $\hat{\alpha}$ in the new basis. They transform to one another as follows:

[TABLE]

The anharmonic hamiltonian $\sum_{i=1}^{N}(\delta_{i}/12)(\hat{a}_{i}-\hat{a}_{i}^{\dagger})^{4}$ can be similarly taken to the normal mode basis –see appendix H for details. The smallness of anharmonicity in transmons allows the nonlinear physical parameters to be evaluated in leading order.

IV Two transmons coupled to a resonator

Two transmons sharing a bus resonator – This is an important circuit for two-qubit gate calibration Blais07 ; Majer07 . Let us denote $\omega_{i}$ with $i=1,2,3$ for the two transmons, and the resonator, respectively –alternatively we sometimes use $r$ (instead of $3$) to emphasize the resonator. The coupling strengths between the transmons and the resonator are $g_{1},g_{2}$. The $\mathbf{M}$-matrix is

[TABLE]

The M-matrix can be taken to a normal mode basis within a wide domain of parameters that includes the superconducting circuits of interest for quantum computation –see below Eq. (8). The cubic equation $\lambda^{3}+b\lambda^{2}+c\lambda+d=0$ determines the eigenvalues $\lambda$ of Eq. (7) with $b\equiv-\sum_{i=1,2,3}\omega_{i}^{2}$, $c\equiv\omega_{1}^{2}\omega_{2}^{2}+\omega_{1}^{2}\omega_{3}^{2}+\omega_{2}^{2}\omega_{3}^{2}-\sum_{i=1,2}4g_{i}^{2}\omega_{i}\omega_{3}$, and $d\equiv 4g_{2}^{2}\omega_{1}^{2}\omega_{2}\omega_{3}+4g_{1}^{2}\omega_{1}\omega_{2}^{2}\omega_{3}-\omega_{1}^{2}\omega_{2}^{2}\omega_{3}^{2}$. Note that the eigenvalues $\lambda_{k}$ determine the circuit dressed frequencies, i.e. $\bar{\omega}_{k}\equiv\sqrt{\lambda_{k}}$, with $k=1,2,3$. By defining $\theta\equiv\lambda+b/3$ the quadratic term is eliminated, i.e. $\theta^{3}-f\theta+h=0$, $f\equiv b^{2}/3-c$ and $h\equiv\left(2b^{3}-9bc+27d\right)/27$. We solve this equation using trigonometric trials functions and find

[TABLE]

A relabelling of indices might be needed to identify corresponding frequencies. Should the inverse cosine be always between -1 and 1, ${h^{2}}/{4}-{f^{3}}/{27}<0$ must be satisfied for real-valued solutions. In circuits suitable for quantum computation, however, since coupling strengths are much smaller than individual frequencies this condition is trivially satisfied Appendix B.

The anharmonic hamiltonian transformed into the normal mode basis in the leading order can be written as

[TABLE]

with self-Kerr $\chi_{i}=\sum_{j=1,2}\delta_{j}U_{ji}^{4}$ and cross-Kerr $\chi_{ik}=\sum_{j=1,2}\delta_{j}U_{ji}^{2}U_{jk}^{2}$ and $U_{ij}$ being define in Eq. (6) —see Appendix C for details. One can evaluate the Kerr cofactors and see that in general there is no simple relation between cross Kerr and self-Kerr coefficients.

In Eq. (9) the $\mathcal{J}$ coupling indicates a direct interaction between two oscillators. The $\mathcal{K}$ and $\mathcal{L}$ terms are multiplied by $\hat{\alpha}^{\dagger}_{i}\hat{\alpha}_{i}$, therefore they are effectively $n_{l}\mathcal{L}_{ikl}$ and $n_{i}\mathcal{K}_{ik}$ with $n$ being integer quantum numbers. These couplings linearly depend on anharmonicity $\delta$ and are stronger in higher excited states. Detailed analysis —see Appendix D and E— show that after block-diagonalization these three interactions appear in the effective hamiltonian of the two-qubit circuit only in the higher orders $\delta^{2}$, thus they are negligible in the leading order.

The dressed frequency of transmons and resonator are

[TABLE]

and the energy levels are $E_{n_{1}n_{2}n_{3}}=\sum_{i=1}^{3}n_{i}\tilde{\omega}_{i}-\chi_{i}n_{i}^{2}/2-2\sum_{k>i}\chi_{ik}n_{i}n_{k}$. Fig. (3a) shows all dressed frequencies given the bare values $\omega_{1}=3$GHz, $\omega_{2}=\alpha\omega_{1}$, and $\omega_{r}=\omega_{1}+\Delta_{r1}$. We obtain perturbative results (dotted) using the formalism explained in Ref. Magesan ; Gambetta13 —more explicitly Eqs. (4.3-4.5) of the first reference. For the fixed coupling strength $g$ and large $\Delta_{r1}\gg g$ the results of perturbation theory and Eq. (10) are in good agreement. Far from the dispersive regime, however, the two become much deviated. For example in a circuit with the frequency of transmon 1 and resonator in resonance perturbation theory diverges, however Eq. (10) predicts a finite dressed frequency gap as shown in Fig. (3a). In the case of $\alpha<1$ and $g/\omega_{1}<(1-\alpha^{2})/3\sqrt{6}$, a series expnasion of Eq. (10) in terms of $g$ shows that in absence of anharmonicity $\tilde{\Delta}_{r1}\approx 2g$ and $\tilde{\Delta}_{12}\approx\Delta_{12}-g+[(1+\alpha)^{-1}+(1+\alpha)/2]g^{2}/\Delta_{12}+o(g^{3})$. Fig. (3b) shows the special case of resonant transmons with $\omega_{1}=\omega_{2}\equiv\omega$. Our exact evaluation indicates that all dressed frequencies become off resonant at the values $\omega$, $\sqrt{({\omega^{2}+\omega_{r}^{2}\pm\Delta\Sigma r^{-1}})/2}$ with $r^{-2}\equiv 1+32g^{2}\omega\omega_{r}/\Delta_{r1}^{2}\Sigma_{r1}^{2}$. The special case of maximal resonance, i.e. $\omega_{r}=\omega$, the dressed frequencies will be found $\tilde{\Delta}_{r1}\approx\tilde{\Delta}_{12}\approx\sqrt{2}g$ in absence of anharmonicity, see Appendix F for details. These are all previously unknown results that can be used for instance to identify bad samples in circuit fabrications.

Before conclusion let us comment on evaluating the effective impedances introduced by the black box quantization Nigg12 . They are assumed to be unknown parameters and can be evaluated in iterative feedbacks between theory and experiment. Our formalism, however, reveals a number of constraints that link between the effective impedances that makes the simpler to be theoretically estimated. For the simple example of a transmon coupled to a resonator, the effective impedances can be found analytically: $Z_{t}^{\textup{eff}}=\tilde{\omega}_{t}{(1+s)Z_{t}}/{2\omega_{t}}$ and $Z_{r}^{\textup{eff}}=\tilde{\omega}_{r}{(1-s)Z_{t}}/{2\omega_{r}}$. Their ratio in the dispersive regime is $Z_{r}^{\textup{eff}}/Z_{t}^{\textup{eff}}\sim(g/\Delta)^{2}$, which indicates the characteristic impedance associated to the transmon exceeds that of the resonator.

V Discussion

We presented a rigorous method to exactly obtain effective qubit parameters from the hamiltonian of superconducting circuits consisting of resonators and JJs at arbitrary coupling strengths and frequency detunings. Using this formalism we exhibited single transmon and two transmons outside of the dispersive regime in closed-form formulas. For complicated circuits finding analytical expressions may not be easy, however, our formalism can determine qubit parameters numerically much easier and more accurate compared to perturbation theory in charge basis, because the M-matrix (defined in the text) linearly scales with the number of qubits and resonators and all we need is to find its eigenvectors. This simplicity and accuracy will play an essential role for scaling up superconducting circuits as it allows to explore the possibilities of new domains of parameters for elevated fidelities.

Acknowledgements.

We thank David DiVincenzo for many useful discussions. Support from Intelligence Advanced Research Projects Activity (IARPA) under contract W911NF-16-0114 is gratefully acknowledged.

Author contributions statement

M.A. conducted theoretical analyses.

Additional information

Competing interests The authors declare no competing interests.

Appendix A Unitary transformation of canonical variables

Consider two $N$ dimensional vectors of canonical variables $\mathbf{q}=(q_{1},q_{2},\cdots,q_{N})$ and $\mathbf{p}=(p_{1},p_{2},\cdots,p_{N})$. These variables satisfy the Poisson bracket relation $\{q_{i},p_{j}\}=\delta_{ij}$ with $i,j=1,2,\cdots,N$ and the definition of $\{f,g\}=\sum_{i=1}^{N}(\partial f/\partial q_{i})(\partial g/\partial p_{i})-(\partial f/\partial p_{i})(\partial g/\partial q_{i})$.

Let us consider the following unitary transformations takes place on these variables: $Q_{i}=\sum_{j=1}^{N}S_{ij}q_{j}$ and $P_{i}=\sum_{j=1}^{N}T_{ij}p_{j}$. In order to have the two new variables $\mathbf{Q}$ and $\mathbf{P}$ to be canonical variables they must satisfy similar Poisson bracket relation as those of old variables: $\{Q_{i},P_{j}\}=\delta_{ij}$. This indicates that $\{Q_{i},P_{j}\}=\sum_{k=1}^{N}(\partial Q_{i}/\partial q_{k})(\partial P_{j}/\partial p_{k})-(\partial Q_{i}/\partial p_{k})(\partial P_{j}/\partial q_{k})$. One can easily simplify these relations into: $\sum_{k=1}^{N}S_{ik}T_{jk}=\delta_{ij}$. Because of the unitarity of the transformation matrices $S$ and $T$ one can see that $\sum_{k=1}^{N}S_{ik}S_{kj}^{\dagger}=\delta_{ij}$. For real matrices we have $S^{\dagger}_{kj}=S_{jk}$, thus $\mathbf{T}=\mathbf{S}$.

Appendix B Constraints within exact formula for 2 transmon circuit

Another condition that can be concluded from Eq. (8) of the main article is the following:

[TABLE]

By definition we have always $b\leq 0$, therefore the condition can be checked in the cases where $\cos$ function is negative, therefore we need to check the following condition: $-2\sqrt{f/3}+|b|/3\geq 0$, which can be further simplified to ${b^{2}}/{3}>c$. Substituting the definitions will introduce the following condition to hold:

[TABLE]

We take first three terms from right side to the left, then simplify left side to arrive at the following condition:

[TABLE]

which trivially holds valid without imposing any limitations on parameters.

Appendix C Unitary transformation for 2 transmons coupled to resonator

The unitary transformation to diagonal basis in the harmonic sector is carried out by the matrix of normalized eigenstates with columns being eigenvectors, which is

[TABLE]

with $V_{i}\equiv g_{i}\sqrt{4\omega_{r}\omega_{i}}$, $\gamma_{ij}=\bar{\omega}_{i}^{2}-\omega_{j}^{2}$, and $N_{i}=\sqrt{V_{2}^{2}\gamma_{i1}^{2}+V_{1}^{2}\gamma_{i2}^{2}+\gamma_{i1}^{2}\gamma_{i2}^{2}}$.

Appendix D Additional interaction terms

In the circuit made of two transmons coupled to a shared resonator, the anharmonic part of Hamiltonian can be simplified to Eq. (9) in the main article. Below are detailed interaction couplings in terms of bare parameters:

[TABLE]

Appendix E Block diagonalization

Let us consider the the Hamiltonians of two harmonic oscillators (labeled as 1, 2) coupled to a resonator (labeled as 3):

[TABLE]

The unperturbed part $H_{0}$ in the eigenbasis of itself id diagonal, however $H_{int}$ is not. In general we may not be able to find a tranformation to fully diagonal matrix, but instead we can separate out a subset of states from the rest of the states. The Schrieffer-Wolff transformation is one way to block diagonalize the interacting Hamiltonian into low energy and high energy sectors. This usually takes place by transforming the Hamiltonian by the anti-hermitian operator $\exp{S}$ in the following way: $H_{BD}=\exp{(-S)}H\exp{S}$, which can be expanded into $H_{BD}=\sum_{n=0}[H,S]_{n}/n!$ with $[H,S]_{n+1}=[[H,S]_{n},S]$ and $[H,S]_{0}=H$. One can in principle assume a geometric series expansion of the transformation matrix: $S=\sum_{i=0}(\epsilon)^{i}S_{i}$; however given that the zeroth or order of $H_{BD}$ is $H_{{BD}_{0}}=[H_{0},S_{0}]=H_{0}$ therefore $S_{0}$ must be diagonal too which is in fact inconsistent with the definition of $S$ to be anti-hermitian and block-off-diagonal, therefore always $S_{0}=0$. In the first order the Hamiltonian is already given by $H_{int}$ which can be made of block-diagonal (bd) and block-off-diagonal (bod) matrices $H_{int}=H_{int}^{\textup{bd}}+H_{int}^{\textup{bod}}$ . Therefore $H_{BD_{1}}=[H_{0},S_{1}]=-H_{int}^{\textup{bod}}$. In the second order: $H_{BD_{2}}=[H_{0},S_{2}]+[H_{int},S_{1}]+(1/2)[[H_{0},S_{1}],S_{1}]$, and so on. Putting all together one can find the effective Hamiltonian up to the second order $H_{BD}=H_{0}+H_{int}^{\textup{bd}}+(1/2)[H_{int}^{\textup{bod}},S_{1}]$.

Using the relations above for the Hamiltonian of Eq. (LABEL:eq._H_SM) in which the interaction is block-off-diagonal one can use the following ansatz

[TABLE]

with $\hat{\alpha^{\prime}}_{k}\equiv\sum_{n=0}^{\infty}\sqrt{n+1}(\omega_{3}-\omega_{k})^{-1}|n\rangle\langle n+1|$ being the modified ladder operator for $k$-th transmon, given that the normal ladder operator for the same transmon is $\hat{\alpha}_{k}\equiv\sum_{n=0}^{\infty}\sqrt{n+1}|n\rangle\langle n+1|$ .

One can explicitly determine the effective Hamiltonian up to the second order of perturbation theory becomes

[TABLE]

Appendix F Resonant transmons

In a circuit with two transmons in resonance $\omega_{1}=\omega_{2}\equiv\omega$ and homogeneous coupling and anharmonicity $g_{1}=g_{2}\equiv g$ and $\delta_{1}=\delta_{2}\equiv\delta$ the harmonic Hamiltonian is $H_{\textup{har.}}=\frac{1}{2}\omega^{2}(Q_{1}^{2}+Q_{2}^{2})+\frac{1}{2}\omega_{r}^{2}Q_{r}^{2}+\frac{1}{2}(X_{1}^{2}+X_{2}^{2}+X_{r}^{2})+g\sqrt{4\omega\omega_{r}}(Q_{1}+Q_{2})Q_{3}$. Defining the vectors $\mathbf{Q}=(Q_{1},Q_{2},Q_{r})^{T}$ and $\mathbf{P}=(P_{1},P_{2},P_{r})^{T}$, this Hamiltonian can be rewritten as $H_{\textup{har.}}=\frac{1}{2}\mathbf{Q}^{\textup{T}}\mathbf{MQ}+\frac{1}{2}\mathbf{X}^{\textup{T}}\mathbf{X}$ with the matrix $\mathbf{M}$ being

[TABLE]

with $V\equiv g\sqrt{4\omega\omega_{r}}$. Because the off diagonal elements are identical, it is easy to find the eigenvalues, which are

[TABLE]

At the exterm resonance with $\omega_{r}=\omega$ the eigenenergies will become

[TABLE]

In the limit of small coupling $g\ll\omega$ this can be simplified to

[TABLE]

Appendix G Anaharmonicity

Consider the following Bogoliubov transformations for transmon ladder operator:

[TABLE]

and using the relation between transmon charge number and phase and the ladder operator $\hat{a}_{n}=\sqrt{\frac{\omega_{n}}{2}}\hat{q}_{n}+i\sqrt{\frac{1}{2\omega_{n}}}\hat{p}_{n}$, and its conjugate as well as similar in the transformed basis $\hat{\alpha}_{n}=\sqrt{\frac{\tilde{\omega}_{n}}{2}}\hat{Q}_{n}+i\sqrt{\frac{1}{2\tilde{\omega}_{n}}}\hat{P}_{n}$, one can find

[TABLE]

in which $\tilde{\omega}$ is the frequency in the transformed basis.

The anharmonicity in Hamiltonian will be $-\frac{\delta_{i}}{12}\left(a_{i}-a_{i}^{\dagger}\right)^{4}$. The operator part can be Bogoliubov transformed to the new basis, keeping terms with as many creations as annihilations, ignoring frequencies:

[TABLE]

Appendix H Bogoliubov transformation for Hamiltonian diagonalization

In this section we use quantum Hamiltonian of a transmon coupled to a resonator is $H=4E_{c}n-E_{J}\cos\phi+H_{\textup{res}}$. Separating the harmonic sector and the anharmonic sector, and using Bogoliubov transformation we diagonalize the interacting harmonic sector into a diagonal quantum harmonic Hamiltonian. We find all Bogoliubov transformation coefficients, which turns out to be similar to the results we took from semiclassical analysis.

Given that charge number operator is proportional to ladder operators $n\sim 2^{-\frac{1}{4}}\left(a+a^{\dagger}\right)$ and phase is the conjugate variable $\phi\sim 2^{\frac{1}{4}}\left(a-a^{\dagger}\right)$, and the resonator Hamiltonian is $H_{\textup{res}}=\omega_{r}b^{\dagger}b$, the circuit Hamiltonian can be written as $H=\omega_{q}a^{\dagger}a-\frac{\delta}{12}\left(a-a^{\dagger}\right)^{4}+\omega_{r}b^{\dagger}b+g\left(a+a^{\dagger}\right)\left(b+b^{\dagger}\right)$ with harmonic part being $H_{\text{har}}=\omega_{q}a^{\dagger}a+\omega_{r}b^{\dagger}b+g\left(a+a^{\dagger}\right)\left(b+b^{\dagger}\right)$.

We would like to Bogoliubov-transform the Hamiltonian into a diagonal Hamiltonian $\mathcal{H}$:

[TABLE]

We use a technique widely used in second quantized QFT, which is to Bogoliubov-transofmation creation and annihilation operators

[TABLE]

Eight equations are needed to determines coefficients; four by enforcing that transformed Hamiltonian preserves eigenvalues, which is equivalent to equating $H_{ho}$ and $\mathcal{H}_{ho}$ and setting coefficients of $\hat{\alpha}\hat{\alpha},\hat{\beta}\hat{\beta},\hat{\alpha}\hat{\beta}$ and $\hat{\alpha}\hat{\beta}^{\dagger}$to zero, respectively:

[TABLE]

The other four are determined by enforcing commutation relations, i.e. $[a,a^{\dagger}]=[b,b^{\dagger}]=1$ and $[a,b]=[a,b^{\dagger}]=0$, respectively, given that $[\alpha,\alpha^{\dagger}]=[\beta,\beta^{\dagger}]=1$ and zero otherwise:

[TABLE]

For simplicity we assume coefficients are real-valued, but the equations are difficult to be analytically solved. A practical simplification can be achieved by defining new variables

[TABLE]

which reformulates equations given above to the followings:

[TABLE]

Given that one may solve the Bogoliubov coefficient equtions, we can determine new frequencies in $\mathcal{H}$:

[TABLE]

One can easily prove that $F_{+}F_{-}=A_{+}A_{-}$, which simplifies equations and helps to find the following two important equalities:

[TABLE]

Substituting them in Eq. (H) we find one equation between $A_{\pm}$:

[TABLE]

This is one of the main equations we need to solve. Another one can be determined taking some non-trivial steps listed below: We use Eq. (H), substitute $B_{\pm}$ from Eqs. (H,H), multiply two side in $E_{+}F_{-}^{2}F_{+}$ and simplify it, magically the final equation is again a second equation that relation $A_{\pm}$:

[TABLE]

Now we solve these two equations together. To do so we first define $x=A_{+}A_{-}$ and substitute in Eq. (27): $a(1-x)x-(2x-1)^{2}=0$ with $a\equiv\frac{\Delta^{2}\Sigma^{2}}{4g^{2}\omega_{r}\omega_{q}}$ and $\Sigma=\omega_{r}+\omega_{q}$ and $\Delta=\omega_{r}-\omega_{q}$. Exact real-valued solution is

[TABLE]

and substituting in Eq. (26) determines exact real-valued $A_{\pm}$:

[TABLE]

In order to find $F_{\pm}$ yet we need to simplify Eq. (H) by multiplying on both sides on $F_{-}F_{+}$ and rewriting $B_{\pm}$ in terms of $A_{\pm}$, $E_{\pm}$ and $F_{\pm}$ as shown in Eqs. (H,H):

[TABLE]

Defining

[TABLE]

then

[TABLE]

We can expand the functions in terms of small coupling $g$ to any order. Below are results up to the fourth order:

Substituting in definition of new frequencies one finds:

[TABLE]

In the weak interaction limit these frequecies turn into Lamb and Stark shifts. Below we evaluate them up to fourth order:

[TABLE]

Anharmonicity can be easily derived using the following relation:

[TABLE]