Proposal for entangling gates on fluxonium qubits via a two-photon transition

We propose a family of microwave-activated entangling gates on two capacitively coupled fluxonium qubits. A microwave pulse applied to either qubit at a frequency near the half-frequency of the $|00\rangle - |11\rangle$ transition induces two-photon Rabi oscillations with a negligible leakage outside the computational subspace, owing to the strong anharmonicity of fluxoniums. By adjusting the drive frequency, amplitude, and duration, we obtain the gate family that is locally equivalent to the fermionic-simulation gates such as $\sqrt{\rm SWAP}$-like and controlled-phase gates. The gate error can be tuned below $10^{-4}$ for a pulse duration under 100 ns without excessive circuit parameter matching. Given that the fluxonium coherence time can exceed 1 ms, our gate scheme is promising for large-scale quantum processors.


I. INTRODUCTION
A programmable quantum computer requires a very low error rate for the two-qubit gate operations, both for quantum error-correction schemes to work [1,2], and for extending the depth of quantum circuits during calculations on noisy intermediate-scale quantum processors [3]. In the superconducting circuits platform [4][5][6], major results were obtained using transmon qubits [7,8], which are much closer to weakly anharmonic oscillators than to two-level systems. Although simplicity and robustness of transmons facilitated the creation of processors with dozens of qubits [9,10], the weak anharmonicity and finite coherence time have been major factors limiting gate errors. These challenges exist in both major families of two-qubit gates realized with these qubits: fluxtunable [11][12][13][14][15][16] and microwave-activated [17][18][19][20][21][22][23] twoqubit gates, where the gate speed is bounded by the anharmonicity.
Fluxonium qubits [24] are architecturally similar devices to transmons but they have a much stronger anharmonicity and considerably longer coherence times [25,26]. Theoretical proposals to realize microwave-activated two-qubit gates with fluxoniums and heavy fluxoniums have previously been based on driving transitions outside of the computational subspace [27,28]. Recent experiments demonstrated fast two-qubit gates on fluxoniums activated by driving close to such transitions [29,30]. Because these noncomputational transitions generally have shorter lifetimes than the computational ones, such a gate scheme is exposed to additional error channels. A flux-tunable entangling gate with fluxonium qubits has also been recently reported [31], but it is subject to extra dephasing errors when a qubit is moved away from its flux sweet spot.
In this work we consider a gate that keeps the state entirely in the computational subspace with qubits parked at their sweet spots, hence benefiting in full from the long coherence of fluxonium qubits. The entangling gate presented here is accomplished by a high-power microwave drive at half the frequency of the |00 − |11 transition, which induces two-photon transitions between |00 and |11 and activates a coherent mixing in the {|00 , |11 } subspace. With weak interaction in the computational subspace, the gate would normally be slow. However, the strong anharmonicity of fluxoniums makes it possible to perform fast gate operations by increasing the drive amplitude without generating leakage to noncomputational levels. For realistic fluxonium parameters, we demonstrate that a 50-ns-long gate with π/2 mixing angle in the {|00 , |11 } subspace can be realized with the leakage error below 10 −4 and total gate error below 10 −3 without using advanced pulse shaping for presently achievable coherence times [25,26]. For longer gates, about 100 ns long, the coherent gate error can be reduced below 10 −4 . With decoherence effects accounted for, the 10 −4 threshold requires some improvement of the best existing lifetimes [25,26] and should be possible in next-generation devices.
The entangling power [32,33] and the local equivalence class of the proposed gates depend on the mixing angle in the {|00 , |11 } subspace and on the magnitude of the effective ZZ coupling. The ZZ coupling originates from both the static level repulsion and an induced ac-Stark shift due to the large drive amplitude. This term modifies entangling properties of the gate. Interestingly, we show that for a half rotation in the {|00 , |11 } subspace, the entangling power of the gate is independent of the contribution due to the total ZZ coupling. In fact, the family of gates with a half-rotation contains gates locally equivalent to √ iSWAP and √ SWAP. Entangling gates activated by two-photon processes were proposed and implemented in trapped-ion systems [34,35]. In superconducting systems, a two-photon gate based on driving the |00 − |11 transition was demonstrated experimentally with transmon qubits [20]. This gate required frequency matching between |0 − |1 transition of one transmon and |1 − |2 transition of another transmon to increase the two-photon Rabi frequency by increasing the hybridization of the |11 state with one of the noncomputational states. In the case of fluxoniums, we can speed up the two-photon Rabi oscillations by increasing the drive amplitude. This trick does not cause leakage outside the computational subspace because of the strong anharmonicity of fluxoniums. Thus, the scheme presented here does not require frequency matching. In general, our gate benefits from higher frequencies of qubit transitions, which lead to a stronger hybridization of states |01 and |10 and, therefore, to a higher two-photon Rabi frequency. From practical considerations, the suggested qubit frequency range is around 0.5 − 1 GHz.
While a single entangling gate combined with individual qubit controls is sufficient to generate a universal set of logical operations, some algorithms may compile more efficiently with a larger two-qubit gate set, especially if this set is hardware efficient [36,37]. For a given algorithm executed on a noisy processor, the maximal depth of a quantum circuit depends on a particular set of gates implemented on the hardware level. The gates based on mixing of |00 and |11 discussed here are locally equivalent to and can be easily mapped by single-qubit X rotations into the operations in the {|01 , |10 } subspace. The family of the proposed gates is equivalent to a complete set of fSim gates for the fermionic-simulation problem [14,38].
The outline of the paper is as follows. In Sec. II, we introduce our model and discuss relevant spectral properties of fluxonium circuits. In Sec. III, we consider coherent two-photon transitions between |00 and |11 twoqubit states and calculate their rate both analytically and numerically. In Sec. IV, we analyze two-qubit gates realized via the two-photon |00 −|11 transition. We discuss local equivalence classes of such gates and their entangling power and simulate the gate error. We conclude in Sec. V.

II. CAPACITIVELY COUPLED FLUXONIUM QUBITS
The circuit diagram of two coupled fluxoniums, labeled as A and B, is shown schematically in Fig. 1(a). We model this system by the Hamiltonian wherê describes individual qubits (α = A, B) [24]. The capacitive interaction between qubits is given bŷ and the coupling to an external microwave drive of frequency ω d and phase γ d -bŷ (a) Circuit diagram of two capacitively coupled fluxonium qubits. (b) Energy levels (thick horizontal lines) of the coupled system are separated into two groups based on whether or not single-qubit indices have the same parity (black solid and magenta dashed lines). Transitions between different parity groups can be performed with a single microwave photon, while transitions within the same group require higher-order multiphoton processes. In the perturbative regime, the coherent two-photon transition between |00 and |11 is dominated by contributions from two virtual states generated by states |10 and |01 , as shown by the wide arrows. These contributions interfere destructively and exactly cancel each other when the coupling vanishes.
Here = h/2π is the reduced Planck constant, f (t) is the time-dependent field envelope, and η A and η B describe the coupling of each qubit to the driving field.
In these equations, the canonical variables are the dimensionless fluxφ α and charge (the number of Cooper pairs)n α , which satisfy the commutation relations [φ α ,n α ] = iδ αα . The kinetic-energy term in Eq. (2) is determined by the charging energy E C,α = e 2 /2C α , where (−e) is the electron charge and C α is the total capacitance of the circuit α. The inductive energy is E L,α = ( /2e) 2 /L α , where L α is the effective linear inductance of a long chain of Josephson junctions, which is a hallmark of the fluxonium. This superinductance is shunted by a small junction, associated with the Josephson energy E J,α . The final term in Eq. (2) depends on φ ext,α = (2e/ )Φ ext,α , where Φ ext,α is the magnetic flux threading the loop formed by the small junction and superinductance. In the limit of C M C A , C B , where C M is the mutual capacitance, the interaction strength in Eq. (3) is given by [27,39].
Below, we assume that fluxonium circuits are at their half-flux-quantum sweet spots defined by φ ext,α = π, where the circuits are first-order insensitive to the external flux noise [40]. We label eigenstates of Hamiltonian (2) as |0 α , |1 α , |2 α , . . . in increasing order of the corresponding eigenenergies E α 0 ≤ E α 1 ≤ . . .. The first two levels of each circuit α form a qubit with transition frequency ω α 01 , where we define single-fluxonium frequencies as ω α kl = E α l − E α k . The qubit transition |0 α − |1 α can display an exceptionally long coherence time exceeding 500 µs [25], which makes it an attractive choice for quantum-information storage. The qubit transition frequency ω α 01 /2π is typically in the 100 MHz -1 GHz range, which is much lower than conventional values in superconducting qubits. In addition, the charge matrix element n α 01 with the notation n α kl = | k α |n α |l α | is suppressed at low frequencies. At the same time, the transition |1 α − |2 α has properties similar to those of the transmon with a typical frequency of several gigahertz [40][41][42]. Because of the potential-energy symmetry at φ ext,α = π, the matrix elements ofn α display parity selection rules, e.g., n α 02 = n α 13 = 0 [27,41,43]. However, n α 03 is not suppressed and can be of the order of n α 12 [41,43]. We label interacting (dressed) two-qubit eigenstates of the Hamiltonian (1) atĤ drive = 0 as |kl implying adiabatic connection to the noninteracting tensor-product states |kl 0 = |k A |l B . The frequency of the two-qubit transition |kl −|k l is denoted as ω kl−k l = E k l −E kl , where E kl is the eigenenergy of |kl . The two-qubit computational subspace {|00 , |01 , |10 , |11 } is well separated from higher levels as illustrated schematically in Fig. 1(b). Two-qubit levels can be divided into two parity groups depending on whether k + l is even or odd, which is shown with solid and dashed lines in Fig. 1(b). Because of the parity selection rules for the charge operatorsn A andn B , the interaction term (3) mixes levels within the same parity group only, while the matrix elements of the drive term (4) are nonzero only between levels belonging to different parity groups [27].

III. TWO-PHOTON RABI OSCILLATIONS
In this section, we consider a continuous drive of the two fluxoniums with a constant amplitude f in Eq. (4) and a drive frequency about half the frequency of the |00 -|11 transition. Understanding this process is essential for the two-qubit gate discussed in Sec. IV.

A. Rabi frequency
Even though two fluxonium excitations cannot be created by a single microwave photon, i.e., 00|Ĥ drive |11 = 0 because of the parity selection rules, the microwave drive can still induce the |00 − |11 transition exchanging two fluxonium excitations with pairs of drive photons. The transition amplitude calculated to the leading (second) order in f can be understood as having contributions from cascaded sequential single-photon transitions via intermediate real states such as |01 and |10 and from coherent two-photon processes via intermediate virtual states [44]. Under certain conditions discussed below, the excitation probabilities of states |01 and |10 can remain low, while the system state oscillates between |00 and |11 with high visibility. Below, we apply a perturbation theory to estimate these probabilities and the frequency of the two-photon Rabi oscillations between |00 and |11 . We approximate fluxoniums as two-level systems, which is reasonable given their strong anharmonicity. Nonperturbative effects and higher fluxonium levels are accounted for in numerical analysis of Sec. III C.
The resonant Rabi frequency for a single-photon transition |kl −|k l such as |00 −|10 is given by the matrix elements of the drive [see Eq. (4)]: where N kl,k l = kl|(η AnA + η BnB )|k l .
Without qubit interactions, the eigenstates of the system are product states of each qubit, |kl 0 = |k A ⊗ |l B , and the Rabi frequencies (5a) reduce to the single-qubit Rabi frequencies for the |0 α − |1 α transitions Here we consider microwave drives with frequencies close toω which reduces to the average of qubit frequencies (ω A 01 + ω B 01 )/2 when J C = 0. When driving with frequency ω d =ω, the probability of single-fluxonium transitions is bounded from above by the contrast of Rabi oscillations in two-level systems. The latter is given by where ∆ AB = |ω A 01 − ω B 01 | is the detuning between qubit frequencies. Taking P α 1, we obtain conditions on the drive amplitudes, which combined with ω d =ω imply that correlated twoqubit oscillations will dominate the dynamics of the system over independent single-qubit excitations. We now consider the transition between |00 and |11 , which is activated by two-photon processes when the drive frequency isω. In this case, the system exhibits full oscillations of its probability being in one of the states, and the frequency Ω of such two-photon Rabi oscillations depends on the matrix element of the two-photon drive between |00 and |11 . We apply second-order perturbation theory together with the rotating-wave approximation (RWA) to obtain Ω = |Ω 00−01 Ω 01−11 − Ω 00−10 Ω 10−11 | where ∆ AB = |ω 00−01 − ω 00−10 | , which differs from ∆ AB by a correction quadratic in J C . This equation is reminiscent of a similar result derived for trapped ions excited with a bichromatic laser [34,35]. By choosing ω d =ω in the derivation of Eq. (10), we neglected the shift of the qubit frequencies due to the Stark shifts, which are quadratic in the drive amplitude f . Equation (10) describes destructive interference between the two contributions corresponding to two paths via virtual states generated by states |01 and |10 , indicated by arrows in Fig. 1(b). Without interaction, J C = 0, we observe that Rabi frequencies Ω kl−k l for single-photon transitions reduce to single-qubit Rabi frequencies Ω A,0 and Ω B,0 , and Eq. (10) yields Ω = 0. This result emphasizes that entanglement is impossible without qubit interactions.

B. Interaction effects
Let us now calculate the two-photon Rabi frequency to the first nonvanishing order in | is typically small for fluxonium qubits. Correction to the denominator of Eq. (10) is quadratic in J C , while corrections to matrix elements are linear. Thus, we find that Ω is finite because of the hybridization of |01 with |10 and of |00 with |11 . The interaction-dressed eigenstates in the first pair of states have the form (12) The mixing amplitude for the pair of states |00 and |11 has a form similar to Eq. (12), but with 2ω in the denominator instead of ∆ AB . Thus, the mixing of states |00 and |11 by interaction is reduced by ∆ AB /2ω compared to that of states |01 and |10 . This factor is not necessarily small in fluxonium qubits. Nevertheless, we ignore it for now to focus on the main principles. We find the following expression for the two-photon Rabi frequency: The two-photon rate increases with hybridization in the computational subspace. One natural way to increase the latter is to reduce the detuning ∆ AB , which, however, enhances the magnitude of spurious single-photon excitations when driving at ω d ≈ω; see Eq. (8). To ensure the predominance of the two-photon process over single-photon excitations, we fix the dimensionless amplitude Ω α,0 /∆ AB rather than f since the dimensionless amplitude fully determines the relative importance of single-photon excitations, see Eq. (8). We conclude that reducing ∆ AB is not practical for increasing the rate of the two-photon |00 − |11 transition. Equation (13) suggests that it is beneficial to have larger values of single-qubit charge matrix elements n α 01 . Because n α 01 = ω α 01 | 0|φ α |1 α |/8E C,α [27], this condition means that having higher qubit frequencies is beneficial to make the two-photon Rabi oscillations faster.
In addition to inducing coherent two-photon oscillations, a strong drive at ω d =ω induces ZZ interactions via the ac-Stark effect (static ZZ interaction is absent in two-level models [45,46]). In particular, we evaluate the relative phase accumulated in the computational subspace for a constant drive during time t: where ∆E kl is the energy shift of level |kl due to the ac-Stark effect. Using second-order perturbation theory for the energy shifts due to the drive combined with firstorder corrections in the coupling rate, Eq. (12), we find that Ωt .
In general, this phase accumulation is not negligible during a full Rabi period t = 2π/ Ω. Note that Eq. (15), obtained to first order in interaction J c , vanishes if the drive is applied to only one qubit, so either η A = 0 or η B = 0 in Eq. (4). A more rigorous analytic treatment of the two-photon process for fluxonium qubits is unnecessarily cumbersome. Below, we present a detailed numerical analysis of the two-photon Rabi oscillations in a system of two fluxoniums.

C. Numerical simulations of Rabi oscillations
Based on Eq. (13), we choose the single-qubit parameters for all numerical simulations in this paper as shown in Table I. With these values, both main qubit transition frequencies are relatively high compared to usual fluxonium qubits, which is accompanied by larger values of 0-1 charge matrix elements. In simulations, we first diagonalize Hamiltonians of single fluxonium circuits (2), and then work with the interacting system taking five lowest levels in each fluxonium. We use the full Hamiltonian (1) in the laboratory frame and, therefore, go beyond the RWA in addition to going beyond the perturbation theory and two-level approximation used in Eq. (10). For simplicity, we only consider the case η A = η B = 1 here.
We plot the Rabi frequency for the |00 −|11 transition as a function of the dimensionless drive amplitude    (10) except that multilevel fluxonium qubits were used to compute energies and matrix elements, and symbols show results of numerical simulations, which were performed as follows. For given λ and J C , we chose ω d that maximizes the contrast between the minimum of P 00→00 (t) and maximum of P 00→11 (t), where P kl→k l (t) is the population of state |k l at time t for the initial state |kl at t = 0. At each ω d , we calculated these probabilities via time-domain simulations at time t ≥ 0. In the drive term (4), we used a pulse with a Gaussian rising edge at 0 < t < t rise = 25 ns, which is given by where σ = t rise /2, and with the amplitude of the flat part (t > t rise ) determined by λ via Eq. (16). An example of such time-domain simulations of the occupation probabilities for ω d that maximizes the contrast is shown in Fig. 2(c) for λ = 0.5 and J C /h = 300 MHz. The optimal contrast in Fig. 2(c) is approximately 80%, which is different from 100% because of a finite t rise in Eq. (17) and because of leaking single-photon transitions such as |00 − |01 . This contrast is at least 75% for all the pairs of λ and J C discussed in Fig. 2 and is close to 100% for λ 0.1. The observed period of oscillations in Fig. 2(c) is 147 ns, which corresponds to a Rabi frequency of 6.8 MHz in agreement with the circles in Figs. 2(a) and 2(b). In a two-level system with qubit A parameters, λ = 0.5 corresponds to the 50% contrast of single-photon Rabi oscillations, see Eq. (8), which agrees with appreciable occupations of states |10 and |01 in Fig. 2(c). Both Figs. 2(a) and 2(b) demonstrate agreement between analytic and numerical calculations at λ 0.25. On the one hand, this agreement is not surprising since Eq. (10) is based on the perturbation theory and RWA in the computational subspace, which are both valid at small λ. On the other hand, Eq. (10) was derived for two-level systems while Fig. 2 presents results for the full fluxonium Hamiltonian (2). To elaborate on the effect of higher levels, we introduce an analog of the dimensionless drive amplitude λ for the |1 A − |2 A transition of qubit A with similar reasoning applied to other transitions: where Ω 1−2 A,0 = 2f n A 12 is the corresponding single-qubit Rabi frequency. For the parameters of Table I, λ 1−2 ≈ 0.16λ, implying that the single-photon transition |1 A − |2 A is suppressed for all the values of λ discussed in Fig. 2(a). For generic fluxonium parameters, we do not expect λ 1−2 /λ to become large. Therefore, even with higher levels taken into account, λ 1 ensures that the |00 − |11 transition is still dominated by coherent twophoton processes via intermediate virtual states with suppressed sequential one-photon transitions via real states. However, a higher-energy noncomputational level |kl (k > 1 or l > 1) formally generates an additional virtual state, producing an extra term in Eq. (10) as long as N 00,kl N kl,11 = 0, where N kl,k l is defined in Eq. (5b). Such terms are exactly zero at J C = 0 and acquire finite values at J C = 0. For the parameters of Table I, they are negligible in comparison to existing contributions because their denominators, which are determined by ω 00−kl −ω, are much larger than ∆ AB . This is not the case in general because possibly large values of n α 12 or n α 03 may make N 00,kl N kl,11 to be sufficiently large and the corresponding additional contribution to Eq. (10) nonnegligible in comparison to the sum of existing terms, which interfere destructively. Equation (10) works for our choice of parameters when the higher-energy states do not contribute significantly, but this approximation may be less accurate for other choices of fluxonium parameters, e.g., for qubits with larger n α 12 /n α 01 or n α 03 /n α 01 .
At λ 0.25, perturbation theory breaks down and the Rabi frequency increases slower than λ 2 . We emphasize that Ω still increases monotonically even far from the perturbative regime, when single-photon transitions can be strongly excited. It can become as large as a few megahertz even for a relatively small J C /h of 100 MHz and can surpass 10 MHz for stronger interaction strengths. We note that Ω is close to being linear in J C , which qualitatively agrees with Eq. (13). Accounting for hybridization between computational and higher-energy states is necessary for a quantitative agreement.
The ability to induce two-photon Rabi oscillations can be used to realize an entangling gate involving the mixing of states |00 and |11 similar to bSWAP of Ref. [20]. A proper pulse shape is required to minimize single-photon transitions [e.g., |00 − |10 in Fig. 2(c)] at the end of the pulse. In addition, driving the |00 − |11 transition with a different ω d creates off-resonant Rabi oscillations that exchange only a fraction of the population between states |00 and |11 . In Fig. 2(d), ω d is chosen to create oscillations of 50% of the state populations. Here, the minima of P 00→01 (t) and P 00→10 (t) occur at times where the state of the system is in an equal superposition of |00 and |11 . This feature is not a generic property of our gate, but it results from our choice of values for the parameters λ and J C . The resulting period of oscillations is 103 ns, which is approximately √ 2 shorter than 147 ns in Fig. 2(c). This behavior is reminiscent of that of a driven two-level system; the period of detuned Rabi oscillations with 50% contrast is exactly √ 2 times shorter than the period of resonant Rabi oscillations at a fixed drive amplitude. It supports our understanding of the high-power drive of |00 − |11 atω in a coupled-fluxonium system in terms of two-photon Rabi oscillations. Note, however, that in a true two-level system, 50% contrast of Rabi oscillations requires a frequency detuning equal to the resonance Rabi frequency. While this Rabi frequency is 6.8 MHz for Fig. 2(c), the difference between the values of ω d in Figs. 2(c) and 2(d) is 4.5 MHz. Therefore, our reasoning in terms of Rabi oscillations in a driven two-level system is correct only qualitatively, while an accurate description of the dynamical behavior requires accounting for other levels.

A. Theoretical concepts
We parameterize the family of gates spanned by the coherent mixing of |00 and |11 states and controlledphase operations by two angles aŝ In the absence of leakage outside of the computational subspace and provided that single-photon processes are negligible, any two-photon process described in the previous section can be reduced to the form (19) by means of single-qubit Z rotations applied from both sides of the operator, which can be implemented as virtual Z rotations in experiments [47].
In addition to the SWAP-like interaction -or XX − Y Y interaction -described by θ in Eq. (19), the microwave drive creates a ZZ interaction between computational states. This coupling leads to a finite ζ in Eq. (19), which cannot be changed to zero by local (single-qubit) operations. The ZZ term has two distinct contributions: the static ZZ coupling, which is caused by the repulsion between computational and noncomputational levels due to interaction (3), and the ZZ coupling induced by the microwave drive used to perform the gate operation. While the effect of the static ZZ coupling is relatively weak (the phase accumulation rate is slightly below 1 MHz for the parameters of Table I with J C /h = 200 MHz) and leads to a small contribution to |ζ| π for short gate durations, the microwave-induced contribution to ζ can be large, as demonstrated in Eq. (15). Thus, one has to include the angle ζ in the definition of the target gate, Eq. (19), to take into account this additional term caused by the drive. As we demonstrate below for three choices of the mixing angle θ, the value of ζ affects the equivalence class of the gate and its entangling properties.
In general, two gates U 00−11 (π, ζ) having different values of ζ in the interval between 0 and π are not locally equivalent. Each class of locally equivalent gates is characterized by special invariants G 1 = G 1 + iG 1 and G 2 , and two gates are locally equivalent if and only if they have the same invariants [48,49]. We calculate them in Appendix A and find the following values for U 00−11 (π, ζ): Another important property of a two-qubit gate is the entangling power P(θ, ζ) [32,33]. It ranges between P = 0 forÛ 00−11 (π, π) (equivalent to a SWAP gate), and P = 2/9 forÛ 00−11 (π, 0) (equivalent to an iSWAP gate). For arbitrary ζ, the entangling power is given by see Appendix B.
The independence of the entangling power on values of ζ for θ = π/2 makes this mixing angle an attractive choice in situations when the induced ZZ coupling is hard to control. Another benefit of implementing a gate with θ = π/2 vs a gate with θ = π is that the former gate can be realized with any off-resonant two-photon Rabi oscillations as long as their contrast V = max t P 00→11 (t) is at least 0.5, while θ = π requires precise swapping of populations via half a period of a resonant Rabi rotation. For example, Fig. 2(d) demonstrates that it is possible to achieve θ = π/2 by choosing half a period of the off-resonant Rabi oscillations with V = 0.5. The drive detuning from the two-photon resonance, ω d = (E 11 − E 00 )/2 , changes the contrast V and the period of Rabi oscillations, which, in turn, affects the gate duration and ζ inÛ 00−11 (π/2, ζ), see, e.g., Eq. (15) for ζ at the resonant drive frequency, when V = 1. Thus, the detuning acts as an additional control, which can be used either in the optimization procedure to improve gate performance when a specific ζ is not needed or in producing a gate with specific ζ.
We also note thatÛ 00−11 (π/2, ζ) is sufficient to realize a gate given by the unitary (19) with any mixing angle θ by combining two gatesÛ 00−11 (π/2, ζ) with single-qubit Z rotations. Some of those Z rotations can be substituted by a change of the microwave-drive phase γ d in the drive term (4) for one of the two-qubit gates. More details are given in Appendix C. This decomposition is similar to the decomposition in Ref. [38] for XY gates, which are excitation-preserving swapping gates activating coherent rotations in the {|01 , |10 } subspace. We note thatÛ 00−11 (π/2, ζ) is not a Clifford gate for any ζ, which requires its characterization via the cross-entropy benchmarking [50] rather than via randomized benchmarking [51].

Mixing angle θ = 0
We finally consider the case when θ = 0. This gate occurs after a full two-photon Rabi oscillation with resonant, θ = π, or off-resonant drive, e.g., θ = π/2. The gate is equivalent to the controlled-phase gate with the invariants determined by ζ and coinciding with invariants of the controlled-phase gate: The entangling power is also ζ dependent:

B. Simulated coherent gate fidelity
Here we demonstrate via numerical simulations that a fast and high-fidelity gate that mixes states |00 and |11 by means of a monochromatic microwave drive is possible. We focus on the gate operation for the mixing angle θ = π/2 and calculate the gate fidelity to an ideal unitaryÛ (π/2, ζ) with a suitable choice of ζ angle. We start from a detailed analysis of coherent gate dynamics and discuss incoherent error in Sec. IV C. As in Sec. III C, we perform simulations for qubit parameters shown in Table I and for fixed and equal drive couplings in Eq. (4), For a given gate duration t gate , we use the Gaussian pulse shape with the rising edge given by Eq. (17), where t rise = t gate /2 and σ = t rise /2. The gate is optimized over the drive amplitude and frequency at fixed tgate = 50 ns for the target operatorÛ00−11(π/2, ζ). The gate fidelity is F ≈ 0.99905 and ζ ≈ 1.02π.
single-qubit Z rotations before and after the gate, we adjust phases of relevant matrix elements ofÛ sim to compare it withÛ 00−11 (π/2, ζ); more details are given in Appendix D. For the target unitaryÛ 00−11 (π/2, ζ), we choose ζ = ϕ 01,01 + ϕ 10,10 − ϕ 00,00 − ϕ 11,11 , where ϕ kl,kl = arg kl|Û sim |kl is the diagonal-matrix-element phase ofÛ sim , so a new value of ζ is chosen each time a newÛ sim is computed. Denoting the simulated gate operator after the application of Z rotations asÛ , we use the standard expression for the two-qubit-gate fidelity [52]: Using this metric, we optimize the coherent gate error over the drive frequency and amplitude. We analyze the dependence of F on the total gate duration t gate and the interaction strength J C .
As an example, we first illustrate the gate operation in time domain in Fig. 3 for such an optimized gate with t gate = 50 ns and J C /h = 200 MHz. We show transition probabilities P kl→k l vs time for all 16 pairs of initial and final states formed from {|00 , |01 , |10 , |11 }. The simulated gate is found to be in the equivalence class determined by θ = π/2 and ζ ≈ 1.02π with the gate fidelity being F ≈ 0.99905. The gate duration of t gate = 50 ns is the shortest time for which 1 − F < 0.001 for the interaction strength used in this simulation. optimized gate of Fig. 3 exhibits multiple fluctuations and large transient excitations of one-photon processes, the gate performance is actually very robust to calibration errors. We illustrate this statement in Figs. 4(a) and 4(b), where we study gate properties around the optimal point by changing the drive frequency and amplitude. We note that 1 − F is below 0.001 in the frequency interval greater than 1 MHz, see the solid line near the vertical arrow in Fig. 4(a). Longer gate durations can be chosen for which this interval is even wider.
We analyze various coherent contributions to the gate error in Fig. 4 with more details given in Appendix D. For the ideal gate operationÛ 00−11 (π/2, ζ), the transition probabilities P kl→k l are either zero, 1/2, or 1, while the actual gate operation contains errors in those probabilities. When errors are small, 1−F is well approximated by the sum of two distinct contributions that are linear in P kl→k l . The first contribution, E comp , is determined by those P kl→k l for which |k l is in the computational subspace; see Eq. (D5). The second contribution, P leak , is the leakage error given by the average probability to end up outside of the computational subspace, see Eq. (D7). We show these two contributions by dashed and dash-dot lines in Fig. 4(a) and 4(b). We find that at the optimal point, the coherent gate error 1 − F is determined by the error in the computational subspace E comp . Leakage errors are below 0.0001 at the optimal point. A more detailed analysis (not shown here) indicates that the remaining leakage errors are mostly coming from excitations of the second excited states of fluxonium circuits via, e.g., transitions |10 − |20 and |11 − |12 .
While E comp and P leak are sufficient to explain the behavior of 1 − F near its minimum, other error mechanisms become dominant away from the optimal point. For instance, we find a large contribution from the mixing error E θ . This error is determined by the differences |P 00→00 −P 00→11 | and |P 11→00 −P 11→11 | and is quadratic in them. Essentially, it can be thought of as the error in the mixing angle θ inÛ 00−11 (θ, ζ). We show this error by the dotted line in Figs. 4(c) and 4(d), which explains well the behavior of 1 − F far away from the optimal point.
In addition to the gate fidelity, we calculate the concurrence C 00 (see Ref. [53]) of a state vector starting in |00 after the application of the gate, see Appendix D. It is shown by the dashed lines in Figs. 4(c) and 4(d). When C 00 = 1, the state is maximally entangled and the gate is thus a perfect entangler. While the mixing error has two almost symmetric minima, which are sharp and deep, the two minima of the concurrence error are asymmetric. This is explained by the dependence of the concurrence not only on the mixing of states |00 and |11 , but also on both amplitudes and phases of states |01 and |10 , see Eq. (D10). The corresponding contributions coming from |01 and |10 have opposite signs in the left and right minima.
In Fig. 5, we show coherent gate error and its budget as a function of gate duration t gate and the interaction strength J C . For J C /h = 200 MHz, we observe that the coherent error can easily go below 10 −4 for a gate duration shorter than 100 ns. In the bottom panels, we show parameter ζ that determines the equivalence class of the gate. Its tendency to decrease with increasing t gate or J C is explained by the contribution to ζ coming from the static ZZ coupling. The effect of static ZZ grows with the gate duration and interaction strength. Finally, we note that combining twoÛ 00−11 (π/2, ζ) gates, we can obtain a controlled-phase gate with phase 2ζ.

C. Incoherent error
In this section, we discuss how qubit decoherence affects the gate error. We consider relaxation and dephasing of only |0 α − |1 α and |1 α − |2 α transitions since the qubit excitation probability above its second excited state is very small. For example, for the parameters of states during 50 ns of gate operation. For states |2 α , these numbers are about an order of magnitude larger (5%, 1.2%, and 0.6 ns) and are thus not large either, but may still result in an important contribution to the gate error since the coherence time of |1 α − |2 α is often significantly shorter than the coherence time of the computational subspace. A small population of state |2 α is consistent with a small λ 1−2 , defined in Eq. (18).
We simulate the gate operation in the presence of decoherence for optimal pulse parameters found in simulation of unitary dynamics in Sec. IV B. The evolution of a density matrix ρ is described by the Lindblad master equation Here T k−l,A This approach to describe relevant noise sources is sufficient for our purpose of providing a crude estimate.
Using these collapse operators, we perform numerical quantum process tomography. We simulate the superoperator describing the evolution of density matrices corresponding to master equation (26), project the operator into the computational subspace, and use it to find the 16 × 16 χ matrix χ real describing the quantum process. We find the ideal χ matrix χ ideal usingÛ 00−11 (π/2, ζ) in Eq. (25) modified by single-qubit Z rotations used to obtainÛ fromÛ sim in unitary simulations; see the text above Eq. (25). We then use which establishes a relation between the gate and process fidelities, where the latter is given by Tr(χ † real χ ideal ) [54]. Using this approach, we study how gate error depends on relaxation and dephasing times. For each transition |k α − |l α , we assume that its relaxation (T k−l,α 1 ) and coherence (T k−l,α 2 ) times are the same, so its pure dephasing time is T k−l,α ϕ = 2T k−l,α 1 . We also assume that these times are the same for both qubits, but different for the two transitions, so we use two lifetime parameters: T 0−1 for relaxation and dephasing of the computational transitions of both qubits and T 1−2 for the |1 α − |2 α transitions. In the top panels of Fig. 6, we show the gate error calculated for Fig. 5(a) parameters at t gate = 50 ns, which are marked by vertical arrows in Fig. 5 and were also discussed in Figs. 3 and 4. In the bottom panels of Fig. 6, we consider Fig. 5(a) parameters at t gate = 93 ns, which is a local minimum of the coherent error. We study separately the effects of decoherence of |0 α − |1 α and |1 α − |2 α transitions. Thus, the left panels of Fig. 6 show 1 − F vs T 0−1 assuming that T 1−2 = ∞ and the right panels discuss T 0−1 = ∞ and finite T 1−2 . Horizontal dashed lines show the coherent error, which was calculated in the previous section, so the difference between symbols and lines is the incoherent contribution coming from either |0 α − |1 α or |1 α − |2 α transitions. The total incoherent error is approximately given by the sum of incoherent errors in the left and right panels. Figure 6 demonstrates that the contribution to the incoherent error coming from |1 α − |2 α transitions is much less important than that coming from the relaxation and dephasing in the computational subspace. This is consistent with small average occupations of the second excited states of fluxoniums during the gate operation. We observe that even a very short T 1−2 ∼ 1 µs results in only about 10 −3 contribution to the gate error, while T 1−2 ∼ 20 µs is sufficient to contribute less than 0.5 × 10 −4 at t gate = 93 ns. Several dozens of microseconds for a coherence time of a transition with frequency in the gigahertz scale is common in modern state-of-the-art superconducting qubits [6] with the best lifetimes exceeding 100 µs [55,56]. Therefore, we do not expect the coherence time of the |1 α − |2 α transition to be a limiting factor for the proposed gate. The contribution from decoherence of the computational transitions is more important. We find that T 0−1 ∼ 100 µs is generally sufficient to bring the gate error below 10 −3 , while the 10 −4 threshold requires T 0−1 > 1 ms. We note that the best fluxonium devices have recently demonstrated lifetimes of 1 ms, although at lower transition frequencies than those discussed here [26]. Nevertheless, we do not see any fundamental obstacles in achieving a millisecond lifetime of the fluxonium with 1 GHz transition, which paves the way towards 10 −4 gate errors.

V. DISCUSSION AND CONCLUSIONS
We demonstrated that fast high-fidelity microwaveactivated two-qubit gates are possible in fluxonium circuits when the system state remains entirely in the lowenergy computational subspace. Despite a relatively weak effect of capacitive interaction between fluxoniums on the computational subspace, the gate time can still be as short as 50 ns due to the strong anharmonicity of the fluxonium spectrum. The anharmonicity typically limits the intensity of the microwave drive. We demonstrated that the microwave amplitude could be large for the proposed two-photon gate without causing noticeable leakage of the state outside of the computational subspace during the pulse, minimizing the effect of decoherence on the gate fidelity. The required amplitude is about 10-20 times larger than its value in schemes utilizing noncomputational levels [29,30] and in single-qubit operations. This strength is on par with the cross-resonance gate, which is activated by microwave fields with resonance Rabi frequencies up to hundreds of megahertz [18] and which has techniques to mitigate cross-talk and spectator errors in transmon processors [22,57]. Strong anharmonicity of the fluxonium and an extra freedom in choosing the qubit frequency will likely make mitigation of these errors even more successful in a fluxonium-based processor.
At a weak drive power, the rate of two-photon transitions is quadratic in the drive amplitude. In this case, the gate would be prohibitively long if the drive were chosen so that single-photon transitions between subspaces {|00 , |11 } and {|01 , |10 } were strongly suppressed. We demonstrated that unintended single-photon transitions between those subspaces could be reduced even for a strong drive by fine tuning the pulse amplitude and frequency together with the gate duration. As a result, the microwave pulse only mixes states |00 and |11 . We focused on the half-mixing angle θ = π/2 for which the entangling power is independent of the phase shift due to the ZZ interaction, which guarantees that the gate is entangling without any control of the Stark-induced phase accumulation. We also note that the half-mixing gate is shorter and is often more robust to pulse imperfections.
In conclusion, we considered a two-qubit gate that is well suited for existing fluxonium devices and is ready to be implemented. The proposed scheme is very generic and can also be realized as a two-color scheme with two microwave drives at two different frequencies ω d1 and ω d2 satisfying ω d1 + ω d2 = ω 00−11 , which provides additional controls to reduce errors. The gate works for fluxoniums parked at the sweet spots of their maximal coherence and does not require any additional hardware beyond microwave control lines necessary to activate single-qubit gates.

ACKNOWLEDGMENTS
We would like to thank Mark Dykman, Ivan Pechenezhskiy, Haonan Xiong, and Long Nguyen for fruitful discussions. We acknowledge the support from NSF PFC at JQI and ARO-LPS HiPS program (grant No. W911NF-18-1-0146). V.E.M. and M.G.V acknowledge the Faculty Research Award from Google and fruitful conversations with the members of the Google Quantum AI team. We used the QuTiP software package [58,59] and performed computations using resources and assistance of the UW-Madison Center For High Throughput Computing (CHTC) in the Department of Computer Sciences. The CHTC is supported by UW-Madison, the Advanced Computing Initiative, the Wisconsin Alumni Research Foundation, the Wisconsin Institutes for Discovery, and the National Science Foundation.

Appendix A: Local invariants
Here we calculate the invariants [48,49] for the gates described by Eq. (19) and demonstrate that the Stark shift ζ makes the gates with different ζ nonequivalent to each other.

Appendix B: Entangling power
Here we provide expressions for calculations of the entangling power P of a two-qubit operator (19) [32]. We use the algebraic technique of Ref. [33] that defines the entangling power as where E(U ) is the operator entanglement (linear entropy) of U : the matrix U R is obtained from U by realignment: (U R ) ij,kl = U ik,jl , and S 12 is the swapping operator We take U (θ, ζ) in the form of Eq. (19) and find that P(θ, ζ) = 1 36 (5 − 4 cos ζ cos θ − cos 2θ) .
We note that the entangling power is independent of the Stark shift ζ only for θ = π/2.