Simulating nonnative cubic interactions on noisy quantum machines

Ideal quantum computers can in principle perform arbitrary unitary operations. However, a lengthy sequence of standard gates is required when the unitary operation results from Hamiltonian evolution that is not native to the hardware. In this case, implementing the operation becomes infeasible when the gate depth exceeds the coherence limit of noisy devices. While this issue may ultimately be resolved by quantum error correction [1, 2, 3, 4, 5, 6], in this letter, we show that efficient compilation plays a key role in realizing quantum computing. When compiling a quantum program using modular gates [7], the end performance of quantum hardware is significantly improved due to the reduction of gate depth. In particular, we show that for cubic interactions, which are essential for nonlinear optics [8], gauge theories [9], as well as plasma and fluid dynamics [10], modular gates enable useful Hamiltonian simulations on present-day noisy hardware, which usually lacks native cubic couplings [11, 12].

To close a major gap in the literature, we consider the following underexamined cubic Hamiltonian ($\hbar=1$) that governs the lowest-order nonlinear interactions in quantum optics and plasma dynamics:

where $g$ is the coupling coefficient and $[A_{j},A_{l}^{\dagger}]=\delta_{jl}$. The resulting Heisenberg equations are special cases of three-wave equations

where $d_{t}=\partial_{t}+\mathbf{v}_{j}\cdot\nabla$ is the convective derivative at the wave group velocity $\mathbf{v}_{j}=\partial\omega_{j}/\partial\mathbf{k}_{j}$. As example applications, Eqs. (2)-(4) are solved to (i) optimize the input seed pulse shape when designing nonlinear optical systems in the pump depletion regime [13, 14, 15]; (ii) determine the turbulence spectrum in the wave-kinetic approach to weak turbulence [16], and (iii) compute the energy delivered to fusion fuel when drive lasers exchange energy in plasmas [17, 18, 19]. More broadly speaking, cubic terms are common for nonlinear problems, for which cubic gates that enacts Eq. (1) serve as basic building blocks. Cubic gates, in conjunction with quartic gates that are usually natively available on hardware, allow one to simulate a majority of problems in physics.

In the usual approach to quantum emulation, terms in the hardware Hamiltonian $H_{0}=\sum_{k=1}^{m}H_{k}$ are turned on or off to generate the unitary evolution due a targeted $H=\sum_{k=1}^{m}a_{k}H_{k}$, where $a_{k}$ are real coefficients. Explicitly, using the Lie-Trotter-Suzuki expansion, $\exp(-iHT)=\lim_{n\rightarrow\infty}[\prod_{k=1}^{m}U_{k}(a_{k}T/n)]^{n}$, where $U_{k}(t)=\exp(-iH_{k}t)$. Difficulty arises when $H$ contains terms that are not natively available in $H_{0}$. Although cubic couplings may be realized natively with special-purpose hardware [20, 21], direct emulations have two major limitations. First, qubits can only be in a superposition of $|0\rangle$ and $|1\rangle$ states, whereas a large amplitude wave, such as that of a laser, involves many photons. Therefore, directly representing wave states by qubit states is inefficient. Second, parameters of qubits may not be easily adjustable [22] to generate Eq. (1). Since quantum hardware likely contains additional terms in $H_{0}$ that cannot be completely turned off, direct emulation is demanding of the hardware tunability and may not be realizable for problems of practical interest.

In comparison, simulation that encodes the interaction at software level with reusable modular gates is more versatile and allows quantum processors to operate as general-purpose computing machines, which do not need to be rebuilt to solve a different problem. In what follows, we construct an efficient algorithm in the action space to enact the nonnative cubic Hamiltonian evolution. Here, the term “action” originates from the action-angle variables of the slow-envelope approximation, which is commonly used to derive the three-wave equations. For Eq. (1), two independent action operators that commute with $H$ are

where $n_{j}=A_{j}^{\dagger}A_{j}$ is the number operator. The simultaneous eigenspace of $H$, $S_{2}$ and $S_{3}$ has dimension $D=\min(s_{2},s_{3})+1$, where the integer $s_{j}\geq 0$ is the eigenvalue of $S_{j}$. Without loss of generality, suppose $s_{2}\leq s_{3}$, which breaks the $2\leftrightarrow 3$ symmetry. Then, in the Fock basis $|n_{1},n_{2},n_{3}\rangle$, any state in the $D$-dimensional subspace $V$ is spanned by

where the expansion coefficients satisfy the normalization condition $\sum_{j}|\alpha_{j}|^{2}=1$. It is easy to check that $S_{j}|\psi\rangle=s_{j}|\psi\rangle$. Moreover, it is important to recognize that $V$ is a closed subspace under the action of $H$, namely, $H|\psi\rangle\in V$ whenever $|\psi\rangle\in V$. Therefore, the total Hilbert space can be decomposed into a direct sum of invariant subspaces, and the Hamiltonian dynamics is confined within the subspaces specified by initial conditions. Notice that using $D=s_{2}+1$ levels, the above mapping can efficiently represent $s_{3}\gg s_{2}$ photons.

Moreover, using the action-space algorithm, the nonlinear three-wave problem is mapped to a finite-dimensional linear Hamiltonian simulation problem. In the Schrödinger picture, $i\partial_{t}|\psi\rangle=H|\psi\rangle$ becomes a system of equations for $\alpha_{j}(t)$

where $h_{j-\frac{1}{2}}=\sqrt{j(s_{2}+1-j)(s_{3}-s_{2}+j)}$ with $h_{-\frac{1}{2}}=h_{D+\frac{1}{2}}=0$. In other words, $H$ is represented by a block tridiagonal matrix with zero diagonal elements. The above mapping utilizes the special form of $H$, but similar algorithms may be constructed for other nonlinear problems when the Hamiltonian has special symmetries.

Once the expansion coefficients are solved for given initial conditions, observables of interest are evaluated by $O(D)$ classical operations during post processing. For example, occupation numbers of the three waves are given by $\langle n_{1}\rangle=\sum_{j=0}^{s_{2}}(s_{2}-j)|\alpha_{j}|^{2}$, $\langle n_{2}\rangle=\sum_{j=0}^{s_{2}}(s_{3}-s_{2}+j)|\alpha_{j}|^{2}$, and $\langle n_{3}\rangle=\sum_{j=0}^{s_{2}}j|\alpha_{j}|^{2}$. Higher-order cumulants, thereby all possible expectation values of interest, can be obtained similarly, and the quantum three-wave problem is then solved.

Before solving the Schrödinger equation, it is instructive to also analyze the cubic interactions in the Heisenberg picture. The number operators satisfy

In comparison, the classical wave action $I_{j}=|A_{j}|^{2}$, where $A_{j}$ is treated as a complex number, satisfies

When identifying $I_{j}\simeq\langle n_{j}\rangle$, the difference between quantum and classical systems is proportional to $3(\langle n_{1}^{2}\rangle-\langle n_{1}\rangle^{2})-\langle n_{1}\rangle$. The term in the parenthesis is exactly zero for states with Poisson statistics, and is usually small compared to other terms when the state is a well-localized semi-classical state. The second term is comparatively small when $s_{2}$ and $s_{3}$ are large, whereby stimulated processes, which depend quadratically on the number of photons, dominate spontaneous emission, which depends linearly on the number of photons. When both terms are small, the system is in the semiclassical regime [23], where the quantum and classical solutions are similar. In the opposite regime, the behaviors of the two systems are disparate. For example, in the thermodynamic limit where occupation numbers satisfy the Boltzmann distribution $|\alpha_{j}|^{2}=1/D$, the quantum system is stationary, whereas the classical system, which is solved by Jacobi elliptic functions [24, 25], is nonstationary.

In the following, we demonstrate the algorithm using three levels, or equivalently, the $|00\rangle$, $|01\rangle$, and $|10\rangle$ states of two qubits. The allowable dimension is $D\leq 3$, and we take $s_{2}=2$ and $s_{3}=s\geq 2$. Then, a natural mapping is $|2,s-2,0\rangle=(1,0,0)^{\text{T}}$, $|1,s-1,1\rangle=(0,1,0)^{\text{T}}$, and $|0,s,2\rangle=(0,0,1)^{\text{T}}$. The normalized Hamiltonian $h=H/|g|$, when restricted to the invariant subspace, is

where $\exp(i\theta)=ig/|g|$. It is worth emphasizing that the general case is associated with a $D\times D$ matrix, where $D\simeq 2^{n}$ with $n$ being the number of qubits needed. In other words, $D$ needs not be three, but the Hamiltonian matrix is always tridiagonal. For the three-level problem, Eq. (11) is exponentiated analytically ¹¹1See Supplemental Material for the unitary matrix of the three-level problem. to determine the unitary time-evolution operator $U(\Delta\tau,\theta,s)=\exp[-ih(\theta,s)\Delta\tau]$. This unitary matrix, for a given $\Delta\tau=|g|\Delta t$, is an input for quantum compilers. Notice that the quantum hardware may not be able to track $U(\tau)$ continuously, because the cubic Hamiltonian has odd parity whereas the native Hamiltonian may not be. Nevertheless, at discrete time steps, $U(\Delta\tau)$ can be realized just like other complex unitary operators.

To compile the three-level unitary operator into a hardware executable form, the standard approach approximates it using a sequence of universal gates. As an example of state-of-the-art quantum cloud services, we implement the unitary operator on Aspen-4-2Q-A of Rigetti Computing [27, 28]. The probabilistic Quil compiler converts $U(\Delta\tau,\theta,s)$ to a sequence of $\sim 20$ native gates, including single-qubit Pauli gates and two CZ gates ²²2See Supplemental Material for an example the Quil program..

As a more hardware-efficient compilation approach, we realize the unitary operator as a single gate using optimal control: Instead of preparing control pulses for standard gates and then using these gates to approximate a targeted operation, we directly optimize a control pulse for the targeted operation. The device we use is a transmon placed inside a 3D superconducting aluminum microwave cavity, whose details are described in [30]. In addition to enhancing the coherence time, this many-qudit single-readout architecture reduces the number of control lines required to access larger computational spaces, and allows simultaneous control of multiple subsystems without relying on their crosstalk [31, 32, 33, 34, 35, 36]. The control Hamiltonian is [37]

where the complex-valued $\epsilon_{j}(t)$ is the slowly-varying envelope of a microwave field whose carrier frequency is $\Omega_{j}$, and $c_{j}$ is its control operator ³³3See Supplemental Material for an explicit expression of the control operator.. The time dependence of $\epsilon_{j}(t)$ may be engineered such that

where $\mathcal{T}$ is the time-ordering operator, $T$ is the length of the control pulse, and $H_{0}$ is the hardware Hamiltonian. The control is not unique but can be implemented using direct digital synthesis [39, 40] and designed using numerical optimization [41, 42, 43, 44] for small problems. For larger problems, the tridiagonal Hamiltonian may be decomposed as $H=\sum_{k=1}^{m}H_{k}$, where each $H_{k}$ is block diagonal. Each block in $H_{k}$, which is small enough and decoupled from other blocks at the hardware level, can now be enacted using numerically optimized control pulse prepared for the subsystem. Subsequently, the full Hamiltonian evolution due to $H$ can be approximated using the Lie-Trotter-Suzuki expansion.

As a test case, we construct a control pulse to achieve $U(0.1,\pi/2,2)$. Based on a thorough characterization that provides parameters of the native Hamiltonian [30], the control pulse is generated in the $|0\rangle\leftrightarrow|1\rangle$ rotating frame using optimize_pulse_unitary in QuTIP [45, 46], which calls an optimizer that uses the gradient ascent algorithms to iteratively update the control pulse to minimize the fidelity error. The optimization uses zero as the initial guess, and is constrained by $|\epsilon|<0.03$. Two forbidden levels beyond the $D=3$ levels are included in the optimization to suppress possible state leakages at the end of the control pulse. The pulse duration ($T=150$ ns) is much shorter than the coherence time but long enough to allow the lowest three levels of the transmon, which has an anharmonicity of $\sim 0.23$ GHz, to be separately addressed. The control pulse $\epsilon(t)$ is shown in Fig. 1(a), and the absolute value of its Fourier transform $\hat{\epsilon}(f)$ is shown in Fig. 1(b). The optimization takes $\sim 10$ iterations to reach convergence with an error target of $10^{-5}$.

When applying to the hardware, optimized control pulses are (i) concatenated in the rotating frame, (ii) transformed to the lab frame by mixing with a carrier wave at the qudit $|0\rangle\leftrightarrow|1\rangle$ transition frequency $f_{01}\approx 4.1$ GHz, and (iii) compensated for spectral filtering effects of the hardware [30]. The final waveform is synthesized using an arbitrary waveform generator at 32-GHz sampling rate, and is sent to the qudit inside a dilution refrigerator through a series of attenuators and a band-pass filter. Occupations of the three levels are measured using dispersive readout [31, 47, 30]. The measured in-phase and quadrature signals are used to classify the states of the qudit, and the classification errors are partially corrected using a confusion matrix [30], which may cause unphysical occupations slightly above 1 or below 0. The estimated occupations during the pulse application are shown in Fig. 1(c) when the qudit is initialized in the ground state, which has minimal state preparation error. By taking into account hardware-specific decay and dephasing ⁴⁴4See Supplemental Material for the noise model based on the Lindblad master equation.. the measurement results (dots) are well explained by numerical solutions (lines) of the Lindblad master equation [49, 50, 51, 52]. At the end of the control pulse, the intended unitary operator is realized, and the targeted transition probabilities (arrows) are attained. Using a modified process tomography [30], we estimate the process matrix of the cubic gate, which yields an average gate fidelity of 99.3%.

To compare the performance of the two compilation approaches, we repeatedly apply the unitary operator to simulate Hamiltonian time evolution. With the quantum processors initialized in the ground state, we read out occupations after repeating $U(0.2,\pi/2,2)$ for $N$ times, and compare the results with the exact solutions (Fig. 2, orange). Using a modular cubic gate ($T=80$ ns), the experimental results (blue) follow the exact solutions up to $N\gtrsim 100$, and the deviations are explained by the master equation solutions (black) using our noise model [48]. In comparison, using a sequence of standard gates, the experimental results (cyan) track the exact solutions only up to $N\sim 10$, which fails to capture even a single period of the nonlinear oscillation. However, it is worth emphasizing that the hardware performance is comparable, and both quantum processors can perform $\sim 100$ gates with high fidelity. The difference in results mostly comes from how the unitary operator is realized: Using the standard compilation approach, each $U$ is realized by $\sim 20$ native gates, whereas using the modular compilation approach, $U$ is realized as a single gate. By reducing the gate depth, the achievable simulation depth is thus significantly increased, which in this case qualitatively changes whether a useful Hamiltonian simulation is feasible or not.

While offering better performance, modular gates are not more costly to compile on classical computers. The standard compilation approach first prepares a universal set of gates that are executable on quantum hardware, and then approximates a targeted $D\times D$ unitary operator $U$ using these native gates. Apart from the cost of preparing the gate set, compiling an arbitrary $U$ requires $O(D^{2})=O(2^{2n})$ native gates [53]. Exponential speedup is possible only when $U=\exp(-iHt)$ and $H$ is native to the hardware. When $H$ is not native, sophisticated quantum algorithms [54, 55, 56, 57, 58] still have favorable complexity scaling, but are not readily implementable on hardware. In comparison, directly compiling $U$ as a modular gate using optimal control requires $O(D^{2})$ classical steps as well. However, once the control pulse is obtained, the unitary operator is enacted much more efficiently as a single pulse. Then, compared with classical computers, which always take $O(D^{2})$ operations to compute $U\psi$, quantum processors may be used as special-purpose hardware accelerators for matrix multiplications. Moreover, the overhead of modular compilation may be amortized when the same $U$ is needed repeatedly or when a family of $U$ can be compiled without expensive numerical optimization.

In fact, high-fidelity parametric gates may be compiled cheaply using interpolation. For example, consider a one-parameter family of the cubic Hamiltonian $h(s)=\sqrt{2s}[(1-\xi)K(2)+(\xi-1/\sqrt{2})K(\infty)]/(1-1/\sqrt{2})$, where $\xi(s)=\sqrt{1-1/s}$, and the nonzero elements of the symmetric $3\times 3$ matrix $K$ are $K_{12}(s)=\xi$ and $K_{23}(s)=1$. Motivated by this form, the interpolated control pulse for $s\geq 2$ at each time slice is taken to be $\epsilon_{\text{I}}(s)=[(1-\xi)\epsilon_{\text{O}}(2)+(\xi-1/\sqrt{2})\epsilon_{\text{O}}(\infty)]/(1-1/\sqrt{2})$, where $\epsilon_{\text{O}}(s)$ denotes the optimized pulse and $\tau\sqrt{2s}=0.2$ is held a constant. In other words, we interpolate the control pulse using the same formula for interpolating the Hamiltonian. To test this scheme, we numerically solve the master equation to obtain the gate fidelity, which is defined by [53] $F(\rho,\sigma)=\text{tr}\sqrt{\rho^{1/2}\sigma\rho^{1/2}}$, where $\rho$ and $\sigma$ are the density matrices attained using optimized and interpolated pulses after one gate application. To optimize control pulses for $\rho(s)$, $\epsilon_{\text{O}}(2)$ is used as the initial guess. Although this compilation shortcut is not proven to work in general, the interpolated pulses are able to drive high-fidelity results for the three-level problem (Fig. 3).

In summary, we demonstrate that nonnative cubic interactions can be programmed entirely at the software level using an efficient action-space algorithm. Moreover, modular compilation reduces the requisite gate depth and qualitatively improves performance on noisy hardware. The resultant high-fidelity gates provide necessary building blocks for Hamiltonian simulations of a large class of nonlinear problems.