Numerical gate synthesis for quantum heuristics on bosonic quantum processors

1 Introduction

Quantum information processing arbitrates controlled interaction of the Hilbert space of a quantum system, for the purpose of generating a target probability distribution expressed in a computational basis defined by an experimental measurement scheme. The Hilbert space of a quantum system generally grows exponentially with the number of degrees of freedom, but for the purpose of quantum information processing it needs to be opportunistically partitioned in order to execute algorithms. The most common encoding exploits a collection of qubits, two-level systems, and it is known that the dimensionality of the Hilbert space is maximized when the states are arranged as a collection of qutrits, three-level systems, for a fixed number of allowed quantum states [1]. Without loss of generality, multiple qudits can be merged into the definition of a new qudit, and a qudit can be mapped via binary encodings into a minimum of log₂(d) qubits and viceversa. For instance, the binary expansion

$\underset{k=0}{\overset{N-1}{\otimes }}|s_k\rangle \underset{\mathit{qubits}}{\overset{\mathit{qudits}}{\rightleftharpoons }}\left|\sum _{k=0}^{N-1}{s}_k{d}^k\right.⟩(1)$

would map the computational basis state of a ququart (i.e. a four-level qudit) onto two qubits |0⟩ → |00⟩, |1⟩ → |01⟩, |2⟩ → |10⟩, |3⟩ → |11⟩.

It is known since the beginning of quantum computing architecture research that universal quantum computing could be achieved by operating constructively on single-qubit and two-qubit at a time [2] via the implementation of quantum gates temporally arranged into quantum circuits. A similar result is known for qudits of arbitrary dimension [3, 4], which can provide hardware-efficient solutions [5] and lower-depth gate compilation and noise improvement compared to qubit-based systems [6–10]. Of particular interest in the current period of technological maturity of quantum processors (the NISQ Era [11]) are variational algorithms such as the Quantum Approximate Optimization Algorithm (QAOA) and the Variational Quantum Eigensolver (VQE) that might achieve some quantum advantage without the fault-tolerance overhead of active error-correction [12]. Typically the quantum circuits of these algorithms feature unitary gates implementing a set of parametrized single-qudit rotations U_M(β) depending on some real angle β. For instance, let us consider the set of SU(2) rotations around the X-axis of the Bloch sphere for qubit systems, and the set of SO(3) rotations that leave invariant the |0⟩ + |1⟩ + |2⟩ state for qutrits. Their matrix representations $U_M^{(2)}(\beta )$ and $U_M^{(3)}(\beta )$ are, respectively:

$U_M^{\left(2\right)}\equiv \left[\begin{array}{cc}\hfill c_{\frac{\beta }{2}}\hfill & \hfill -i{s}_{\frac{\beta }{2}}\hfill \\ \hfill i{s}_{\frac{\beta }{2}}\hfill & \hfill c_{\frac{\beta }{2}}\hfill \end{array}\right](2)$

$U_M^{\left(3\right)}\equiv \frac{1}{3}\left[\begin{array}{ccc}\hfill 1+2c_{\beta }\hfill & \hfill 1-c_{\beta }-\sqrt{3}{s}_{\beta }\hfill & \hfill 1-c_{\beta }+\sqrt{3}{s}_{\beta }\hfill \\ \hfill 1-c_{\beta }+\sqrt{3}{s}_{\beta }\hfill & \hfill 1+2c_{\beta }\hfill & \hfill 1-c_{\beta }-\sqrt{3}{c}_{\beta }\hfill \\ \hfill 1-c_{\beta }-\sqrt{3}{s}_{\beta }\hfill & \hfill 1-c_{\beta }+\sqrt{3}{s}_{\beta }\hfill & \hfill 1+2c_{\beta },\hfill \end{array}\right]$

Where c_x, s_x indicate cos(x) and sin(x) and the computational basis states are ordered in the canonical ascending way. The two-qudit gates of interest for QAOA/VQE ansätze are often diagonal in the computational basis. For instance, the following two-qudit and two-qutrit unitary gates U_C(γ) introduce a phase shift by the angle γ if the two qudits have the same computational state:

$U_C^{\left(2\right)}\equiv \left[\begin{array}{cccc}\hfill e^{i\gamma }\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill e^{i\gamma }\hfill \end{array}\right](3)$

$U_C^{\left(3\right)}\equiv \left[\begin{array}{ccccccccc}\hfill e^{i\gamma }\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill e^{i\gamma }\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill e^{i\gamma }\hfill \end{array}\right],$

where the canonically ordered basis for the matrix representation is used¹. Note that for the most common case of qubits $U_C^{(2)}\propto \exp (i(\gamma /2){\sigma }_z\otimes {\sigma }_z)$, where σ_z are the standard Pauli matrices. For circuit quantum electrodynamics (cQED) systems, note also that there are ways to find effective spin models, which is generally used for encoding of quantum heuristic algorithms [13].

Implementing parametrized gates such as (2–3) starting from the elementary interactions provided by a NISQ processor is a non-trivial problem of synthesis [14], which often can be tackled only via heuristic numerical approaches and online experimental calibration [15]. In this work, we consider the problem of synthesis of gates of the type (2–3) by driving with carefully optimized time-dependent interactions in a system of interacting states. More specifically, the Hilbert space we are considering is spanned by a truncated set of anharmonic bosonic modes, defined with second quantized operators a_m, coupled in a density-density fashion. The corresponding many-body Hamiltonians and their truncated diagonal first quantization representations are:

$H_m={\omega }_m{a}_m^{†}{a}_m+{\xi }_m{\left(a_m^{†}{a}_m\right)}^2(4)$

$\begin{array}{ccc}\hfill {\left.\right|}_{n_m}& \to \hfill & |0\rangle \langle 0|+\sum _{n=1}^{n_m-1}\left[{\omega }_{m}n+{\xi }_m{n}^2\right]|n\rangle \langle n|,\hfill \\ \hfill H_{mm\prime }^{\mathit{int}}& =\hfill & {\xi }_{m{m}^{\prime }}{a}_m^{†}{a}_m{a}_{m^{\prime }}^{†}{a}_{m\prime }\hfill \end{array}(5)$

${\left.\right|}_{n_m{n}_{m^{\prime }}}\to |00\rangle \langle 00|+\sum _{n=1}^{n_m-1}\sum _{k=1}^{n_{m\prime }-1}{\xi }_{m{m}^{\prime }}nk|nk\rangle \langle nk|,$

Where n_m, n_m′ are the number of levels considered for each mode. In photonic implementations, ξ_m is called the self-Kerr coefficient for mode m and ξ_mm′ is called the cross-Kerr coefficient between modes m and m′.We are considering two illustrative setups in order to describe how quantum information could be manipulated in systems featuring Hamiltonians of the type (4–5). In particular we consider that there is one “control qubit mode” T whose Hilbert space is truncated to the first two computational states:

$H_T=|{0\rangle }_T\langle 0{|}_T+\left({\omega }_T+{\xi }_T\right)|{1\rangle }_T\langle 1{|}_T,(6)$

and either one or two computational modes (C) interacting with the control mode, respectively truncated to the first 8 and 3 computational states:

${\mathcal{H}}^{\left(A\right)}=H_T{\left.+H_m+H_{Tm}^{\mathit{int}}\right|}_{\begin{array}{c}{n}_T=2\\ n_m=8\end{array}}(7)$

${\mathcal{H}}^{\left(B\right)}={\left.\begin{array}{c}{H}_T+H_l+H_m\hfill \\ +H_{Tl}^{\mathit{int}}+H_{Tm}^{\mathit{int}}\hfill \end{array}\right|}_{\begin{array}{c}{n}_T=2\\ n_l=3\\ n_m=3\end{array}},$

Where the dependence over the ω and ξ parameters of the Hamiltonians is implied. This setup is a specific case of a generalized Jaynes-Cummings model [16]. Note that for ${\mathcal{H}}^{(B)}$, each C mode is naturally a qutrit, while as noted in Eq. 1, the quantum occupation numbers of the cavity modes could be directly associated to qubit registers via binary expansion [17].

In Figure 1, we show the energy spectrum of two specifications of ${\mathcal{H}}^{(A)}$ and ${\mathcal{H}}^{(B)}$ as well as a pictorial representation of a possible experimental setup that could be described by such effective Hamiltonians: A transmon circuit is embedded into a multimode 3D superconducting cavity, driven by the field of a coupled antenna. Indeed, our reference Hamiltonians can be derived by considering the superconducting transmon to be coupled with the cavity resonator in a dispersive way, i.e., by considering the effective interaction derived by perturbation theory assuming that the ratio of the transmon-cavity coupling and the difference between the transmon and the cavity fundamental frequencies is small, and neglecting the small effective couplings (i.e. cross-Kerr) between the cavity modes [16, 18]. The quantum control drive can be introduced in the model by adding a time-dependent term that allows to create and destroy excitations of a mode m:

$H_m^{\mathit{drive}}\left(t\right)=d_m\left(t\right)a_m+{\bar{d}}_m\left(t\right)a_m^{†},(8)$

where d_m(t) are complex functions. This control Hamiltonian could be related to the (comparatively slowly varying) field generated by the antenna via phenomenologically justified approximations [19].

FIGURE 1

FIGURE 1. Top: System figure with colored waves representing different cavity electromagnetic modes. Eigenspectrum for ${\mathcal{H}}^{(A)}$ (left) and ${\mathcal{H}}^{(B)}$ (right). System parameters and resonant frequencies are given in Section 2. Arrows indicate the transition frequencies. Dashed (continuous) arrows represent transitions between different energy levels with |0⟩_T and |0⟩_T (|1⟩_T).

Having introduced the main definitions and the systems under study, we outline the rest of the paper. In Section 2 we present the synthesis problem from a numerical point of view, following the implementation of quantum optimal control numerics in the open source package Juqbox.jl [20]. Subsection IIA will present results for the synthesis of simple QAOA proof-of-concept circuits based on the parallel execution of gates (2)–(3). Finally, in Section 3 we will discuss future work, including improvements and generalizations of our case study to larger and realistic systems and what is needed for this method to be applied in practice for compilation of variational quantum algorithms in realistic bosonic quantum processors based on 3D cQED technology.

2 Pulse engineering approach

The gate synthesis problem that we are facing could be framed as the task of discovering the functions d_m(t) that allow the Schrödinger evolution for a time τ of $\mathcal{H}+H^{\mathit{drive}}$ to match as close as possible a target unitary operation U:

$U\simeq \mathcal{U}\left(\tau \right)=\mathcal{T}\exp \left[-\frac{i}{\hslash }{\int }_0^{\tau }dt\left(\mathcal{H}+H^{\mathit{drive}}\left(t\right)\right)\right](9)$

In particular, as discussed in the previous section, we will be considering Eq. 7 for $\mathcal{H}$ and Eqs (2-3) as target unitary matrices. In order to solve synthesis numerically, the problem is cast into an optimization challenge over a finite number of real parameters, which can be tackled following the theory of quantum optimal control (QOC) [21]. There are multiple strategies currently implemented for gate sythesis via QOC or machine learning, all with respective benefits and tradeoffs. However, these methods are currently tested on specific limited cases, and insights are difficult to generalize, e.g. see [22–24]. In this paper we follow the techniques described in Ref. [25], targeting specifically cQED models, which we will now briefly review and contextualize for the system under study.

We leverage a key simplification of the QOC problem, consisting in the decomposition of the d_m(t) control functions into a truncated basis spanned by a linear combination of N_b B-spline quadratic polynomials, S_b(t), corresponding to wavelets modulated with N_f resonant frequencies, i.e.

$\begin{array}{ccc}\hfill d_m\left(t\right)& =\hfill & \sum _k^{N_f}{e}^{i{\mathrm{\Omega }}_{m,k}}{W}_{m,k}\left(t\right)\hfill \\ \hfill W_{m,k}\left(t\right)& =\hfill & \sum _b^{N_b}{\alpha }_{m,k,b}{S}_b\left(t\right),\hfill \end{array}(10)$

where αs are complex coefficients, representing the unknowns of the optimization problem. The choice of B-splines as a basis for expansion is motivated by computational efficiency of parametrization of the control functions. The resonant frequencies Ω_m,k are defined by considering the energy differences between the states corresponding to the creation or annihilation of a boson, leaving the remaining occupations unchanged. Signals tuned at these frequencies initiate transitions as it can be proven by first order time-dependent perturbation theory.

We show in Figure 1 the resonant frequencies for our illustrative systems: for ${\mathcal{H}}^{(A)}$, we count 8 transitions related to T-bosons and 14 transitions for C-bosons for a total of 22 resonant frequencies. For ${\mathcal{H}}^{(B)}$, there are 9 resonant frequencies in total that trigger T transitions, and 24 transitions related to the C modes. However, some transitions are degenerate–only 17 different frequencies are required.

We consider the following values of parameters, with reference to a perspective reference cQED potential implementation: ω_T/2π = 5 GHz; ω_m/2π = 3 GHz, ω_l/2π = 4 GHz; ξ_m/2π = 0.6 MHz, ξ_l/2π = 0.9 MHz; ξ_T/2π = 200 MHz. In line with our inspiration of a cavity-transmon systems in the dispersive regime [26], we assign interaction parameters to be the geometric means of the local self-interactions ${\xi }_{Tm}/2\pi =\sqrt{{\xi }_m\times {\xi }_T}/2\pi =$ 10.95 MHz and ${\xi }_{Tl}/2\pi =\sqrt{{\xi }_l\times {\xi }_T}/2\pi =$ 13.42 MHz. The parameters that we used for ${\mathcal{H}}^{(A)}$ and ${\mathcal{H}}^{(B)}$ are inspired from expectations of results that would be obtained applying black-box quantization to Tesla-cavity systems [27] coupled dispersively to transmons with coherence times ≃ 100μs [28]. Following that inspiration, we assume small linewidth for the cavity mode compared to their separations and we set the minimum frequency difference between the transmon and cavity mode frequencies to be of the order of the GHz, in order to justify independent access of the control pulses to the transmon and for each cavity modes.

Following Juqbox.jl [20], the pulse engineering algorithm attempts to discover the best α_m,k,b coefficients (i.e. 2 × N_f × N_b real parameters), which works as follows. Initially, a random pulse is selected by initializing the vector of parameters using random positive numbers uniformly distributed within [0, 0.2 MHz). Then, an objective function is calculated (see Subsection IIA) and the pulse is iteratively updated by computing the Schrödinger evolution and gradients efficiently by symplectic time-integration of adjoint equations (29). Note that due to the B-spline parametrization, the number of control parameters does not depend directly on the pulse total duration τ. However, the number of B-splines N_b defines the design of the temporal structure of the pulses, so one needs to choose large enough τ and N_b to allow the method to converge to a numerically robust solution. In particular, the slowest frequency resolution of the pulses is given by 1/τ. We choose to vary τ in the 500–8000 ns range for our numerical experiments on ${\mathcal{H}}^{(A)}$ and ${\mathcal{H}}^{(B)}$, allowing for a frequency resolution of 0.125–2 MHz. The B-splines vary on the time scale τ/N_b. Hence we choose N_b = 10 to allow resolution at the scale of ξ_Tm, which controls multiple energy separations in the spectrum. The values of ξ_l, ξ_m define the smallest resonant frequencies.

2.1 Evaluation test case: QAOA

Our numerical prototype experiment is based on the synthesis of QAOA-like quantum circuits, which in their basic implementation consist of the layered alternated application of phase-separation unitary gates and mixing gates [30]. With reference to the known Max-k-Cut qudit mapping of QAOA [31], where k corresponds to the dimensionality of the qudits, we can craft the phase-separation layers using U_C(γ)_ij gates and we have the freedom of designing the mixing layers using the U_M(β)_i gates in Eqs 2-3, where i, j indicate the distinguishable qudits that are targeted by the specific gate execution. Other choices would also be appropriate [32]. For clarity, in Figure 2A, we show the two toy-model circuits that we are going to synthesize, respectively via pulse engineering on ${\mathcal{H}}^{(A)}$ and ${\mathcal{H}}^{(B)}$. For completeness, the test circuit include the initialization operation, which is usually taken to be a generalized Hadamard gate (although it could be substituted by a mixing over the |0⟩^⊗N state).

FIGURE 2

FIGURE 2. (A) Prototype circuits for the synthesis of Max-k-Cut QAOA algorithms. Single C-mode represent 8 computational states (equivalent to 3 qubits). (B) Illustrative Fourier spectrum of a high-fidelity engineered pulse via Juqbox.jl. Top row shows results for d_T(t) while the bottom row shows the control of the computational modes (d_m(t) and d_l(t)). Darker tones (black, blue, orange) indicate the pulses that synthesize mixing layers, while light tones (gray, cyan, yellow) refer to phase-separation layers. (C) Fidelity for pulse engineered QAOA layers of the prototype circuits. Black lines indicate the mean across angles, individually plotted in gray. Each line is the mean of 10 random restarts (20–80 percentiles across restarts is plotted as shaded area). Leakage plots are presented in the Supplementary Material.

We note that quantum processor programmers have formally the freedom to execute gates sequentially or in parallel, and to exchange them in temporal execution order if they commute. However, in a real world implementation, if the processor is not fault-tolerant, under reasonable assumptions we expect decoherence and dephasing errors to be roughly proportional to execution time, so a compiler for NISQ algorithms often tries to parallelize gate execution as much as possible [33]. Moreover, considering the mapping of the computational variables to the spectrum of the Hamiltonians (Figure 1), the possible qudit identity assignments are inequivalent with respect to pulse engineering, although it would be inconsequential if the synthesis was perfect. SWAP operations could restrict the number of active qudits, by relegating some states to be just memory storage and not participate in processing. However, these operations and controls for our Hamiltonians need to be synthesized as well, increasing the complexity of the entire compilation significantly. Bearing in mind these considerations, in our case study we choose to implement the single-qudit gates in parallel when possible, without implementing SWAPs but directly synthesizing all required two-body interactions instead across the entire Hilbert space. We will discuss in Section 3 the scalability issues associated to this approach.

Noting that in cQED implementations, the Hamiltonians in Eq. 7 are defined on truncated versions of a physically infinite Hilbert space, it is customary to include a few additional guard states corresponding to high occupation of boson modes to help the robustness of the numerical optimization, i.e., the following parameters are renormalized $n_T\to {\tilde{n}}_T=n_T+\delta n_T$, $n_m\to {\tilde{n}}_m=n_m+\delta n_m$ and $n_l\to {\tilde{n}}_l=n_l+\delta n_l$, where δn represent guard states with values in Table 1.

TABLE 1

TABLE 1. Parameters used for prototype (See Figure 2).

Following [25], the optimization objective to be minimized is chosen to be a sum of the infidelity and average leakage. The infidelity is a measure of a similarity score between the synthesized unitary matrix and the target, which can be defined as $O_F=1-|\mathrm{T}\mathrm{r}(\mathcal{U}{(\tau )}^{†}U)/E{|}^2$, where E is a normalization constant. The average leakage is defined as $O_L=(1/\tau ){\int }_0^{\tau }\mathrm{T}\mathrm{r}({\mathcal{U}}^{†}(t)W\mathcal{U}(t))dt$, where W is a diagonal matrix which is non-zero only on the indices corresponding to the guard levels. The weights in W are set to be 1.0 for the highest guard state and then decrease exponentially in powers of 10 for each lower state. The objective of the numerics is to minimize O = O_F + O_L by solving the related optimization problem on the α parameters by using the IPOPT L-BFGS optimizer [34] and using the efficient Juqbox.jl numerical integration scheme to compute the required O and ∇_αO.

The optimization heuristics has a stopping condition based on either the achievement of a target threshold fidelity (1-O_F) or the execution of a maximum number of iterations. As mentioned, we perform multiple restarts initializing the optimization with different random pulses (see Table 1 for a summary of some of the parameters used for the numerical experiments). Computations have been performed allowing an optimization time in the order of days. See Supplemental Material for computational details.

To give a sense of the resulting control signals that generate the QAOA circuit layers, we show the resulting Fourier transform of one engineered d_T(t), d_m(t), d_l(t) functions in Figures 2B, for one random seed and pulse time τ = 8000 ns, which in retrospect we know guaranteeing high fidelity of the synthesis. The angle parameters β and γ have been set to a fixed arbitrary value of π/5 for illustration but the qualitative features of the pulses that we are describing are preserved for different τ and angles. As evident from the plots, the scheme and parameters described above clearly generate peaks around the identified resonant frequencies corresponding to the single-boson transitions in Figure 1. In particular, for the d_T(t) controls, the highest peak corresponds to ω_T while the other equispaced peaks are centered among multiples of ξ_Tm, ξ_Tl or integer combinations of the two energy values for the ${\mathcal{H}}^{(B)}$ system. For the C-mode controls, the peaking frequencies are topped by ω_m, ω_l and generally ξ_Tm, ξ_Tl plus multiples of ξ_m, ξ_l respectively. Clearly, each reported spectrum corresponds to a real-time microwave combination of pulses that can be crafted via an arbitrary waveform generator (AWG) in an experimental setup.

In Figures 2–c, we provide the aggregated performance of the pulse engineering approach, plotting the fidelity between the final pulse $\mathcal{U}(\tau )$ and target circuit layers (phase separation and mixing) and Hadamard gates, for different pulse times. We show the mean fidelity, estimated averaging 10 random initializations, i.e. restart of the L-BFGS optimizer (the default optimizer for Juqbox.jl), for QAOA layers parametrized with 11 different γ and β (from -π to π in fractions of π/5). As expected, notwithstanding outliers, the statistics is sufficient to indicate that the method can reach the target 0.99 fidelity if the pulse is allowed to be sufficiently long.

3 Discussion and outlook

In the previous section, we described a proof-of-concept of numerical synthesis for simple quantum circuits describing the building blocks of Max-k-Cut QAOA algorithms using qubits (mapped onto qudits) and qutrits, on bosonic quantum processors. The main question that is left to be addressed is if the synthesis approach we employed is sufficiently robust to be applied at application-scale. We break down the question in a discussion of three scalability challenges: Computational effort, realistic implementation, and circuit fidelity.

3.1 Computational effort

As mentioned, the computational effort required by numerical packages to obtain high-fidelity in our case study is already very significant and scales both with the Hilbert space size and with the pulse duration. This means that the proposed methodology will most certainly not be viable if straightforwardly applied to systems at large scale, although larger synthesis can be achieved if the code is optimized to leverage GPU clusters. The envisioned practical synthesis of larger circuits will necessarily need to be broken down in modules, each of which working on a subspace of the entire Hilbert space. The requirement for this modularization is that the gate synthesized numerically in a system with few modes will have to be applied in a system with several modes and levels. The optimal gate from numerics should ideally act as an identity on the degree of freedoms that were not considered in the synthesis in order not to cause the crosstalk problem [35]. Scaling up the single-mode case ${\mathcal{H}}^{(A)}$ that we used will not likely be viable, since the non-local mapping onto qubits would require any gate to address the entire level structure independently from the locality of the gates, which is why we opted to synthesize the entire phase-separation circuit as opposed to the individual two-qubit gates independently. However, it is envisionable to generalize the ${\mathcal{H}}^{(B)}$ system adding more C-modes, i.e., considering the Hamiltonian

${\mathcal{H}}_{\mathit{multi}}^{\left(B\right)}\left[N\right]=H_T+\sum _{j=1}^N\left[H_{m_j}+H_{T{m}_j}^{\mathit{int}}\right],(11)$

which is ${\mathcal{H}}^{(B)}$ for N = 2. If the ξ parameters of each C-mode are sufficiently separated, the peaked frequency structure of the engineered pulses suggests that it is possible that none of the peaks in the final pulses would correspond to resonances with single-boson excitations that we don’t want to trigger, which would likely induce very small leakage outside the two-mode target computational space. This needs to be verified theoretically or numerically in future work. Ultimately, frequency crowding will be an issue and more sophisticated numerics or frequency spacing and bandwidth engineering will be required.

It should be noted that if the modularization works as expected, the computing time spent synthesizing algorithmic primitives would be an offline una tantum cost to be paid to populate a lookup table (LUT) that would be accessed at runtime by the perspective user of the quantum solver. Indeed, similarly as in other domains, it is envisioned that the LUT would be computed for a large grid of parameters (angles γ and β in our QAOA example) and then machine learning algorithms would learn and return an interpolation of the engineered pulses if the compiler is called for a parameter that was not pre-computed, or would use nearby known points to initialize a fast optimization round to engineer a new pulse on the fly [36].

3.2 Realistic implementation

While the described technique is generically applicable to any bosonic interacting system, our case study has a specific 3D cQED implementation in mind, as illustrated in the inset of Figure 1.

It should be noted that the general framework that we employed, pulse engineering via QOC, while proven powerful [37] is not the only known approach to achieve universal synthesis of unitary quantum gates defined in the Fock space for these kind of systems. For instance, the use of selective number-dependent arbitrary phase (SNAP) protocol [38, 39] or echoed conditional displacement [40] are strong candidates for the universal control of a single-mode system. Qudits have potential to be affected by noise less so than qubits [7] but working with large photon-number states comes with additional complications in terms of decoherence, which are still theoretically not entirely understood [41].

The multiqudit system (Eq. (11)) could be viable but its practical implementation will likely suffer from the aforementioned quantum and classical crosstalk problems whose handling is currently one of the main active research topics of the 3D multimode cQED domain [42]. Even assuming that the bandwidth of the control pulses and the level spacing has sufficient resolution, there is a need for the co-design of a NISQ cQED architecture that would allow two-mode gates to operate in large Hilbert space with a controllable effect over spectator modes that are subject to an always-on interaction [43]. Theory results on quantum adiabatic protocols [44, 45] on bosonic systems could provide an initial reference point to be generalized [46, 47].

3.3 Fidelity

The fidelity target we used in our prototype (0.99) is in line with the fidelity of native gates in industrial grade quantum processors but it is of course somewhat arbitrary. In accordance with conservative models of uncorrelated errors, we could estimate the final fidelity of the entire circuits in Figures 2–a as the product of the fidelities of each synthesized layer, which means that ultimately the fidelity decreases exponentially with the number of layers. Hence, quantum-volumetric tests [48] would fail rather fast if we were to scale our circuits beyond few variables. However, it should be noted that for quantum optimization algorithms of the variational type, it is not clear if high fidelities are required, considering that the underlying computational principle is preserved for Lindblad evolution [49]. The degree of freedom of parameter setting might contribute to mitigate the misspecification of the gates due to poor synthesis. The non-requirement of exact synthesis is intuitive, since for optimization tasks we are not necessarily trying to reproduce a quantum process but rather to drive the system towards a probability distribution, which might be achievable also with partially coherent systems or in the presence of spurious unknown interactions that give rise to systematic coherent errors. So, as long as the nature of the errors is not specifically adversarial against the optimization tasks, there is still reasonable hope that a low-fidelity circuit could deliver speedup in the NISQ era. An important contribution that we are considering to improve the fidelity would be to generalize the technique of Juqbox.jl to open systems, and fit the experimental noise to solve for a more realistic model. Fortunately, there has already been active development in that direction, including enabling quantum optimal control and pulse-level programming in XACC [50, 51] with QuaC plugin [52], and a recently released open-source package for high-performance optimal control, Quandary [53].

In conclusion, we investigated the application of quantum optimal control techniques to design unitary gates for a class of physical systems that could be programmed to act as qudit-based quantum computers. We used variational algorithms such as QAOA for qubits (mapped onto a single qudit) and qutrits as targets for our case-study. Our current results, similar to other applied quantum computing works for multimode cQED [54], are still limited on small proof-of-concept models, due to limitations in computational effort, realistic implementation and achievable fidelity. While we identified pathways to overcome such limitations, we should note that for the purpose of variational optimization there are multiple recent attempts to employ co-designed digital-analog approaches that are directly related to QOC as optimization algorithms [14, 55, 56], and might not require the burdens of high-fidelity gate synthesis. We envision that our work could also contribute to those innovative methods that have already been delivering promising results.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

ABÖ run the simulations. Both authors analyzed the results and contributed to the manuscript.

Funding

This material is based upon work supported by the U.S. Department of Energy, Office of Science, National Quantum Information Science Research Centers, Superconducting Quantum Materials and Systems Center (SQMS) under contract number DE-AC02-07CH11359.

Acknowledgments

We thank Jens Koch, Srivatsan Chakram, Taeyoon Kim, Joshua Job, Matthew Reagor, Matthew Otten, Keshav Kapoor, Silvia Zorzetti, Sohaib Alam, Doga Kurkcuoglu and the SQMS 3D Algorithms Group and SQMS Codesign Group for discussions and feedback. We thank Adam Lyon, Jim Kowalkowski, Yuri Alexeev and Norm Tubman for their assistance on computing aspects, including support through XSEDE computational Project no. TG-MCA93S030 providing compute time at Bridges-2 of the Pittsburgh Supercomputer Center. ABÖ thanks Gabriel Perdue, Adam and Jim for their guidance during his early career years. We thank Anders Petersson for his support in configuring Juqbox.jl. DV acknowledges support via NASA Academic Mission Service (NNA16BD14C). We gratefully acknowledge the computing resources provided on Bebop, a high-performance computing cluster operated by the Laboratory Computing Resource Center at Argonne National Laboratory.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2022.900612/full#supplementary-material

Footnotes

¹|00⟩, |01⟩, |10⟩, |11⟩ for qubits, and |00⟩, |01⟩, |02⟩, |10⟩, |11⟩, |12⟩, |20⟩, |21⟩, |22⟩ for qutrits.

References

1. Greentree AD, Schirmer S, Green F, Hollenberg LC, Hamilton A, Clark R, et al. Maximizing the hilbert space for a finite number of distinguishable quantum states. Phys Rev Lett (2004) 92:097901. doi:10.1103/physrevlett.92.097901

PubMed Abstract | CrossRef Full Text | Google Scholar

2. DiVincenzo DP. Two-bit gates are universal for quantum computation. Phys Rev A (Coll Park) (1995) 51:1015–22. doi:10.1103/physreva.51.1015

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Bullock SS, O’Leary DP, Brennen GK. Asymptotically optimal quantum circuits ford-level systems. Phys Rev Lett (2005) 94:230502. doi:10.1103/physrevlett.94.230502

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Wang Y, Hu Z, Sanders BC, Kais S. Qudits and high-dimensional quantum computing. Front Phys (2020) 479: 589504. doi:10.3389/fphy.2020.589504

CrossRef Full Text | Google Scholar

5. Liu Y, Sinanan-Singh J, Kearney MT, Mintzer G, Chuang IL. Constructing qudits from infinite-dimensional oscillators by coupling to qubits. Phys Rev A (Coll Park) (2021) 104:032605. doi:10.1103/physreva.104.032605

CrossRef Full Text | Google Scholar

6. Gokhale P, Baker JM, Duckering C, Brown NC, Brown KR, Chong FT. Asymptotic improvements to quantum circuits via qutrits. In: Proceedings of the 46th International Symposium on Computer Architecture; 2019 June 22–26; Phoenix Arizona (2019). p. 554–66.

CrossRef Full Text | Google Scholar

7. Otten M, Kapoor K, Özgüler AB, Holland ET, Kowalkowski JB, Alexeev Y, et al. Impacts of noise and structure on quantum information encoded in a quantum memory. Phys Rev A (Coll Park) (2021) 104:012605. doi:10.1103/physreva.104.012605

CrossRef Full Text | Google Scholar

8. Blok MS, Ramasesh VV, Schuster T, O’Brien K, Kreikebaum J-M, Dahlen D, et al. Quantum information scrambling on a superconducting qutrit processor. Phys Rev X (2021) 11:021010. doi:10.1103/physrevx.11.021010

CrossRef Full Text | Google Scholar

9. Gustafson EJ. Prospects for simulating a qudit-based model of (1+1)D scalar QED. Phys Rev D (2021) 103:114505. doi:10.1103/physrevd.103.114505

CrossRef Full Text | Google Scholar

10. Gustafson E, Noise improvements in quantum simulations of sqed using qutrits, arXiv preprint arXiv:2201.04546 (2022).

Google Scholar

11. Preskill J. Quantum computing in the nisq era and beyond. Quantum (2018) 2:79. doi:10.22331/q-2018-08-06-79

CrossRef Full Text | Google Scholar

12. Cerezo M, Arrasmith A, Babbush R, Benjamin SC, Endo S, Fujii K, et al. Variational quantum algorithms. Nat Rev Phys (2021) 1:625–44. doi:10.1038/s42254-021-00348-9

CrossRef Full Text | Google Scholar

13. Miyazaki R, Effective spin models of kerr-nonlinear parametric oscillators for quantum annealing, arXiv preprint arXiv:2202.07160 (2022).

CrossRef Full Text | Google Scholar

14. Magann AB, Arenz C, Grace MD, Ho T-S, Kosut RL, McClean JR, et al. From pulses to circuits and back again: A quantum optimal control perspective on variational quantum algorithms. PRX Quan (2021) 2:010101. doi:10.1103/prxquantum.2.010101

CrossRef Full Text | Google Scholar

15. Klimov PV, Kelly J, Martinis JM, Neven H, The snake optimizer for learning quantum processor control parameters, arXiv preprint arXiv:2006.04594 (2020).

Google Scholar

16. Blais A, Grimsmo AL, Girvin S, Wallraff A. Circuit quantum electrodynamics. Rev Mod Phys (2021) 93:025005. doi:10.1103/revmodphys.93.025005

CrossRef Full Text | Google Scholar

17. Sawaya NP, Menke T, Kyaw TH, Johri S, Aspuru-Guzik A, Guerreschi GG, et al. Resource-efficient digital quantum simulation of d-level systems for photonic, vibrational, and spin-s Hamiltonians. Npj Quan Inf (2020) 6:49. doi:10.1038/s41534-020-0278-0

CrossRef Full Text | Google Scholar

18. Ma W-L, Puri S, Schoelkopf RJ, Devoret MH, Girvin S, Jiang L, et al. Quantum control of bosonic modes with superconducting circuits. Sci Bull (2021) 66:1789–805. doi:10.1016/j.scib.2021.05.024

CrossRef Full Text | Google Scholar

19. Gerry C, Knight P, Knight PL. Introductory quantum optics. Cambridge, United Kingdom: Cambridge University Press (2005).

Google Scholar

20. Petersson NA, Garcia F, Juqbox.jl, GitHub, Available at: https://github.com/LLNL/Juqbox.jl (Accessed September, 2021)(2021).

Google Scholar

21. Palao JP, Kosloff R. Quantum computing by an optimal control algorithm for unitary transformations. Phys Rev Lett (2002) 89:188301. doi:10.1103/physrevlett.89.188301

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Riaz B, Shuang C, Qamar S. Optimal control methods for quantum gate preparation: A comparative study. Quan Inf Process (2019) 18:100. doi:10.1007/s11128-019-2190-0

CrossRef Full Text | Google Scholar

23. Niu MY, Boixo S, Smelyanskiy VN, Neven H. Universal quantum control through deep reinforcement learning. Npj Quan Inf (2019) 5:33. doi:10.1038/s41534-019-0141-3

CrossRef Full Text | Google Scholar

24. Baum Y, Amico M, Howell S, Hush M, Liuzzi M, Mundada P, et al. Experimental deep reinforcement learning for error-robust gate-set design on a superconducting quantum computer. PRX Quan (2021) 2:040324. doi:10.1103/prxquantum.2.040324

CrossRef Full Text | Google Scholar

25. Petersson NA, Garcia F, Optimal control of closed quantum systems via b-splines with carrier waves, arXiv preprint arXiv:2106.14310 (2021).

Google Scholar

26. Nigg SE, Paik H, Vlastakis B, Kirchmair G, Shankar S, Frunzio L, et al. Black-box superconducting circuit quantization. Phys Rev Lett (2012) 108:240502. doi:10.1103/physrevlett.108.240502

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Romanenko A, Pilipenko R, Zorzetti S, Frolov D, Awida M, Belomestnykh S, Posen S, Grassellino A. Three-Dimensional superconducting resonators at T<20 mK with photon lifetimes up to τ=2 s 20 mk with photon lifetimes up to τ= 2 s. Phys Rev Appl (2020) 13:034032. doi:10.1103/physrevapplied.13.034032$\lessbin $

CrossRef Full Text | Google Scholar

28. Nersisyan A, Poletto S, Alidoust N, Manenti R, Renzas R, Bui C-V, et al. IEEE (2019). p. 31–1.Manufacturing low dissipation superconducting quantum processors2019 IEEE International Electron Devices Meeting (IEDM)7-11 Dec. 2019; San Francisco, CA

CrossRef Full Text | Google Scholar

29. Petersson NA, Garcia FM, Copeland AE, Rydin YL, DuBois JL, Discrete adjoints for accurate numerical optimization with application to quantum control, arXiv preprint arXiv:2001.01013 (2020).

Google Scholar

30. Hadfield S, Wang Z, O’Gorman B, Rieffel EG, Venturelli D, Biswas R, et al. From the quantum approximate optimization algorithm to a quantum alternating operator ansatz. Algorithms (2019) 12:34. doi:10.3390/a12020034

CrossRef Full Text | Google Scholar

31. Fuchs FG, Kolden HØ, Aase NH, Sartor G. Efficient encoding of the weighted MAX $$k$$-CUT on a quantum computer using QAOA. SN Comp Sci (2021) 2:89. doi:10.1007/s42979-020-00437-z

CrossRef Full Text | Google Scholar

32. Deller Y, Schmitt S, Lewenstein M, Lenk S, Federer M, Jendrzejewski F, et al. Quantum approximate optimization algorithm for qudit systems with long-range interactions, arXiv preprint arXiv:2204.00340 (2022).

Google Scholar

33. Venturelli D, Do M, Rieffel E, Frank J. Compiling quantum circuits to realistic hardware architectures using temporal planners. Quan Sci Technol (2018) 3:025004. doi:10.1088/2058-9565/aaa331

CrossRef Full Text | Google Scholar

34. Wächter A, Biegler LT. On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Math Program (2006) 106:25–57. doi:10.1007/s10107-004-0559-y

CrossRef Full Text | Google Scholar

35. Özgüler AB, Job JA. Dynamics of qudit gates and effects of spectator modes on optimal control pulses. arXiv preprint [arXiv:2207.14006] (2022).

Google Scholar

36. Xu D, Özgüler AB, Di Guglielmo G, Tran N, Perdue GN, Carloni L, et al. Neural network accelerator for quantum control. arXiv preprint [arXiv:2208.02645] (2022).

Google Scholar

37. Heeres RW, Reinhold P, Ofek N, Frunzio L, Jiang L, Devoret MH, et al. Implementing a universal gate set on a logical qubit encoded in an oscillator. Nat Commun (2017) 8:94. doi:10.1038/s41467-017-00045-1

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Heeres RW, Vlastakis B, Holland E, Krastanov S, Albert VV, Frunzio L, et al. Cavity state manipulation using photon-number selective phase gates. Phys Rev Lett (2015) 115:137002. doi:10.1103/physrevlett.115.137002

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Fösel T, Krastanov S, Marquardt F, Jiang L, Efficient cavity control with snap gates, arXiv preprint arXiv:2004.14256 (2020).

Google Scholar

40. Eickbusch A, Sivak V, Ding AZ, Elder SS, Jha SR, Venkatraman J, et al. Fast universal control of an oscillator with weak dispersive coupling to a qubit, arXiv preprint arXiv:2111.06414 (2021).

Google Scholar

41. Hanai R, McDonald A, Clerk A, Intrinsic mechanisms for drive-dependent purcell decay in superconducting quantum circuits, arXiv preprint arXiv:2106.05179 (2021).

CrossRef Full Text | Google Scholar

42. Chakram S, Oriani AE, Naik RK, Dixit AV, He K, Agrawal A, et al. Seamless high-q microwave cavities for multimode circuit qed, arXiv preprint arXiv:2010.16382 (2020).

Google Scholar

43. Alam MS, Belomestnykh S, Bornman N, Cancelo G, Chao Y-C, Checchin M, et al. Quantum computing hardware for hep algorithms and sensing, arXiv preprint arXiv:2204.08605 (2022).

Google Scholar

44. Özgüler AB, Joynt R, Vavilov MG Steering random spin systems to speed up the quantum adiabatic algorithm. Phys Rev A (2018) 98(6):062311.

Google Scholar

45. Das A, Chakrabarti BK. Colloquium: Quantum annealing and analog quantum computation. Rev Mod Phys (2008) 80:1061–81. doi:10.1103/revmodphys.80.1061

CrossRef Full Text | Google Scholar

46. Pino M, García-Ripoll JJ. Quantum annealing in spin-boson model: From a perturbative to an ultrastrong mediated coupling. New J Phys (2018) 20:113027. doi:10.1088/1367-2630/aaeeea

CrossRef Full Text | Google Scholar

47. Starchl E, Ritsch H. Unraveling the origin of higher success probabilities in quantum annealing versus semi-classical annealing. J Phys B: Mol Opt Phys (2022) 55:025501. doi:10.1088/1361-6455/ac489a

CrossRef Full Text | Google Scholar

48. Blume-Kohout R, Young KC. A volumetric framework for quantum computer benchmarks. Quantum (2020) 4:362. doi:10.22331/q-2020-11-15-362

CrossRef Full Text | Google Scholar

49. Yang Z-C, Rahmani A, Shabani A, Neven H, Chamon C. Optimizing variational quantum algorithms using pontryagin’s minimum principle. Phys Rev X (2017) 7:021027. doi:10.1103/physrevx.7.021027

CrossRef Full Text | Google Scholar

50. Nguyen T, Santana A, McCaskey A. Extending xacc for quantum optimal control. In: 2020 IEEE International Conference on Quantum Computing and Engineering (QCE); 12-16 Oct. 2020. IEEE (2020). p. 391–401.

CrossRef Full Text | Google Scholar

51. Nguyen T, McCaskey A. Enabling pulse-level programming, compilation, and execution in xacc. IEEE Trans Comput (2021) 71:547–58. doi:10.1109/tc.2021.3057166

CrossRef Full Text | Google Scholar

52. Otten M. Quac: Open quantum systems in C, a time-dependent open quantum systems solver (2017).

Google Scholar

53. Günther S, Petersson NA, Dubois JL, Quandary: An open-source c++ package for high-performance optimal control of open quantum systems (2021), arXiv:2110.10310 [quant-ph] .

Google Scholar

54. Kurkcuoglu DM, Alam MS, Li AC, Macridin A, Perdue GN, Quantum simulation of ϕ4 theories in qudit systems, arXiv preprint arXiv:2108.13357 (2021).

Google Scholar

55. Gokhale P, Ding Y, Propson T, Winkler C, Leung N, Shi Y, et al. Partial compilation of variational algorithms for noisy intermediate-scale quantum machines. In: Proceedings of the 52nd Annual IEEE/ACM International Symposium on Microarchitecture; Columbus OH USAOctober 12 - 16, 2019 (2019). p. 266–78.

CrossRef Full Text | Google Scholar

56. Choquette A, Di Paolo A, Barkoutsos PK, Sénéchal D, Tavernelli I, Blais A, et al. Quantum-optimal-control-inspired ansatz for variational quantum algorithms. Phys Rev Res (2021) 3:023092. doi:10.1103/physrevresearch.3.023092

CrossRef Full Text | Google Scholar