Fast and differentiable simulation of driven quantum systems

The controls enacting logical operations on quantum systems are described by time-dependent Hamiltonians that often include rapid oscillations. In order to accurately capture the resulting time dynamics in numerical simulations, a very small integration time step is required, which can severely impact the simulation run time. Here, we introduce a semianalytic method based on the Dyson expansion that allows us to time-evolve driven quantum systems much faster than standard numerical integrators. This solver, which we name Dysolve, efficiently captures the effect of the highly oscillatory terms in the system Hamiltonian, significantly reducing the simulation's run time as well as its sensitivity to the time-step size. Furthermore, this solver provides the exact derivative of the time-evolution operator with respect to the drive amplitudes. This key feature allows for optimal control in the limit of strong drives and goes beyond common pulse-optimization approaches that rely on rotating-wave approximations. As an illustration of our method, we show results of the optimization of a two-qubit gate using transmon qubits in the circuit QED architecture.

Figure 1

Tree diagram showing the different branches of the time-ordered integral, with two example branches highlighted.

Figure 2

An example set of pulse amplitudes. The black bars indicate the chosen drive amplitudes $u_j$, where $\mathrm{\Delta }t=1$ ns. The red bars indicate the subpixels, which provide an approximate interpolation of the true pulse delivered to the system (blue), with a bandwidth ${\omega }_0/2\pi =851$ MHz.

Figure 3

Dysolve benchmarks. Frobenius norm distance from the propagator ${\widehat{U}}_{\mathrm{ref}}(0,T)$ for the Dysolve algorithm for a $T=500$ ns evolution with (a) one input drive, (b) two input drives, and (c) three input drives. Contraction time of the Dysolve algorithm as a function of the subpixel number $N_s$ for (d) one, (e) two, and (c) three input drives. The slope of the data is precisely 1, meaning that the computational time scales linearly with the number of subpixels. Preparation stage computation times are (10 ms, 18 ms, 24 ms) for (Dysolve-2, 3, 4) and one input drive, (22 ms, 68 ms, 136 ms) for two input drives, and (45 ms, 200 ms, 442 ms) for three input drives. In panels (a) and (d), we compare the error and computation time, respectively, of Dysolve with qutip bdf and adams solvers using high-precision settings (see Appendix pp4 for more details).

Figure 4

Superconducting circuit for cross-resonance gates between transmon qubits. The leftmost transmon, of frequency ${\omega }_{\mathrm{c}}/2\pi$, plays the role of the control qubit and is strongly driven by the voltage source $V_{\mathrm{c}}$ at the frequency ${\omega }_{\mathrm{t}}/2\pi$ of the target qubit (rightmost transmon). The target qubit is also driven by $V_{\mathrm{t}}$ using a relatively weak tone that serves to give additional control. ${\widehat{\phi }}_{\mathrm{c}}$ and ${\widehat{\phi }}_{\mathrm{t}}$ correspond to the phase degrees of freedom associated with the control and target qubits, respectively.

Figure 5

Infidelity of the cross-resonance gate as a function of the chosen gate time when the gate is operated at the qubit-qubit detuning $\mathrm{\Delta }/2\pi =({\omega }_c-{\omega }_t)/2\pi =210$ MHz. The different curves demonstrate the result of optimization of the pulse shape under RWAs at different threshold frequencies ${\omega }_{\mathrm{thres}}$, followed by a subsequent calculation of the true fidelity in the RWA-free case.

Figure 6

Infidelity of the cross-resonance gate as a function of the detuning of the qubits when optimized under GRAPE for 5000 iterations. The control qubit frequency is fixed to ${\omega }_c/2\pi =5.10$ GHz.

Figure 7

Linear interpolation of the pulse sequences (solid yellow lines) with two subpixels per pixel. The black bars indicate the original Gaussian filtering, with the blue dashed line the “true” filtered pulse.

Figure 8

Comparison of the Dysolve algorithm and qutip's propagator over an evolution of 20 Rabi oscillations, with the same parameters as in Fig. 3 of the main text. (a) Frobenius norm distance between the Dysolve algorithm with different subpixel numbers and qutip's propagator with adams (circles) and bdf (x's) methods. The squares are the norm distance between Dysolve-n and Dysolve-4 with ${10}^4$ subpixels. The dashed line indicates the error between propagator with the adams and bdf methods. (b) Computational time of Dysolve and propagator for results of the left panel. The dashed lines refer to the computational time for propagator with adams and bdf with the high-precision settings described in Appendix pp4.