Fast elementary gates for universal quantum computation with Kerr parametric oscillator qubits

Kerr parametric oscillators (KPOs) can stabilize the superpositions of coherent states, which can be utilized as qubits, and are promising candidates for realizing hardware-efficient quantum computers. Although elementary gates for universal quantum computation with KPO qubits have been proposed, these gates are usually based on adiabatic operations and thus need long gate times, which result in errors caused by photon loss in KPOs realized by, e.g., superconducting circuits. In this work, we accelerate the elementary gates by experimentally feasible control methods, which are based on numerical optimization of pulse shapes for shortcuts to adiabaticity. By numerical simulations, we show that the proposed methods can achieve speedups compared to adiabatic ones by up to six times with high gate fidelities of 99.9%. These methods are thus expected to be useful for quantum computers with KPOs.

Applications of KPOs to fault-tolerant quantum com- * taro.kanao@toshiba.co.jp puting [56] have also been studied [57].Quantum gates preserving the bias of errors mentioned above have been proposed [58], which can be utilized for hardware-efficient quantum error correction [59].Analytically engineered control methods for shortening the gate times of the biaspreserving gates have recently been proposed [13].Furthermore, for noisy intermediate-scale quantum (NISQ) applications [60], variational quantum algorithms [61,62] for KPOs have been proposed, such as quantum supervised machine learning [63] and a quantum approximate optimization algorithm [64].
For KPO qubits, elementary gates for universal quantum computation have been proposed [8,9], which are based on adiabatic evolution and consist of Z, X, and ZZ rotations denoted by R z , R x , and R zz , respectively.Experimentally, a study [19] has demonstrated adiabatic R z and nonadiabatic R x , and another study [71] has adiabatically performed both R z and R x .Theoretically, other kinds of gate implementations have been proposed [15,[73][74][75][76][77].
Shorter gate times are desirable, because they can reduce errors caused by photon loss in KPOs and also enable faster computation.However, the previous adiabatic elementary gates [8,9] need long gate times and otherwise diabatic transitions out of a qubit space cause leakage errors.To reduce leakage errors, in this work, we focus on control methods called shortcuts to adiabaticity (STA) [78].For KPOs, STA have been proposed for cat-state generation [9,40] and R zz with a phase rotation of a parametric drive [15].Also, a variant of the derivative removal by adiabatic gate (DRAG) technique, which is related to STA, has been proposed for the bias-preserving gates [13].
To accelerate the elementary gates for universal quantum computation with KPO qubits, our approach is based on an STA called counterdiabatic terms (or counter terms for short) [79,80], but does not use the exact counter terms, which are often experimentally infeasible.Instead, we first approximate the counter terms by experimentally feasible terms [81] and then numerically optimize the pulse shapes for the gate operations.As a result, we successfully shorten gate times, keeping high gate fidelities.By this approach, the gate operations become faster by 2.6 times for R z , 6.0 times for R zz , and 2.6 times or higher for R x .Interestingly, the states of KPOs during the optimized gate operations for the shortest gate times are not necessarily instantaneous eigenstates, which indicates that the numerical optimization explores gate operations beyond the STA.We also numerically show that the optimized gate operations are robust against systematic errors in the amplitudes of gate pulses, and the shortened gate times can suppress errors caused by single-photon loss.We expect that these optimized elementary gates for KPO qubits will be useful for NISQ applications in a near term and fault-tolerant quantum computation in a long term.

II. APPROXIMATE STA A. Elementary gates for KPO qubits
We first introduce the model of the KPO and elementary gates for universal quantum computation with the KPO qubits [7,9].In a rotating frame and within the rotating-wave approximation, the Hamiltonian for a KPO is given by [6] where a, K, and p are the annihilation operator, the Kerr coefficient, and the amplitude of a parametric drive, respectively.In this study, the reduced Plank constant ℏ is set to 1.The two degenerate eigenstates of the Hamiltonian corresponding to effective ground states of the KPO [73] are written as where |±α⟩ are coherent states with an amplitude α = p/K.In this work, we use the following computational basis [58,73], which are exactly orthogonal.Equations ( 3) and ( 4) are approximately equal to |±α⟩, respectively, for p/K = 4 used in this study.
For the KPO qubits, elementary gates for universal computation can consist of Z, X, and ZZ rotations, which are expressed respectively as [56] R z (ϕ) = e −iϕ/2 0 0 e iϕ/2 , (5) where ϕ, θ, and Θ are respective rotation angles.For universal computation, arbitrary ϕ, θ = π/2, and Θ = π/2 are enough [8,9,56].For KPOs, these elementary gates can be implemented based on adiabatic control with a single-photon drive, a detuning, and a linear coupling, respectively.The Hamiltonians corresponding to the single-qubit gates are where g 0 (t) is the amplitude of a gate pulse.A linear coupling necessary for R zz can be realized with beam-splitter coupling [8,9] or two-mode squeezing [13], described by where a i and H KPOi are the annihilation operator and the Hamiltonian in Eq. (1) for the ith KPO.

B. Approximate counter terms for STA
An ideal counter term H 1 (t) for STA exactly reproduces adiabatic evolution with H 0 (t) by finite-time evolution with H 0 (t)+H 1 (t) [78][79][80], but is often experimentally infeasible.Under certain assumptions, we approxi-mate H 1 (t) by where the dots denote the time derivative (See Appendix A for the details of the assumptions and derivations).Importantly, these H 1 (t) are experimentally feasible as follows.
• H 1 (t) for R z in Eq. ( 14) can be implemented with a single-photon drive with its phase shifted by π/2 from the original single-photon drive in Eq. ( 9).
• H 1 (t) for R x in Eq. ( 15) can be realized by a twophoton drive with its phase shifted by π/2 from the original parametric drive.
• H 1 (t) for R zz in Eq. ( 16) is another two-mode squeezing than the original one in Eq. ( 13), which can be realized in a previously proposed superconducting circuit for R zz [74].
Note that the counter term in Eq. ( 16) can be derived from both R zz by the beam-splitter coupling in Eq. ( 12) and R zz by the two-mode squeezing in Eq. ( 13).However, we numerically find that the two-mode squeezing in Eq. ( 13) gives better results with the counter term in Eq. ( 16), which can be understood from the matrix elements of A 0 and H 1 (t) (see Appendix B 1 for details).We thus use the two-mode squeezing in Eq. ( 13) in the following.
Here we also comment on another candidate of a counter term for R zz , We numerically found that this term does not work as a counter term (a similar result has been mentioned in Ref. [13]).Equation (17) does not cancel unwanted transitions out of the qubit space, because H 1 (t) in Eq. ( 17) and A 0 in Eqs. ( 12) and ( 13) have different permutation symmetry, namely, symmetry with respect to the interchange of KPO1 and 2 (see Appendix B 1 for details).

C. Numerical optimization
To go beyond the analytic approximate H 1 (t) in Eqs. ( 14)-( 16), our proposed approach uses arbitrary waveforms for the amplitudes of the counter pulses g 1 (t) as and numerically optimizes g 1 (t) as well as g 0 (t) in Eqs. ( 8) and (11).Total Hamiltonians are then given by, for the single-and two-qubit gates, respectively, where A 0 are given in Eqs. ( 9), (10), and (13).Here, the two-mode squeezing Hamiltonian in Eq. ( 13) is used as mentioned above.We expect that this approach, which numerically optimizes pulse shapes for STA, will be useful for other qubit systems.
To optimize g 0 (t) and g 1 (t) numerically, we express the waveforms of the pulse amplitudes by [82] where T is a gate time and N f determines the number of frequency components.We choose the symmetric g 0 (t) and antisymmetric g 1 (t) with respect to time reversal t → T −t, because an exact counter term is antisymmetric when the other term is symmetric (see Appendix B 2).In g 0 (t), we include the sine terms to allow for nonzero ġ0 (t) at t = 0, T [8].Since the highest frequencies in g j (t) are limited to N f /T and g j (t) are zero at initial and final times, these waveforms are expected to be experimentally feasible.
We numerically optimize g j,n in Eqs ( 24) and ( 25) to maximize an average gate fidelity F [83,84] given in Eq. (C1) in Appendix C, using the quasi-Newton method with the BFGS formula [85].We set the initial values of g j,n for the optimization to the ones corresponding to analytic waveforms for adiabatic elementally gates without and with the counter terms in Eqs. ( 14)-( 16) (see Appendix D).

III. NUMERICAL SIMULATIONS
In the present simulations, we regard the Kerr coefficient K as the unit of the frequency and set the amplitude of the parametric drive to p = 4K.We express states and operators in the photon-number basis with the largest photon number of 40, which is large enough.We simulate the time evolution of states by numerically solving the Schrödinger equation unless otherwise stated.We use the fourth-order Runge-Kutta method with the step size of 10 −4 /K.We compare the following four cases depending on the waveforms and the counter terms: 1. Analytic waveforms without the counter terms.
In this work N f in Eqs. ( 24) and ( 25) is set to 10.The analytic waveforms and initial values of g j,n for the numerical optimization are given in Appendix D.

A. Simulation results for Rz
Figure 1(a) shows the average infidelities 1 − F for R z (π) as functions of the dimensionless gate time KT .The infidelities decrease with increasing gate time, indicating the adiabaticity of the gate, where the errors are mainly due to the leakage of population to the states outside the qubit space.We define a minimum gate time T min by minimal T satisfying 1 − F < 10 −3 , and compare T min for the above four cases.With analytic waveforms, KT min are 1.3 and 1.2, respectively without and with the counter term.By the numerical optimization, KT min are shortened to 0.9 and 0.5, respectively.Thus the numerically optimized R z with the counter term is 2.6 times faster than the original analytic R z without the counter term.These results show that the counter term is effective and the improvement is enhanced by the numerical optimization.
We examine the optimized gate operation with the counter term at KT min = 0.5.The optimized waveforms of g j (t) are shown in Fig. 1(b).Figure 1(c) shows the resulting time evolutions of the mean photon number and population in the qubit space with the initial state |C + ⟩, where P is a projector onto the computational basis states, It is notable that despite the large amplitudes of g j (t), the mean photon number and the population in the qubit space are almost unchanged.
To see the state in more detail, we use the Wigner function W (x, y), which is a quasiprobability distribution for x = a + a † /2, y = a − a † /(2i) [86] and is calculated by the technique in Ref. [7]. Figure 2 shows W (x, y) during the gate operation with the optimized g j (t) in Fig. 1 can preserve the coherent states when the effective potential of the KPO is well approximated by the double well [12,57].Interestingly, we numerically found that the cat states in Fig. 2 are not instantaneous eigenstates of H 0 (t), which indicates that our proposed approach is beyond STA.We next show that the optimized g j (t) with the counter term for R z (π) can be used for R z (ϕ) with arbitrary ϕ by introducing only one time-independent scaling parameter λ.The pulse amplitudes are set to λg 0 (t) and λg 1 (t).Resulting ϕ is determined by maximizing F . Figure 3(a) shows 1 − F as a function of ϕ at KT min = 0.5, which demonstrates that this method gives high-fidelity R z (ϕ) for arbitrary ϕ in 0 ≤ ϕ ≤ π.An exact counter term suggests that this continuous gate by the one parameter λ is possible because the changes in the states are small during the gate operation as shown in Fig. 2 (see Appendix B 3 for details).On the other hand, this continuous gate does not hold for R x (θ) as also mentioned later in Sec.III C, where the states largely change during the gate operation.
To examine the optimality and robustness of the optimized g j (t), we evaluate R z (π) with (1 + δ j ) g j (t) for given relative errors δ j , which can model systematic er- rors in the pulse amplitudes [15].Figure 3(b) shows 1 − F as a function of δ j .First, at δ 0 = δ 1 = 0, the gradient of 1 − F with respect to δ j vanishes, implying that δ 0 = δ 1 = 0 is an optimal point.Second, the ellipse in Fig. 3 0.2.Thus, the numerically optimized gate operation with the counter term provides a speedup by 6.0 times compared with the original analytic waveform without the counter term.The gate-time dependence of the infidelities of R zz in Fig. 4(a) is similar to that of R z in Fig. 1(a), which may be because for one KPO the other acts like a single-photon drive as in R z .Optimized waveforms g j (t) with the counter term at KT min = 0.2 are shown in Fig. 4(b).Figure 5(a) shows that these optimized g j (t) can be used for continuous R zz (Θ) by λg j (t) with the time-independent scaling parameter λ as in the case of R z (ϕ).Also, the optimality and robustness of R zz (π/2) are evaluated with (1 + δ j )g j (t). Figure 5(b) shows that the gradient of 1 − F is zero at δ 0 = δ 1 = 0, indicating its optimality.the original analytic waveform without the counter term.However, we find that for KT ≤ 0.9, the maximum value of |g j (t)|/K can be larger than 50, which might be infeasible because the rotating-wave approximation would be no longer valid [50].Thus, in the following, we examine optimized g j (t) at KT = 1, which is an acceleration by 2.6 times compared with the analytic waveform without the counter term.Then, the pulse amplitudes with |g j (t)|/K < 20 are obtained as shown in As mentioned in Sec.III A, we find that R x (θ) for continuous θ is not obtained by λg j (t) with the optimized g j (t) for θ = π/2.This might be because the states largely change from the two coherent states during the optimized R x (π/2) unlike R z and R zz (see Appendix B 3).

C. Simulation results for Rx
The optimality and robustness of the optimized g j (t) are again evaluated by the gate operation by (1 + δ j ) g j (t). Figure 8 shows the average infidelity as a function of δ j , indicating that the optimality and robustness hold.Compared with R z and R zz , the average infidelity for R x shows larger correlation between δ 0 and δ 1 , that is, the optimized R x is more robust against the relative errors with δ 0 ≃ δ 1 .This result suggests that the counter pulse g 1 (t) plays a more important role in R x than in the others.

D. Effect of single-photon loss
Finally, we evaluate errors in the presence of singlephoton loss, which we choose as a representative of decoherence sources in KPOs.We solve the master equation for a density operator ρ, where is the commutation relation for operators and κ is the loss rate.H(t) are given in Eqs. ( 22) and ( 23), and the optimized g j (t) are used.
Here, we evaluate an average gate fidelity calculated with a finite number of initial states, Floss , which is defined by Eq. (C3) in Appendix C. Figure 9 shows 1 − Floss as functions of the loss rate.For all of R z , R zz , and R x , the errors can be suppressed below 1 − Floss < 10 −3 for the loss rate as large as κ/K ≤ 3×10 −4 , which can be achieved by optimizing the gate times KT as follows.For small κ/K, the errors are dominated by the leakage errors, which depend little on κ/K and decrease with increasing KT .For large κ/K, the errors are mainly due to dephasing caused by the single-photon loss [9], because 1 − Floss with and without the gate operations overlap for the same KT , especially for R z and R zz in Figs.9(

IV. SUMMARY
We have shown that adiabatic elementary gates for universal qunatum computation with KPO qubits can be accelerated by utilizing feasible counter terms derived from STA and numerically optimizing the pulse shapes for them.The optimized gate operations are feasible in experiments with, specifically, superconducting circuits.We thus expect that the proposed methods are useful for quantum computers with KPOs.

Appendix A: Approximate counter terms
We start from an exact counter term H 1 (t) given by [78,80] where |E k (t)⟩ is the kth eigenstate of H 0 (t) with its eigenvalue E k (t), The first term in the right-hand side of Eq. (A1) cancels transitions between different |E k (t)⟩, while the second term corrects phase factors.We ignore the second term because we find that ⟨E k (t)| Ėk (t)⟩ vanishes in the approximation below.We then express H 1 (t) in an explicitly Hermitian form using the time derivative of k |E k (t)⟩⟨E k (t)| = I (I is the identity operator) as To cancel transitions from a qubit space, we restrict the summation to k for computational basis states.From Eq. (A3), we derive approximate counter terms for KPOs.We approximate |E k (t)⟩ corresponding to the qubit states by a variational method using a coherent state |β⟩ as a trial state [28].Its amplitude β is determined by seeking an extremum with For R z , we approximately solve Eq. (A4) by assuming g 0 (t)/ 2Kα 3 ≪ 1 and obtain the following expression for the qubit states, where α = p/K.For real β(t), holds.Equation (A3) can then give Eq.( 14) as where | is a projector onto the qubit space.P (t) is ignored in Eq. (A9), because the difference between Eqs. (A8) and (A9) is proportional to a † [I − P (t)] − [I − P (t)]a, which has matrix elements mainly between states outside the qubit space, and therefore is negligible.

FIG. 1 .
FIG. 1.(a) Average infidelities for Rz(π) as functions of gate time for an analytic waveform of a pulse amplitude without a counter term (Analytic, H1) and with it (Analytic, H1 + H2), and for numerically optimized waveforms without the counter term (Optimized, H1) and with it (Optimized, H1 +H2).The line indicates 1 − F = 10 −3 .(b) Waveforms of the amplitudes of a gate pulse g0(t) and a counter pulse g1(t) for Rz(π), which are numerically optimized for KTmin = 0.5.(c) Mean photon number and population in the qubit space during Rz(π) with gj(t) in (b).The initial state is |C+⟩.

FIG. 4 .
Figure 4(a) shows 1 − F for R zz (π/2) as functions of KT .With analytic waveforms, the minimum gate time satisfying 1 − F < 10 −3 is KT min = 1.2 and 0.6 without and with the counter term, respectively.With the numerical optimization, corresponding KT min are 0.8 and

Figure 6 shows 1 −FIG. 5 .
Figure6shows 1 − F for R x (π/2) as functions of KT .With analytic waveforms, KT min are 2.6 both with and without the counter term.With the numerical optimization, KT min are 1.7 and 0.6 without and with the counter term, respectively, which means that our approach can achieve a 4.3 times faster gate operation than that with
a) and 9(b).The dephasing errors increase with KT [cf.Eq. (C5) in Appendix C].For intermediate κ/K, the leakage and the dephasing compete, and thus 1− Floss are the smallest at larger KT than KT min , achieving 1 − Floss < 10 −3 for κ/K ≤ 3 × 10 −4 as above.In Fig. 9(c), for large κ/K, R x gives smaller 1 − Floss than no gate operation, because the mean photon number is decreased during R x as shown in Fig. 7(b) and the effect of single-photon loss becomes smaller [cf.Eq. (C5) in Appendix C].

1 .
Approximate counter term for Rz
Figure10shows the Wigner function during the optimized R x in Fig.7(a).Figure10(b) shows that for |ψ init ⟩ = |C + ⟩ the intermediate state looks like a vacuum state, which agrees with the small mean photon number in Fig.7(b).The vacuum state may be realized because for large |g 0 (t)|/K, H 0 (t) can be approximated by g 0 (t)a † a and the vacuum state becomes an eigenstate.On the other hand, Fig.10(e)shows that for |ψ init ⟩ = |C − ⟩ the intermediate state resembles |C − ⟩, which is consistent with the large population in the qubit space for Kt = 0.5 in Fig.7(c).
Figure10shows the Wigner function during the optimized R x in Fig.7(a).Figure10(b) shows that for |ψ init ⟩ = |C + ⟩ the intermediate state looks like a vacuum state, which agrees with the small mean photon number in Fig.7(b).The vacuum state may be realized because for large |g 0 (t)|/K, H 0 (t) can be approximated by g 0 (t)a † a and the vacuum state becomes an eigenstate.On the other hand, Fig.10(e)shows that for |ψ init ⟩ = |C − ⟩ the intermediate state resembles |C − ⟩, which is consistent with the large population in the qubit space for Kt = 0.5 in Fig.7(c).