Non-perturbative analytical diagonalization of Hamiltonians with application to coupling suppression and enhancement in cQED

In this subsection, we introduce the NPAD for symbolic diagonalization of Hermitian matrices and discuss how it can be applied to obtain effective models.

In this algorithm, a Givens rotation is defined in each iteration to remove one specifically targeted off-diagonal term. By iteratively applying the rotations, the transformed matrix converges to the diagonal form. The rotation keeps the energy structure when the off-diagonal coupling is small while always exactly removing the coupling even when it is comparable to or larger than the energy gap. Compared to the Jacobi method used in numerical diagonalization [12, 14, 13], we truncate the iteration at a much earlier stage. As each iteration consists only of a few algebraic expressions, the algorithm produces a closed-form, parametric expression of the transformed matrix.

We start from a two-by-two Hermitian matrix and define a complex Givens rotation that diagonalizes it. Then, we generalize the rotation to higher-dimensional matrices, discuss the convergence of the iteration, and how to use it as a symbolic algorithm. In Section III.1, we show a concrete application where we apply NPAD with only two rotations to approximate the energy spectrum of a near-resonant quantum system which can not be studied perturbatively.

In the previous subsection, we introduced NPAD that produces a closed-form, parametric expression of an approximately diagonalized matrix. Here, we show that in the perturbative limit, where the coupling is much smaller than the bare energy difference, the Jacobi iteration reduces to a Schrieffer-Wolff-like transformation. Interestingly, the recursive nature of the Jacobi iteration is preserved in this limit. Instead of looking for one generator that diagonalizes the full matrix as in the traditional Schrieffer-Wolff transformation (SWT), an iteration is constructed such that every time only the leading-order coupling is removed. We refer to it as recursive Schrieffer-Wolff transformation (RSWT) because of the recursive expression it produces. We also show that RSWT demonstrates an exponential improvement in complexity compared to SWT for perturbation beyond the leading order. In Section III.2, we demonstrate an application of RSWT in estimating the ZZ interaction between two Transmon qubits in a dispersive regime.

Both the NPAD and the RSWT methods introduced in the previous sections can be designed to only target a selected set of off-diagonal terms and, hence, used for block-diagonalization. This is especially useful to engineer transversal coupling in a subsystem and leave the remaining levels as intact as possible. Here, we briefly discuss these generalizations. Notice that it is always possible to first diagonalize the matrix and then reconstruct the block diagonalized form that satisfies certain conditions, for instance as in Ref. [44]. In the following, we discuss only methods that do not diagonalize the matrix first.

In NPAD, by construction, each rotation removes one off-diagonal element. With Givens rotations only applied to the inter-block elements, an iteration for block diagonalization can be defined. The norm of all off-diagonal entries, $\norm{H}_{F}$, is still monotonously decreasing according to Eq. 11. Hence, a limit exists and its convergence is also the convergence of the block diagonalization. However, the convergence is not always monotonous with respect to the norm of all inter-block terms. This is because a Givens rotation may rotate a large intra-block term into an inter-block entry. Therefore, the algorithm may not always converge faster than the full diagonalization would. Nevertheless, if the dominant coupling terms in the Hamiltonian are known, the Jacobi iteration can be designed to target at those to realize an efficient block-diagonalization. In Section III.3, we show an example of this in computing the cross-resonance coupling strength through NPAD.

For perturbation, RSWT can be applied as a block diagonalization method under the constraint that both the inter-block and the intra-block coupling are much smaller than the inter-block energy gap. This can be achieved by slightly modifying the RSWT iterations: We first separate the diagonal, the intra-block and the inter-block terms: $H_{n}=D_{n}+V^{\rm{intra}}_{n}+V^{\rm{inter}}_{n}$. Next, in algorithm 2 we only define $S$ for those non-zero entries in $V^{\textnormal{inter}}_{n}$, i.e. the couplings we wish to remove. And in the last step we replace Eq. 20 with

In this definition, the leading inter block coupling is of the order $\mathcal{O}([S_{n+1},V^{\rm{intra}}_{n}])$. As we do not remove the intra-block coupling, we still get $V^{\rm{intra}}_{n}=\mathcal{O}(\lambda)$. Therefore, the remaining coupling is $\mathcal{O}(\lambda V^{\rm{intra}}_{n})$, i.e. the perturbation order is increased by one, instead of being squared as in the case of full diagonalization. Therefore, the exponential reduction of the number of commutators does not always apply in the case of block diagonalization. However, notice that the small parameter $\lambda$ here is defined as the (largest) ratio between the inter-block couplings and gaps, which is usually much smaller than those within the block. Hence, if carefully designed, the convergence can still surpass the full diagonalization in the first few perturbative orders.

To help understand the proposed methods, we here discuss the difference between them and the traditional methods. We first compare RSWT with traditional SWT and then NPAD with the perturbation methods.

For RSWT, with the same target accuracy, e.g., $\mathcal{O}(\lambda^{K})$, it should provide the same expression as from SWT, up to the error $\mathcal{O}(\lambda^{K+1})$. However, compared to the SWT, RSWT requires a much smaller number of iterations and commutators. To reach $\mathcal{O}(\lambda^{K})$, SWT needs $K-1$ iterations, while RSWT only needs $\lfloor\log_{2}(K)\rfloor$ because of the quadratic convergence rate. More importantly, the total number of commutators grows only linearly for RSWT, compared to the exponentially fast growth for SWT [7].

Intuitively, this is because RSWT uses the recursive structure and avoids unnecessary expansions of the intermediate results. Mathematically, this can be seen from the following two aspects: First, in RSWT, each iteration improves the perturbation level from $\lambda^{k}$ to $\lambda^{2k}$, instead of $\lambda^{k+1}$. Hence, the number of iterations increases only logarithmically with respect to the perturbation order, as seen in the definition of $n_{\rm{max}}$ in algorithm 2. This is because we always treat the transformed matrix as a new one and remove the leading-order coupling. It is consistent with the quadratic convergence rate of the Jacobi iteration with small off-diagonal terms. Second, in RSWT, the generator $S_{n}$ is only used at the current iteration. Hence, there are no mixed terms such as $[S_{2},[S_{1},V_{0}]]$, in contrast to SWT.

The total number of commutators required to reach level $\lambda^{K}$ is shown in Table 1, where we have taken into consideration that if $\mathcal{C}_{t}(A,B)$ is known, computing $\mathcal{C}_{t+1}(A,B)$ only requires one additional commutator. The detailed calculation is presented in appendix B.

The NPAD method, on the other hand, uses non-linear rotations to replace the linear perturbative expansion. More concretely, in the Jacobi iteration, by targeting only one coupling in each recursive iteration, the unitary transformation can be analytically expressed as a Givens rotation, thus avoiding the BCH expansion in Eq. 16. Therefore, it efficiently and accurately captures the non-perturbative interactions in the system.

To compare it with the perturbation methods, we estimate the number of operations required for NPAD in the perturbative regime. Assume we construct the Jacobi iteration from the $G$ coupling terms used in generating an $S$ in RSWT. Applying those unitaries is, to the leading order, the same as applying one RSWT iteration. A single Givens rotation on a Hamiltonian takes $\mathcal{O}(N)$ operations, where $N$ is the matrix size. Thus, the cost for computing the effective Hamiltonian after $G$ rotations is the same as computing one commutator $[S,\boldsymbol{\cdot}]$, up to a constant factor. Because the Givens rotation avoids the BCH expansion, there is no nested commutators and the total number operations is $\mathcal{O}(n_{\rm{max}}NG)$ with $n_{\rm{max}}$ the number of iterations in algorithm 2. Hence the number of terms scales logarithmically with respect to $K$ instead of linearly as for RSWT, i.e., a super-exponential reduction compared to SWT (Table 1). However, the non-linear expressions provided by NPAD are usually harder to simplify and evaluate by hand compared to the rational expressions obtained from perturbation.

From the above discussion, one can see that it is also straightforward to combine NPAD with perturbation. Instead of fully diagonalizing the matrix, the Jacobi iteration can be designed to remove only the dominant couplings and combined with perturbation methods to obtain simplified analytical expressions. In fact, this is often used implicitly in the analysis when, e.g., a strongly coupled two-level system is perturbatively interacting with another quantum system. The Jacobi iteration suggests that this can be generalized systematically to more complicated scenarios.

In this first application, we apply the NPAD method described in Section II.1 to study a model consisting of two directly coupled qubits in the near-resonant regime, where the ZZ interaction can be used to construct a control-Z (CZ) gate (see Figure 1 block B) [21, 22, 23, 24]. We show that, with two rotations, NPAD provides an improvement on the estimation of the interaction strength for at least one order of magnitude, compared to approximating the system as only a single avoided crossing between the strongly interacting levels, as is standard in the literature. In addition, if one of the non-computational bases is comparably further detuned than the other, the correction takes the form of a Kerr nonlinearity, with a renormalized coupling strength accounting for the near-resonant dynamics.

We consider the Hamiltonian of two superconducting qubits that are directly coupled under the rotating-wave [45, 46] and Duffing [47] approximations:

where $b_{q}$, $\omega_{q}$, $\alpha_{q}$ are the annihilation operator, the qubit bare frequency and the anharmonicity, respectively. The parameter $g$ denotes the coupling strength. In this Hamiltonian, the sum of the eigenenergies is always a constant $E_{10}+E_{01}=\omega_{1}+\omega_{2}$ because of the symmetry. Hence, the ZZ interaction comes solely from the interaction between the state $\ket{11}$ and the non-computational basis $\ket{20}$ and $\ket{02}$. If the frequency is tuned so that the state $\ket{11}$ is close to one of the non-computational states, the coupling will shift the eigenenergy, leading to a large ZZ interaction (Figure 3a).

For simplicity, we consider the Hilbert subspace consisting of $\ket{20}$, $\ket{11}$, $\ket{02}$ and write the following Hamiltonian

The parameters in the diagonal elements are given by $\delta=(\omega_{1}-\omega_{2}+\alpha_{1})/2$ and $\Delta=3(\omega_{1}-\omega_{2})/2-\alpha_{2}+\alpha_{1}/2$. To keep the result general, we use two different coupling strengths $g_{1}$ and $g_{2}$, although according to Eq. 24 they both equal $\sqrt{2}g$. Without loss of generality, we assume the state $\ket{02}$ is comparably further detuned from the other two, i.e. $\Delta>g_{j},\delta$. If in contrast $\ket{20}$ is further detuned, one can exchange the $\ket{02}$ and $\ket{20}$ in the matrix and redefine $\delta$ and $\Delta$ accordingly. Notice that this Hamiltonian is different from a $\Lambda$ system [3], where coupling exists only between far detuned levels.

To implement the CZ gate, one tunes the qubit frequency $\omega_{1}$ so that the states $\ket{11}$ and $\ket{20}$ are swept from a far-detuned to a near-resonant regime. Hence, the perturbative expansion diverges and cannot be used. A naive approach is to neglect the far-detuned state $\ket{02}$ and approximate the interaction as a single avoided crossing. In this case, $\zeta$ is approximated by

However, the interaction $g_{2}$ results in an error that, in the experimentally studied parameter regimes, can be as large as 10%, as shown in Figure 3b.

In the following, we show that with only two Givens rotations, one can obtain an analytical approximation, with the error reduced by one order of magnitude. The correction can be understood as a Kerr non-linearity with a renormalized coupling strength.

To get an accurate estimation of the ZZ interaction $\zeta$, we need to calculate the eigenenergy of $\ket{11}$ by eliminating its coupling with the other two states. Therefore, we will make two rotations sequentially on the entry $(0,1)$ and $(1,2)$, given by

where $U_{1}$ and $U_{2}$ are Givens rotations (Eq. 3) constructed for eliminating the entries $(0,1)$ and $(1,2)$. Because the matrix is real symmetric, the phase $\phi$ in Eq. 3 is 0.

The first transformed Hamiltonian, $H^{(1)}=U_{1}HU_{1}^{\dagger}$, takes the form

where $E_{2}=\delta\sqrt{1+\frac{g_{1}^{2}}{\delta^{2}}}$ is the eigenenergy for diagonalizing the two-level system of $\ket{20}$ and $\ket{11}$, consistent with Eq. 26. The notations used is the same as in Section II.1.3. In this frame, the coupling between $\ket{11}$ and $\ket{02}$ is reduced to $c_{01}g_{2}$, where $c_{01}$ is given by the non-linear expression

This non-linearity is crucial for the accurate estimation of the eigenenergy.

The second rotation further removes this renormalized coupling $c_{01}g_{2}$, giving

Including the new correction, $g_{2}c_{01}t_{12}$, the eigenenergy of state $\ket{11}$ reads

In Figure 3b, we plot the error of the estimated interaction strength $\zeta=H^{\prime}_{1,1}+\delta$ using typical parameters of superconducting hardware, compared to the numerical diagonalization $\tilde{\zeta}$. An improvement of at least one order of magnitude is observed compared to traditional methods.

Following the assumptions that $\Delta\gg\delta,g_{j}$, Eq. 31 simplifies to

We see that the correction takes the form of a Kerr non-linearity [48], but with a renormalized coupling strength $c_{01}g_{2}$. This non-linear factor $c_{01}$ accounts for the dynamics between $\ket{20}$ and $\ket{11}$ in the near-resonant regime. The same effect can be observed in higher levels where similar three-level subspaces exist. This approximation is plotted as a dashed curve in Figure 3b.

The error of this estimation comes both from the expansion of the square root in Eq. 31 as well as from the remaining coupling in $H^{(2)}$. The former can be approximated by the next order expansion

For the latter, we consider the remaining coupling in $H^{(2)}$ between $\ket{20}$ and $\ket{11}$, which reads $g_{2}s_{01}s_{12}$. In the limit $\Delta\gg\delta,g_{j}$, we have $s_{12}\leq\frac{\theta_{12}}{2}\leq\frac{g_{2}}{\Delta-E_{2}}\ll 1$, indicating that this coupling is much smaller than the energy difference. Hence, further correction can be estimated by

The contribution of the other remaining coupling between $\ket{20}$ and $\ket{02}$ is much smaller due to the large energy gap. Since $\epsilon_{2}$ is one order smaller than the $\epsilon_{1}$, $\epsilon_{1}$ will be the dominant error. We plot the region below this error in Figure 3b as a shaded background.

For the more general cases without assuming $\Delta\gg\delta,g_{j}$, it is hard to provide an error estimation due to the non-linearity. However, the result in Figure 3b indicates that Eq. 31 still shows a good performance in other parameter regimes commonly used in superconducting hardware, with an error smaller than 3%. We also observe that an improvement for another order of magnitude can be achieved by introducing a third rotation again on the entry $(0,1)$.

In this second example, we use the two methods to investigate the suppression of ZZ cross-talk with the qubit-resonator-qubit setup in the dispersive cQED regime, which corresponds to Figure 1 block A. We demonstrate that in the traditional setup without direct inter-qubit coupling, the ZZ interaction defined in Eq. 23 can still be zeroed in a quasi-dispersive regime by engineering the two parameters of qubit-resonator detuning. The zero points are described by an equation of a circle in the $\lambda^{4}$-perturbation. To accurately capture the interaction strength in the quasi-dispersive regime, we also compute with RSWT the $\lambda^{6}$-perturbation and show that the NPAD method with only 8 Givens rotations provides an expression with similar accuracy.

We consider a Hamiltonian of two superconducting qubits connected by a resonator:

Due to the finite detuning between the resonator and the qubits, a static ZZ interaction exists even if there is no additional control operation on the system. In order to implement high-quality quantum operations, this interaction needs to be sufficiently suppressed.

Several approaches have been developed to suppress the ZZ interaction. One way is to add a direct capacitive coupling channel in parallel with the resonator [34, 27, 33, 32, 30, 36, 49, 28, 31, 35, 29]. By engineering the parameters, the two interaction channels cancel each other. The interaction can either be turned on through a tunable coupler or through the cross resonant control scheme. The second approach is to choose a hybrid qubit system with opposite anharmonicity, which allows parameter engineering to suppress the ZZ interaction. One implementation is using a transmon and a capacitively shunt flux qubit (CSFQ) [25, 26]. Other methods include using additional off-resonant drive [39, 40, 41] and different types of qubits have also been proposed [38].

Most of the above works are based on the strong dispersive regime, where the resonator is only weakly coupled with the qubits. In this regime, the ZZ interaction strength $\zeta$ is only determined by the effective interaction with the two non-computational qubit state, $\ket{20}$ and $\ket{02}$ [7]

where $\Delta_{q}=\omega_{q}-\omega_{c}$ is the qubit-resonator frequency detuning, $\alpha_{q}$ the anharmonicity and $J_{jk,j^{\prime}k^{\prime}}$ the effective coupling strength between the physical qubit state $\ket{jk}$ and $\ket{j^{\prime}k^{\prime}}$. They are obtained by performing a leading order SWT and effectively decouple the resonator from the two qubits. In this regime, it is impossible to achieve zero ZZ interaction unless the two anharmonicity $\alpha_{q}$ adopt different signs.

However, Eq. 36 is only valid when ignoring the higher level of the resonator. If we reduce the qubit detuning $\Delta_{q}$ so that it becomes comparable with the anharmonicity $\alpha_{q}$, the second excited state of the resonator comes into the picture and can be used to suppress the ZZ interaction, also known as the quasidisq regime [10, 37]. We identify this regime as the quasi-dispersive regime because $g/\Delta_{q}$ is manufactured larger than $0.1$, e.g. in superconducting qubits with weak anharmonicity such as Transmons, though we show the same analysis can also hold for stronger anharmonicities. As a result, the calculation of $\zeta_{\textnormal{disp}}$ cannot be treated by only the leading-order SWT. In particular, we will see that, in the straddling regime, where $|\Delta_{1}-\Delta_{2}|<\alpha$, the interaction with the second excited resonator state leads to a $\lambda^{4}$-perturbative correction that can be used to suppress the ZZ interaction.

In the following, we first use the $\lambda^{4}$-perturbation to qualitatively understand the energy landscape and then investigate the higher-order corrections. For the $\lambda^{4}$-perturbation, using RWST, we only need 2 iterations and evaluate 4 commutators instead of 3 iterations and 11 commutators, as for traditional perturbation (Table 1).

In fact, the traditional approach that first approximates the system as an effective qubit-qubit direct interaction and then applies another perturbation to obtain the ZZ strength is also a two-step recursion [7]. However, for simplicity, it neglects the resonator states in the second perturbation. As detailed in appendix C, adding the resonator states, we obtain a better estimation for the quasi-dispersive regime. The result is consistent with the diagrammatic techniques used in [50, 25].

To illustrate the energy landscape, we write the interaction strength as

with $\Delta_{-}=\Delta_{1}-\Delta_{2}$. The first two terms coincide with Eq. 36 in the strong dispersive regime, up to $\mathcal{O}(\frac{g^{4}}{\Delta^{3}})$.

Assuming $\alpha=\alpha_{1}=\alpha_{2}$ and set $\zeta^{(4)}=0$ in Eq. 37, we obtain an equation of a circle that describes the location of the zero points

where $\Delta_{+}=\Delta_{1}+\Delta_{2}$ and $\Delta_{-}=\Delta_{1}-\Delta_{2}$. In this $\lambda^{4}$-perturbation, the zero-points depend only on the anharmonicity $\alpha$ but not on the coupling strength $g_{q}$. Equation 38 indicates that the ZZ interaction can be suppressed by varying the sum and difference of the two qubit-resonator detunings, as illustrated in Figure 4a. Because the perturbative approximation is only valid away from the resonant lines, the useful part of the parameter regime is the half-circle with $\Delta_{+}<\alpha$, in particular, the region marked by the grey box in Figure 4a.

[figure]style=plain,subcapbesideposition=top

In addition, we also studied different contributions to the ZZ interaction. In Figure 5, we plot the strong dispersive approximation, the $\lambda^{4}$-perturbation as well as the contribution of second excited qubit and resonator state to $\zeta^{(4)}$ (see appendix C for analytical expressions). One observes from the plot that, in the quasi dispersive regime, the increasing virtual interaction with the second excited resonator state acts against the interaction with the second excited qubit states. Notice that all contributions to $\zeta^{(4)}$ are virtual interactions of the second excited state, i.e., $\zeta^{(4)}=\zeta_{\textnormal{t}}+\zeta_{\textnormal{r}}$, as illustrated in Figure 5.

Although the $\lambda^{4}$-perturbation gives insight into the different contributions to the energy shift, perturbation beyond the order $\lambda^{4}$ also has a non-negligible contribution in the quasi-dispersive regime. Since RSWT requires considerably fewer commutators, we are able to compute the $\lambda^{6}$-perturbation, with only two iterations and 7 commutators (see Table 1). The $\lambda^{6}$-perturbation captures the location of the minimum more accurately, but still shows a false minimum close to the resonant regime, as shown in Figure 4b.

Apart from perturbation, we also apply NPAD to compute the interaction strength. We first define 4 Givens rotations with respect to the direct qubit-resonator coupling terms from the original Hamiltonian. The rotations are then applied sequentially to obtain the first effective Hamiltonian. Next, we apply another 4 rotations targeted at the two-photon couplings, such as the effective qubit-qubit coupling. The indices of those 8 rotations are listed in the first two columns of Table 2. These two steps are equivalent to the two iterations in RSWT. However, the recursive Givens rotations replace the BCH expansion, resulting in a much simpler calculation. Illustrated in Figure 4b, the approximation with those 8 rotations is as good as the $\lambda^{6}$-perturbation, but without the false minimum. Both capture the zero points very well compared to the numerical diagonalization, where the 4 lowest levels are included for each qubit and the resonator.

With those calculations, we can then investigate the effect of the high order corrections. We find that, for instance, $g_{q}$ shifts the zero point to the regime of smaller frequency detuning, corresponding to shrinking the half-circle in the numerical calculation in Figure 4a. In addition, for stronger coupling strength, the dip becomes narrower, which indicates a trade-off between the interaction strength and feasibility of qubit fabrication [51]. A detailed description of the effect of higher-order perturbation in the quasi-dispersive regime is presented in appendix D.

Overall, our investigation reveals different contributions to the ZZ interaction and provides tools to study the energy landscape in this quasi-dispersive regime. Because of the comparably smaller detuning, operations on this regime provide stronger interactions for entangling gates, and hence may achieve a better quantum speed limit for universal gate sets, i.e. without sacrificing local gates [10].

Following the previous examples, we here study superconducting qubits under an external cross-resonance drive. The cross-resonance interaction is activated by driving the control qubit with the frequency of the target qubit, which has been studied intensively and demonstrated in various experiments [52, 53, 54, 55, 56, 33]. In the two-qubit subspace, the dominant Hamiltonian term is written as a Pauli matrix $ZX$, which generates a CNOT gate up to single-qubit corrections. Therefore, ideally, only the population of the target qubit will change after the gate operation. The effective model is usually derived by block diagonalizing the non-qubit leakage levels as well as the population flip of the control qubit [7, 8, 57]. The coupling strength is then characterized by the coefficient of the $ZX$ Hamiltonian term.

The analytical block diagonalization of the Hamiltonian is only possible when neglecting all the non-qubit levels. Hence, perturbative expansion is often used, where the small parameter is defined as $\Omega/\Delta_{-}$, i.e., the ratio between the drive amplitude and the qubit-qubit detuning. However, to achieve fast gates, the qubit-qubit detuning is often designed to be small, ranging from 50 MHz to 200 MHz. Therefore, the perturbative diagonalization only works well for a weak drive.

In the rest of this subsection, we show that with only 4 two-by-two Givens rotations on the single-photon couplings, we can block-diagonalize the drive term and obtain an estimation of the coupling strength as good as the numerical result and far above the perturbative regime.

We start from the static Hamiltonian $H$ in Eq. 24 and define a driving Hamiltonian in the rotating frame

The full Hamiltonian is then written as $H+H_{\mathrm{d}}-H_{\mathrm{R}}$ where $H_{\mathrm{R}}=\omega_{d}(b_{1}^{\dagger}b_{1}+b_{2}^{\dagger}b_{2})$ with $\omega_{d}$ the driving frequency [7]. To compute the interaction strength, both the qubit-qubit effective interaction $g$ and the drive on the control qubit $\Omega$ need to be diagonalized. In particular, the second one can be as large as the energy gap and dominant in the unwanted couplings [57]. For simplicity, we assume $g$ is small and diagonalize it with a leading-order perturbation, discarding all terms smaller than $\mathcal{O}(g^{2})$. In this frame, one obtains a $ZX$ interaction that increases linearly with the drive strength [7]. This is equivalent to moving to the eigenbases of the idling qubits and allows us to focus on applying NPAD to the drive $H_{\mathrm{d}}$. The same method used in Section III.2 can be applied here to improve this approximation.

Targeting the dominant drive terms listed in the right column of Table 2, we construct 4 Givens rotations. The rotations are constructed with respect to the same Hamiltonians and then applied iteratively as separate unitaries. The obtained $ZX$ interaction strength reads

with $c_{j}=\cos(\theta_{j}/2)$, $s_{j}=\sin(\theta_{j}/2)$ and $\Delta_{-}=\omega_{1}-\omega_{2}$. The rotation angles are defined by the drive strength $\theta_{1}=\arctan(\frac{\Omega}{\Delta_{-}})$ and $\theta_{2}=\arctan(\frac{\sqrt{2}\Omega}{2\Delta_{-}+a_{1}})$.

This analytical coupling strength is plotted in Figure 6, compared with the perturbative expansions in Ref. [7] and numerical block-diagonalization. The result matches well with the numerical calculation, even when the ratio $\frac{\Omega}{\Delta_{-}}$ is approaching one. On the contrary, the perturbative expansion shows a large deviation as the driving power increases. The numerical block-diagonalization is implemented using the least action method [44, 7, 29]. To our surprise, although no least action condition is imposed on the Jacobi iteration, the method automatically follows this track and avoids unnecessary rotations. This suggests that the Jacobi iteration chooses an efficient path of block-diagonalization.

Notice that in the above example, no rotations are performed for levels beyond the second excited state because they are not directly coupled to the qubit subspace. In other parameter regimes, more couplings terms may become significant and need to be added to the diagonalization. For instance, the two-photon interaction between $\ket{0}$ and $\ket{2}$ of the control qubit will be dominant in the regime where $\Delta_{-}\approx-\alpha_{2}/2$ [8]. The fact that high precision can be achieved with only rotations on the single-photon couplings in this example also indicates that the dominant error of perturbation lies in the BCH expansion used in diagonalizing the strong single-photon couplings, rather than in higher levels or high-order interactions.