Accurate methods for the analysis of strong-drive effects in parametric gates

A possible realization of this three-mode system is shown in Fig. 1b) and consists of two fixed-frequency transmon qubits interacting via a capacitively shunted flux qubit whose two branches contain one and $N$ Josephson junctions, respectively You et al. (2007); Steffen et al. (2010); Yan et al. (2016). In a single-mode approximation, this generalized flux qubit plays the role of coupler mode and the parametric drive is realized by modulating the reduced external flux $\varphi_{\text{ext}}=2\pi\Phi_{\mathrm{ext}}/\Phi_{0}$, with $\Phi_{\mathrm{ext}}$ the flux threading the coupler’s loop and $\Phi_{0}$ the flux quantum. For certain values of the static external flux, the coupler has a positive anharmonicity, which is important in obtaining gates with a vanishing $ZZ$ interaction Ku et al. (2020); Zhao et al. (2020a); Xu and Ansari (2020); Zhao et al. (2020b). We stress that we use this specific circuit implementation for illustration purposes only, and that the methods presented here apply beyond the weakly-anharmonic regime.

Quantizing the circuit of Fig. 1c) using the standard approach Reynaud et al. (1997); Vool and Devoret (2017) yields the Hamiltonian (see Appendix B for a detailed derivation)

where the transmons and the coupler are described by

These expressions use pairs of canonically conjugate superconducting phase difference and Cooper pair number for the bare modes, $\left[\hat{\bm{\varphi}}_{j},\hat{\bm{n}}_{k}\right]=i\delta_{jk}$ for the mode indices $j,k=a,b,$ or $c$, and we set $\hbar=1$. The Josephson energies are denoted $E_{\text{J}a}$, $E_{\text{J}b}$ for the transmon modes, whereas $\beta E_{\text{J}c}$ is the Josephson energy of one of $N$ array junctions in the coupler, and $\alpha$ is a factor parametrizing the anisotropy between the two branches. The parameter $\beta$ is a renormalization of the superinductance due to disorder in the junction array and finite zero-point fluctuations (see Appendix B). Moreover, the parameter $\alpha$ accounts for a renormalization of the small junction energy due to hybridization with the modes in the junction array. Furthermore, $E_{\text{C}a}$, $E_{\text{C}b}$ and $E_{\text{C}c}$ are charging energies. In the transmon regime, $E_{\text{J}a}/E_{\text{C}a}$ and $E_{\text{J}b}/E_{\text{C}b}\gtrsim 50$ Koch et al. (2007).

The coupler loop is threaded by an external flux $\varphi_{\text{ext}}$ which can be modulated in time with a modulation amplitude $\delta\varphi$, taken to be small compared to the flux quantum

As discussed by You et al. (2019), quantization of the coupler loop under time-dependent flux imposes that the external flux be included in both branches of the potential energy in $\hat{H}_{c}(t)$, with weighting factors $\mu_{\alpha,\beta}$ determined by the capacitive energies of the two branches (see Appendix B for a detailed derivation). This subtlety is important, as the details of the flux modulation determine the parametric interactions between the two qubit modes.

Finally, the three bare modes interact through linear terms induced by the capacitive coupling

The introduction of normal modes will eliminate this linear coupling Hamiltonian.

In this subsection, we introduce a simple model which captures the essential qualitative features of the Hamiltonian of Eq. 1. Our toy model consists of three linearly coupled Kerr-nonlinear oscillators and has the form given in Eq. 1 now with

Comparing to the full-circuit model, note that $\bm{\omega}_{a(b)}\approx\sqrt{8E_{\text{J}a(b)}E_{\text{C}a(b)}}-E_{\text{C}a(b)}$ whereas the anharmonicities of the transmon qubits are negative and amount to $\bm{\alpha}_{a(b)}\approx-E_{\text{C}a(b)}$. In an experimental implementation, the parameters defining the coupler—the anharmonicity $\bm{\alpha}_{c}$ and the frequency $\bm{\omega}_{c}(t)$—can be varied by applying a time-dependent external flux to activate a chosen gate.

The parametric drive resulting from the flux modulation of Eq. 3 is modeled by a modulation of the coupler frequency at a frequency $\omega_{\text{d}}$

In a more detailed analysis of the coupler (see Sec. IV), we take into account the time dependence of the anharmonicity $\bm{\alpha}_{c}$, but we choose to neglect it in this toy model.

Note that we have reduced the complexity of the problem in a few ways: We have truncated the Josephson expansion to include only quartic terms. All photon number nonconserving terms have been dropped. Higher harmonics of the drive of Eq. 6 are neglected, and we have not considered the ac-Stark shifts of the various coupling constants. All of these refinements are taken into account in the analysis of the full circuit Hamiltonian in Sec. IV. Thus the toy model is significantly simpler than the full circuit theory, but nonetheless still contains the necessary ingredients that allow us to illustrate the general method introduced in this paper.

As usual, our starting point is a decomposition of the system Hamiltonian into an unperturbed, exactly solvable part, and a perturbation:

Here, we have introduced the dimensionless power-counting parameter $\lambda$ to keep track of the order in perturbation theory, to be set at the end of the calculation to unity, $\lambda\to 1$. Now we move to the interaction picture with respect to $\hat{H}^{(0)}$. Letting $\hat{U}_{0}(t)=\mathcal{T}e^{-i\int_{0}^{t}dt^{\prime}\hat{H}^{(0)}(t^{\prime})}$, where $\mathcal{T}$ is the time-ordering operator, we find for the interaction-picture Floquet Hamiltonian

In the above we only assume that the unperturbed time-evolution operator $\hat{U}_{0}(t)$ is known. Equation 8 can be seen as a unitary transformation between two Floquet Hamiltonians Sambe (1973). Thus, the Floquet quasienergies corresponding to $\lambda\hat{H}_{I}^{(1)}(t)-i\partial_{t}$ must be identical to those of $\hat{H}^{(0)}+\lambda\hat{H}^{(1)}(t)-i\partial_{t}$, while the eigenstates are related by $\hat{U}_{0}(t)$.

In an iterative Schrieffer-Wolff approach, we treat the operator $\lambda\hat{H}^{(1)}{(t)}$ as a small perturbation from which we derive corrections to the known Floquet quasienergies of $\hat{H}^{(0)}$ Malekakhlagh et al. (2020b); Petrescu et al. (2020). To this end, we consider a unitary transformation on the interaction-picture Floquet Hamiltonian $\hat{H}_{I}(t)-i\partial_{t}\equiv\lambda\hat{H}_{I}^{(1)}(t)-i\partial_{t}$, and the corresponding Baker-Campbell-Hausdorff (BCH) expansion in powers of the generator of this unitary, that is

This equation defines the effective Hamiltonian, whose spectrum is equal (up to a desired precision in $\lambda$) to that of the original driven theory. The generator $\hat{G}_{\text{I}}(t)$ can be solved for iteratively in powers of $\lambda$ (see Appendix A), which allows us to perform the rotating-wave approximation order by order

where the terms on the right-hand side are defined below.

To obtain a lowest-order term of the effective Hamiltonian, $\lambda\overline{\hat{H}}_{I}^{(1)}$, we separate the interaction picture Hamiltonian into oscillatory and non-oscillatory terms with the notation

where we define the constant part of a time-dependent operator $\hat{O}(t)$ by Mirrahimi and Rouchon (2015)

whereas the oscillatory part of the operator is

Since the time-averaging operation removes all terms that are oscillatory in time, $\overline{\hat{H}}_{I}^{(1)}$ is the first-order RWA Hamiltonian Mirrahimi and Rouchon (2015), whereas $\widetilde{\hat{H}}_{I}^{(1)}$ is canceled by an appropriate choice of the corresponding term at order $\lambda$ in the generator.

One can iterate this procedure at every order, collecting terms that are oscillatory and then canceling them. The second-order RWA Hamiltonian (for a derivation, see Appendix A) reads

This form becomes analogous to the second term in the Magnus expansion Goldman and Dalibard (2014); Mirrahimi and Rouchon (2015) when the perturbation has a vanishing mean, i.e. $\overline{\hat{H}}_{I}^{(1)}=0$.

Expressing the toy-model Hamiltonian as the sum of static quadratic terms, $\hat{H}^{(0)}$, and of time-dependent and Kerr terms, $\hat{H}^{(1)}(t)$, the first step in deriving parametrically activated interactions between the transmon modes is to diagonalize the former, which we write as

This diagonalization is achieved with an orthonormal transformation $\hat{\bm{\alpha}}=\sum_{\beta=a,b,c}u_{\alpha\beta}\hat{\beta}$, for $\alpha=a,b,c$, and which is chosen such that $\hat{H}^{(0)}$ takes the form

where $\hat{a}$, $\hat{b}$ and $\hat{c}$ are the normal modes and $\omega_{a,b,c}$ the corresponding mode frequencies. The $u_{\alpha\beta}$ are hybridization coefficients encoding the connectivity of the three modes through the capacitive couplings $\hat{H}_{g}$ entering in $\hat{H}^{(0)}$. In this normal-mode basis, the remainder of the Hamiltonian reads

The expression above illustrates that coupling between the normal modes arises from the nonlinearity and the parametric drive.

Our choice of unperturbed Hamiltonian $\hat{H}^{(0)}$ and perturbation $\lambda\hat{H}^{(1)}$ in Eqs. 23 and LABEL:Eq:H1Toy, respectively, is guided by black-box quantization Nigg et al. (2012): the unperturbed Hamiltonian is linear and diagonal in the normal-mode basis, whereas the perturbation consists of Kerr-nonlinear terms, on the one hand, and quadratics appearing from the normal-mode expansion of the parametric drive, on the other hand. As we show below, while better choices are possible (see Sec. III.3), this choice leads to a simple and intuitive form for the effective Hamiltonian.

As an example of the many common types of interactions that can be activated by a parametric drive Blais et al. (2021), an $i$SWAP interaction between the transmon modes arises if we set the modulation to be at the frequency difference between the two transmon modes

This choice yields the first-order RWA Hamiltonian

The first row of this equation contains the $i$SWAP interaction of amplitude $J^{(1)}_{ab}$. The second row contains the mode anharmonicities, and the third row contains cross-Kerr interactions, the first of which is the $ZZ$ term.

The coupling constants in the above effective Hamiltonian result from the normal-mode transformation of the quadratic part of the toy model and take the form

for all $j,k=a,b,c$, and $j\neq k$. In practice, one wants to maximize $J_{ab}^{(1)}$ to obtain a fast gate, while minimizing the cross-Kerr interactions $\chi^{(1)}_{jk}$ to avoid the accumulation of coherent errors. Cross-Kerr interactions are a source of infidelity for a $i$SWAP-type gate, as well as in other gate implementations Ku et al. (2020); Sung et al. (2020); Noguchi et al. (2020); Kandala et al. (2020); Zhao et al. (2020a); Krinner et al. (2020b). In the first-order RWA Hamiltonian, to cancel the cross-Kerr interaction between the two transmons, we use a coupler with a positive anharmonicity Zhao et al. (2020b) $\bm{\alpha}_{c}>0$, together with $\bm{\alpha}_{a},\bm{\alpha}_{b}<0$, which is distinct from using qubits of opposite anharmonicities Ku et al. (2020); Xu and Ansari (2020), or couplers with negative anharmonicity Zhao et al. (2020a); Sete et al. (2021). Equation 27 forms the basis for the optimization of the gate parameters. Before pursuing this further, we first derive important corrections to the gate Hamiltonian from the oscillatory part of the Hamiltonian, $\widetilde{\hat{H}}_{I}^{(1)}(t)$. Finally, note that the first-order term $\chi^{(1)}_{jk}$ is only a static, i.e. $\delta$-independent, cross-Kerr interaction.

At second order in perturbation theory, there is no correction to the $i$SWAP gate frequency $J_{ab}^{(2)}=0$. In the regimes of interest, where the coupler frequency is close enough to the qubit frequencies for the interaction between the coupler and the qubits to be non-negligible, the dominant contribution to the second-order RWA correction to the cross-Kerr interaction $\chi_{ab}^{(2)}$ is

The full expression for $\chi_{ab}^{(2)}$ can be found in Appendix C. Inspecting the hybridization coefficients $u_{\alpha\beta}$ and the denominators in Eq. 28, we deduce that the second-order correction to the static cross-Kerr interaction, corresponding to the first term in Eq. 28, arises from a virtual two-photon excitation of the coupler (generated by the commutator $[\hat{a}\hat{b}\hat{c}^{\dagger 2},\hat{a}^{\dagger}\hat{b}^{\dagger}\hat{c}^{2}]$). This correction would not be present in a two-level approximation McKay et al. (2016). On the other hand, the second term in Eq. 28 is the lowest-order contribution to the dynamical cross-Kerr interaction.

As mentioned in the previous subsection, other choices for $\hat{H}^{(0)}$ and $\lambda\hat{H}^{(1)}$ are possible which give better accuracy in comparisons with exact Floquet numerics. In this subsection, we take the unperturbed Hamiltonian $\hat{H}^{(0)}(t)$ to consist of the Fock-space diagonal part of $\hat{H}(t)$, namely

where the quartic couplings in the last two rows are exactly those defined in Eq. 27, with the superscript now changed from $1$ to $0$ to reflect their presence in the unperturbed Hamiltonian.

We expect this starting point, Eq. 29, to lead to more precise results, because of two reasons. Firstly, the perturbation $\hat{H}^{(1)}$ is now off-diagonal in Fock-space. The effects of the anharmonicities of the modes are now included at the level of $\hat{H}^{(0)}$, and, in particular, we expect a dressing of contributions corresponding to two-photon excitations, such as Eq. 28. Secondly, due to the second row of Eq. 29, we can derive the effect of harmonics of the drive frequency through a Fourier expansion of the time-evolution operator $\hat{U}_{0}(t)$, defined below.

Following the steps of the previous subsections, we evaluate the time evolution operator with respect to the unperturbed Hamiltonian $\hat{H}^{(0)}(t)$, that is

with $[\hat{H}^{(0)}(t),\hat{H}^{(0)}(t^{\prime})]=0$. The time-dependence in the exponent is handled via the Jacobi-Anger identity Abramowitz and Stegun (1964)

for any operator $\hat{O}$, where $J_{n}(z)$ is the $n^{\textit{th}}$ Bessel function of the first kind Abramowitz and Stegun (1964). This expansion allows us to keep track of all harmonics of the drive. In practice, since the modulation amplitude is small, $\delta/\omega_{\text{d}}\ll 1$, only a few terms will be necessary.

We obtain to first order (and truncating after the first Bessel function)

which agrees, up to linear terms in $\delta$, with the expression in Eq. 27. On the other hand, the cross-Kerr interaction at this order is vanishing $\chi_{ab}^{(1)}=0$.

To second order in perturbation theory, the dominant contributions to $J_{ab}^{(2)}$ are

We do not reproduce here the full form containing 20 terms. The second-order contribution $\chi_{ab}^{(2)}$ is non-vanishing, and contains approximately 450 terms in expanded form. Despite the complexity of these full expressions, they are easy to derive and manipulate with symbolic computation tools Žitko (2011). Focusing on the static cross-Kerr interaction, i.e. in the $\delta\to 0$ limit, the dominant correction resulting from the above changes amounts to our previous Eq. 28, but replacing the denominator of the first term of that expression by a form that faithfully includes the contribution from the anharmonicities, as expected in the case of a virtual two-photon excitation of the coupler mode. That is, approximately

The full expression for the corrected static cross-Kerr interaction can be found in Appendix C.

In Fig. 2 we compare analytical results to numerical results obtained from exact diagonalization (at $\delta=0$), or a solution of the Floquet eigenspectrum at $\delta\neq 0$ (see Sec. V). We find that the agreement between numerics and analytics is excellent for the gate interaction strength, as well as for the static $\delta=0$ cross-Kerr interaction. However, we find that second-order perturbation theory is insufficient to reproduce the effects of the drive on the cross-Kerr interaction, even for modest drive amplitudes. We expect that higher-order perturbation theory should correctly capture the drive-amplitude dependence of the anharmonicities, but these contributions have been inaccessible in our study due to the large memory demands of the computer algebra manipulations.

To define the bare modes, we begin with the full circuit model Hamiltonian of Eqs. 1 and 2. We normal-order expand the Josephson cosine potentials in this Hamiltonian over a set of creation and annihilation operators, which we define as follows

with analogous equations for modes $\hat{\bm{b}}$ and $\hat{\bm{c}}$. The coefficients $\eta_{a,b,c}$ are chosen such that terms proportional to $\hat{\bm{a}}^{2},\hat{\bm{b}}^{2}$, and $\hat{\bm{c}}^{2}$ vanish in the time-averaged Hamiltonian. This amounts to three transcendental equations:

for $j=a,b,c$, where we have defined the form factors

Two remarks are in order: First, in the transmon limit $\mathcal{F}(\eta_{a(b)})\approx 1$ where we recover the usual expression $\eta_{a(b)}\approx\sqrt{8E_{\text{C}a/b}/E_{\text{J}a/b}}$. Second, the parameter $\eta_{c}$ depends on the parametric drive amplitude $\delta$, which indicates that the mode $c$ impedance is drive-dependent. This has important consequences for the precision of the calculation of coupling constants dressed by the parametric drives. In particular it allows us to capture the ac-Stark shift of the coupler mode at the lowest order in perturbation theory. In what follows, sine and cosine functions of the phase are normal-order expanded according to Eq. 83 of Appendix D. In turn, trigonometric functions of the flux modulation are expanded in Jacobi-Anger series over the harmonics of the frequency of the drive.

Using the above definitions, the transmon Hamiltonian $\hat{H}_{a}$ takes the familiar form

The second term on the right-hand side contains the nonlinear part of the Josephson potential, i.e. the inductive part is subtracted. Up to quartic order, $\hat{H}_{a}$ takes the form

The first row of this expression is a Kerr oscillator Hamiltonian as in the toy model of Sec. II, whereas the second row contains corrections from quartic counter-rotating terms. Here, we have introduced the mode frequency and anharmonicity which take the forms Koch et al. (2007)

Note that for the approximate equality in the first row we have used a Taylor expansion of Eq. 36 for $\eta_{a}$. The equations for mode $\hat{\bm{b}}$ are identical from the above with a change of subscripts and operators $a\to b$.

The coupler Hamiltonian differs from that of the transmon modes in two fundamental ways: It breaks parity symmetry due to the external flux, and it is time-dependent. Following Eqs. 13 and 12, we write this time-dependent Hamiltonian as

The creation and annihilation operators of the coupler mode can then be defined by extracting the quadratic part of the time-averaged coupler Hamiltonian. Using Eq. 35 where $a\rightarrow c$ together with

and a similar relation for the second branch of the coupler [see Eq. 2], we find in analogy to Eq. 38 for the Hamiltonian of the transmon mode

Crucially, in this first-order rotating-wave approximation of the parametric drive, the Josephson energy is renormalized by the factor $J_{0}(\mu_{\alpha,\beta}\delta)$, see also Verney et al. (2019). We interpret this as an effective reduction of the Josephson potential barrier, and consequently an increase of phase fluctuations, in the presence of drives. Moreover, the presence of the non-zero external flux results in the parity breaking sine terms in Eq. 43.

The second term of $\hat{H}_{c}(t)$ in Eq. 41, the oscillatory part, take the form

which can be expanded in a Jacobi-Anger series in harmonics oscillating at the frequency $n\omega_{\text{d}}$, where $n$ is an integer.

As above, the next step is to expand the coupler Hamiltonian up to quartic terms in the creation and annihilation operators. In contrast to the transmon Hamiltonian of Eq. 39, parity breaking leads to the appearance of monomials of odd order. The non-oscillatory part is

The first row of the above expression takes the form of the coupler Hamiltonian in the approximation of the toy model of Sec. II, while the remaining rows contain counter-rotating terms to quartic order. Here, the parametric drive-dependent mode frequency and anharmonicity read

while the prefactors of the counter-rotating terms are

The contribution from the oscillatory part $\widetilde{\hat{H}}_{c}(t)$ is too lengthy to be reproduced here, and is given up to the second harmonic of the parametric modulation frequency $\omega_{\text{d}}$ in Table 2 of Appendix D.

Finally, the last term of the full circuit Hamiltonian to consider is the linear interaction $\hat{H}_{g}$ induced by the capacitive coupling. Using Eq. 35, this Hamiltonian takes the form

where the ellipsis represents two more terms corresponding to the cyclic permutations of the mode indices.

The Hamiltonian specified by Eqs. (39) and its equivalent for the $b$ transmon mode, Eq. 45, the terms summarized in Table 2 of Appendix D, and Eq. 48, form the basis of both the normal-mode analysis, and of the full-circuit numerical simulation performed in Sec. V.

Building upon the above results, we now follow the procedure developed with the toy model in Sec. III.2 to obtain effective gate Hamiltonians under parametric modulations. To do so, we first collect under $\hat{H}^{(0)}$ the time-independent quadratic terms derived above. We then eliminate the linear coupling $\hat{H}_{g}$ of Eq. 48 from $\hat{H}^{(0)}$ through a normal-mode transformation

The linear transformation is determined by a set of 18 hybridization coefficients that relate bare mode coordinates to normal mode coordinates (see Appendix E for the procedure to compute these coefficients)

for $\alpha=a,b,c$. We stress that the hybridization coefficients $u_{\alpha\beta}$, $v_{\alpha\beta}$ depend on the amplitude of the parametric drive. As a result, drive effects such as the ac-Stark shift of the nonlinear oscillators are accounted for already at the level of the normal-mode decomposition of the circuit Hamiltonian.

In analogy with our treatment in Sec. II of the toy model, we have taken $\hat{H}^{(0)}$ to be the unperturbed Hamiltonian with respect to which the interaction picture is defined. Primarily in order to keep the expressions more concise, we opt to neglect the corrections analyzed in Sec. III.3. With $\hat{H}^{(0)}$ the unperturbed Hamiltonian, the remaining interaction terms are $\lambda\hat{H}^{(1)}(t)\equiv\hat{H}-\hat{H}^{(0)}$. Expressing these in the interaction picture with respect to $\hat{U}_{0}=e^{-i\hat{H}^{(0)}t}$ as in Eq. 8, we find